不同降温模式下高加索三叶草的转录组比较分析

牟 丹， 马戈亮， 卡着才让， 李宗仁， 忠克吉， 谢久祥*， 李希来*

(1.青海大学农牧学院， 西宁 青海 810016； 2. 河南蒙古族自治县自然资源局， 黄南 青海 811599)

高加索三叶草(TrifoliumambiguumBieb.)是豆科三叶草属中唯一具有地下根蘖和利用根蘖进行强克隆生长的多年生长寿命草种[1]。该物种原产于高纬度的东欧和西亚的高海拔(2 750 m)地区[2]，当地比较恶劣的生态条件使其具有很强的抗寒、抗旱、耐水淹、耐牧等特性[3]。据报道，在高加索山脉海拔3 200 m的草甸中也能发现这种植物[3]。

青藏高原具有显著的高原大陆性气候，昼夜温差大，年均温低，无霜期短，极端低温、倒春寒等灾害天气时有发生[4-5]。低温对分布在高海拔或较高纬度的多年生植物来说是一种不可避免的环境胁迫[6]，严重影响并制约着高寒地区牧草的生长发育和优良豆科牧草的引种工作[7]。

2015年，本课题组将高加索三叶草从新西兰引种到青海省河南蒙古自治县(河南县)海拔3 500～3 700 m的3处试验地，发现其在野外均能自然越冬，并形成了克隆植株[8-9]。河南县地处青藏高原东部，高原大陆性气候较明显。高加索三叶草在河南县的成功越冬说明其抗寒性很强，足以应对河南县极端低温等恶劣天气。目前，对植物响应低温胁迫的转录组研究多集中在单一降温模式处理，如适温培养后经冷驯化缓慢降温处理或直接低温处理[10-12]，而对同一物种在遭遇不同降温模式的低温胁迫转录组比较研究鲜见报道。本研究以高加索三叶草为材料，利用低温人工气候箱模拟缓慢降温和骤然降温2种降温模式对高加索三叶草植株进行-6℃低温胁迫4 h，比较不同降温模式下的基因表达差异性和关键通路富集情况，并进行抗寒相关基因和转录因子筛选，为高加索三叶草抗寒分子育种奠定理论基础。

1 材料与方法

1.1 试验材料

供试材料为长势一致、生长健康的高加索三叶草植株，种子由内蒙古农业大学王明玖教授提供。

1.2 试验设计

试验于2019年5—8月在青海大学校内牧草试验地和草业科学系实验室进行。将种子播种到含混合基质(草炭∶蛭石∶珍珠岩=3∶1∶1)的花盆中。60 d后将长势一致的幼苗转移到低温人工气候箱(光照12 h、昼夜温度20℃/18℃、光照强度12 000 lx、湿度60%)适应性培养30 d。试验分3组：常温对照组(CK)、缓慢降温组(GC)和骤然降温组(SC)，每组均为3个重复。常温对照组的幼苗始终置于适应性培养下的培养箱中生长，2个降温组除温度外，其余条件均与适应性培养下的条件相同，具体过程为：(1)GC组：由20℃以2℃/h降温速率降到4℃/4℃(昼/夜)处理48 h，然后以2℃/h降温速率降到2℃/2℃(昼/夜)处理44 h，最后再直接降温到-6℃进行低温胁迫4 h，期间有光照；(2)SC组：由20℃直接降温到-6℃进行低温胁迫4 h，期间有光照。从10∶00开始，GC和SC 2个处理组均于-6℃低温胁迫4 h，在14∶00对2个处理组和CK组同时采样，剪取叶片并立即置于液氮速冻，后放入-80℃超低温冰箱中保存备用。

1.3 高通量测序及数据分析

使用Trizol试剂盒(Invitrogen,Carlsbad,CA,USA)从高加索三叶草叶片中提取总RNA。RNA质量在Agilent 2100生物分析仪(Agilent Technologies,Palo Alto,CA,USA)上进行评估，并使用无RNase的琼脂糖凝胶电泳进行检测。使用NEB#7530试剂盒(#E7530,New England Biolabs)进行RNA高通量测序文库构建，由广州基迪奥生物科技有限公司使用Illumina HiSeqTM4000测序平台进行测序。对高通量测序下机的Raw reads通过fastp[13]进行质控(去除含adapter的reads、含N比例大于10%的reads和低质量的reads)，得到高质量Clean reads。利用Trinity[14]软件对高质量的Clean reads进行序列组装，然后通过N50 数值和基因的数量进行组装完整性评估。组装出来的Unigene序列通过blastx比对到非冗余蛋白库(Non-Redundant Protein Sequence Database,NR)，京都基因与基因组百科全书(Kyoto Encyclopedia of Genes and Genomes,KEGG)，同源蛋白簇(Clusters of Orthologous Groups of proteins/euKaryotic Ortholog Groups,COG/KOG)和SwissProt 4大数据库。

1.4 差异表达基因分析

采用DESeq2[15]软件对基因差异表达进行分析，基于差异分析结果，筛选错误发现率(False discovery rate,FDR)<0.05，且归一化差异倍数|log2FoldChange|>1的基因作为差异表达基因(Differentially expressed genes,DEGs)。然后通过基因本体论(Gene ontology,GO)[16]和KEGG[17]2大数据库对DEGs进行富集分析，符合FDR<0.05的GO条目和KEGG通路被定义为显著富集的条目和通路。

1.5 实时荧光定量PCR分析

为了验证RNA-seq数据，以高加索三叶草Actin(unigen 0052143)基因为内参基因，从对照组和骤然降温组中随机挑选的10个DEGs进行实时荧光定量PCR (Quantitative Real-time PCR,qRT-PCR)分析，引物设计见表1，各处理3次重复，采用相对定量法(2-ΔΔCt)[18]，计算基因的相对表达量。

表1 qRT-PCR所用引物

2 结果与分析

2.1 测序数据评价与功能注释

从转录组测序看来，各组的Raw reads均在4 292万条以上，Clean reads在4 233万条以上，过滤后总共得到67.31 G的Clean bases，GC含量在43.14%～44.64%，Q20超过 96.87%，Q30超过 91.67%(表2)。通过组装，得到59 257个Unigene，平均长度为961 bp，N50长1 678 bp(表3)。这些结果表明，测序数据的质量可用于后续分析。

表2 数据产出质量情况

表3 组装结果统计分析

由表4所示，对组装得到的59 257个Unigene进行基本功能注释，共有37 248个(62.86%)Unigene在NR，KEGG，KOG，SwissProt 4大数据库中注释成功，分别为37 168个(62.72%)，32 706个(55.19%)，17 636个(29.76%)，21 790个(36.77%)。其中，与NR数据库进行比对后(图1)，发现高加索三叶草与红三叶(Trifoliumpratense)相似度最高，同源基因相似率达(24.90%)；其次为蒺藜苜蓿(Medicagotruncatula)(相似率15.67%)和鹰嘴豆(Cicerarietinum)(相似率5.04%)。

表4 Unigene功能注释情况

图1 同源物种分布图

2.2 差异表达基因(DEGs)筛选与比较

为了鉴定高加索三叶草在不同降温模式下的DEGs，我们通过FDR<0.05且|log2FoldChange|>1的阈值进行筛选。如图2所示，与CK相比，在GC组和SC组中分别选出8 020个DEGs(4 831个上调表达，3 189个下调表达)和6 289个DEGs(3 390个上调表达，2 899个下调表达)。另外，对2处理组的DEGs进行韦恩分析(图3)，发现GC组特有5 314个DEGs，SC组特有3 583个DEGs；2处理组共有 2 706个DEGs，推测这些共有的DEGs是高加索三叶草响应低温胁迫的潜在抗寒基因。

图2 不同降温模式下DEGs的火山图

图3 不同降温模式下DEGs的韦恩图

2.3 DEGs的GO功能富集分析

为了研究高加索三叶草在不同降温模式下的DEGs涉及的生物学功能，我们对GC和SC 2处理组的DEGs分别进行了GO富集分析。由图4和表5可知，GC和SC组的DEGs均被富集到生物过程(Biological process)、分子功能(Molecular function)、和细胞组分(Cellular component)3个功能类别。生物过程中，2处理组的DGEs均集中在代谢过程(Metabolic process)、细胞过程(Cellular process)和单有机体过程(Single-organism process)；在GC组中，显著性排名前3的GO条目是刺激反应(Response to stimulus)、化学反应(Response to chemical)、内源性刺激响应(Response to endogenous stimulus)；在SC组中，排名前3的是酸性化学响应(Response to acid chemical)、外界刺激的响应(Response to external stimulus)、内源性刺激响应(Response to endogenous stimulus)。分子功能中，2处理组的DGEs均集中在催化活性(Catalytic activity)和结合(Binding)；在GC组中，排名前3的是氧化还原酶活性(Oxidoreductase activity)、调节区DNA结合(Regulatory region DNA binding)和调节区核酸结合(Regulatory region nucleic acid binding)；在SC组中，排名前3的是氧化还原酶活性(Oxidoreductase activity)、转移酶活性，转移己糖基(Transferase activity,transferring hexosyl groups)和葡糖基转移酶活性(Glucosyltransferase activity)。细胞组分中，2处理组的DGEs均集中在细胞(Cell)、细胞部分(Cell part)和细胞器(Organelle)；在GC组中，排名前3的是叶绿体(Chloroplast)、光合膜(Photosynthetic membrane)和光合系统(Photosystem)；在SC组中，排名前3的是细胞壁(Cell wall)、外部封装结构(External encapsulating structure)和细胞外周(Cell periphery)。

图4 不同降温模式下DEGs的GO富集分析

表5 不同降温模式下DEGs排名前3的GO条目

2.4 DEGs的KEGG富集分析

为了探究高加索三叶草响应不同降温模式的代谢通路情况，我们对GC和SC 2处理组的DEGs分别进行了KEGG通路富集分析。如图5所示，GC组的DEGs在光合作用(Photosynthesis)、光合作用天线蛋白(Photosynthesis-antenna proteins)和淀粉和蔗糖代谢(Starch and sucrose metabolism)等通路中显著富集；SC组的DEGs在淀粉和蔗糖代谢(Starch and sucrose metabolism)、玉米素生物合成(Zeatin biosynthesis)、倍半萜和三萜生物合成(Sesquiterpenoid and triterpenoid biosynthesis)、二萜生物合成(Diterpenoid biosynthesis)、类黄酮生物合成(Flavonoid biosynthesis)和苯丙烷生物合成(Phenylpropanoid biosynthesis)等通路中显著富集。

图5 不同降温模式下DEGs的KEGG富集分析

2.5 抗寒相关基因分析

由韦恩分析结果推测，共有的2 706个DEGs可能是不同降温模式下高加索三叶草潜在的抗寒基因，因此对其进行GO与KEGG富集分析(图6)。GO富集分析表明这些DEGs在生物过程、分子功能和细胞成分3类别中共有91条显著富集的GO条目，包括激素应答(Response to hormone)、信号转导(Signal transduction)、二糖响应(Response to disaccharide)、非生物刺激响应(Response to abiotic stimulus)和二糖代谢过程(Disaccharide metabolic process)等。KEGG富集分析表明这些DEGs仅在淀粉和蔗糖代谢通路中显著富集；此外，植物激素信号转导通路中所富集的DEGs数目最多。结合GO与KEGG富集分析结果推测淀粉和蔗糖代谢、植物激素信号转导2条通路可能是参与高加索三叶草低温胁迫响应的重要代谢途径。为进一步研究其响应低温胁迫的分子机制，本研究对这2条通路中特异上调的DEGs进行了重点分析。经过筛选，得到16个潜在的抗寒相关基因(表6)。

图6 抗寒相关基因的GO和KEGG富集分析

表6 高加索三叶草抗寒关键基因

2.6 抗寒相关转录因子分析

对由韦恩分析得到的共有DEGs进行转录因子分析，筛选到了属于29个转录因子家族的150个转录因子(Transcription factors,TFs)(图7)。其中，含转录因子数量最多的家族为AP2/ERF，其次为ZFP，MYB，WRKY，NAC等。进一步对GC和SC 2处理组中的这150个转录因子的表达量进行了统计分析，发现在GC组中有105个上调，45个下调；在SC组中有102个上调，48个下调；而在2处理组中均上调的有93个，均下调的有36个。另外，通过GO富集分析筛选出7个响应低温(Response to cold)且上调表达的转录因子，其中，AP2/ERF家族的DREB1/CBF亚家族有6个(3个DREB1A、2个DREB1C、1个DREB1F)，ZFP家族中1个(1个ZAT11)。

图7 抗寒相关转录因子分析

2.7 实时荧光定量PCR验证

利用实时荧光定量PCR技术检测随机挑选的 10 个DEGs的表达模式结果可知(图8)，这些基因的表达趋势与RNA-Seq结果一致，说明RNA-Seq数据真实可靠。

3 讨论

温度是影响所有生物的重要环境因素，在低温胁迫下可以正常生长的植物种类较少，低温胁迫严重改变了植物的生长、发育、生产力和分布[19]。高加索三叶草是一种优良豆科牧草，具有强大的抗寒能力[20]，可在接近海平面以及高海拔地区正常生长[21]。本研究通过转录组测序比较了高加索三叶草在不同降温模式下的差异表达基因数量和关键通路，还分析了潜在的抗寒基因和转录因子,这些结果将有助于了解高加索三叶草响应低温胁迫的分子调控机制。

3.1 高加索三叶草响应不同降温模式的转录差异

本研究由高加索三叶草的RNA-seq数据组装得到59 257个Unigene，其中62.86%的Unigene在4大数据库中注释成功。同时，与NR数据库比对后，发现高加索三叶草与红三叶的同源相似率最高(24.90%)。与CK相比，在GC和SC 2处理组中均筛选到大量的DEGs，而GC组的DEGs数量高于SC组，且前者特异表达的DEGs也高于后者。这一结果说明不同降温模式的低温胁迫均可诱导高加索三叶草表达大量的基因来缓解低温伤害，而缓慢降温过程可诱导表达更多的基因和特异基因来抵御低温。菊花(Dendranthemagrandiflorum)经冷驯化后再降温处理的差异表达和特异表达基因数量均高于直接降温处理的数量[22]，这与本研究结果一致。

富集分析表明高加索三叶草面临不同降温模式的低温胁迫时，在代谢通路方面会作出不同的响应。相关研究表明短时间的低温会破坏光能的吸收与转运[23]，本研究中GC组的DEGs在光合作用通路中显著富集，推测低温胁迫严重影响高加索三叶草在缓慢降温过程中的光合作用。另一方面，SC组的DEGs在倍半萜和三萜生物合成、二萜生物合成通路中显著富集。植物对短期和长期的环境温度变化均能迅速调节萜类化合物合成途径[24]，由此推测高加索三叶草能通过快速调节萜类化合物生物合成通路来应对骤然降温。此外，SC组的DEGs在类黄酮生物合成和苯丙烷生物合成通路中显著富集，这和冷冻胁迫下油菜(BrassicanapusL.)的DEGs富集结果相似[12]。苯丙烷生物合成途径在植物面临低温胁迫时发挥重要作用[25]，在该途径和类黄酮生物合成途径中，某些次生代谢物质的积累可以增强植物清除低温胁迫产生的活性氧的能力，是一种积极地自我保护机制[26]，据此推测苯丙烷、类黄酮等次生代谢产物合成途径在高加索三叶草响应骤然降温过程中起着关键作用。值得注意的是，2处理组的DEGs均在淀粉和蔗糖代谢中显著富集，结果与Wang等[22]对菊花响应冷驯化或非冷驯化后的冷冻胁迫转录组分析结果类似，可看出淀粉和蔗糖代谢在高加索三叶草响应不同降温模式的低温胁迫时均具有重要作用。

3.2 高加索三叶草关键抗寒基因

糖在植物中不仅用于能量代谢，也对植物生长发育、代谢调控及抵抗胁迫等方面具有重要调节作用[27]。植物面临低温胁迫时，糖可以提供碳源和底物来诱导其他相关的抗寒性和生理生化过程，有助于增强抗寒性，还可以避免低温引起的蛋白质凝固，进一步提高植物的抗寒性[28-29]。本研究通过对GC和SC 2处理组共有DEGs进行富集分析，结果发现这些DEGs仅在淀粉和蔗糖代谢通路中显著富集。在此通路中，筛选出一个编码β-淀粉酶(Beta-amylase,BAM)的基因BAM3，其表达量在低温胁迫后显著上调。BAM是一种重要的淀粉水解酶，分解淀粉生成麦芽糖[30]。温度胁迫下，烟草(Nicotianatabacum)过度表达PtrBAM1能提高BAM活性、促进淀粉降解、增加麦芽糖和可溶性糖的含量[31]。本研究还筛选出一个上调的编码β-呋喃果糖苷酶(Beta-fructofuranosidase,β-FFase)或称转化酶(Invertase,Inv)的基因BFRUCT1。Inv能将蔗糖水解为葡萄糖和果糖[32]，也能通过调节渗透压改变植物细胞的生理状态[33-34]。此外，蔗糖合成酶(Sucrose synthase,SuSy)能将蔗糖分解为果糖和尿苷二磷酸葡萄糖[35]；Wang等[36]研究发现低温胁迫下，SUS基因在烟草中上调表达，这与我们的研究结果一致。由此推测，高加索三叶草通过上调表达这些基因来促进可溶性糖的合成，从而提高自身抗寒性。

植物激素不仅能够调控植物正常的生长发育过程，在植物抵御低温等外界环境胁迫中也发挥着重要作用[37]。本研究中，2处理组的DEGs在植物激素信号转导通路富集数量最多，这表明激素信号通路在高加索三叶草对低温的响应中起着关键作用。这些DEGs主要富集在生长素和脱落酸2种信号途径中。其中，涉及生长素信号通路的4个生长素早期响应基因(IAA12，IAA26，GH3.1，SAUR32)的表达量均上调，在水稻(Oryzasativa)响应低温胁迫的研究中也出现类似结果[38]。据报道，拟南芥(Arabidopsis)的一个GH3的功能获得型突变体wes1-D表现明显的抗冻表型[39]。另外，本研究还发现2个上调的涉及脱落酸信号通路的SAPK2基因，均编码丝氨酸/苏氨酸蛋白激酶(Serine/threonine protein kinase,SAPK/SnRK)。SnRK2蛋白激酶家族成员的上调表达促进AREB/ABFs转录因子的磷酸化，从而促进低温胁迫信号转导和抵御环境胁迫[40-41]。可以看出，高加索三叶草可能通过上调表达这些与激素相关的基因来促进信号转导，从而抵御低温胁迫。

3.3 高加索三叶草主要抗寒转录因子

转录因子在植物低温应答过程中起着关键的调控作用[42]。目前已经鉴定了一系列参与调控植物低温胁迫应答的转录因子，包括AP2/ERF，MYB，bHLH，ZFP，NAC，WRKY等转录因子家族成员[43]。本研究中，从GC和SC 2处理组共有的DEGs中筛选到150个转录因子。在这些转录因子中，AP2/ERF家族成员最多(28个)，推测它们在高加索三叶草对低温胁迫的反应中起重要作用。AP2/ERF家族转录因子，尤其是DREB1/CBF亚家族转录因子被认为是植物冷驯化的重要调节因子，它们能大量调控冷诱导基因如COR的表达，从而提高植物抗寒性[44-45]。研究表明，拟南芥DREB1/CBF基因的表达受到低温胁迫诱导[46]。本研究筛选到与低温响应有关的6个上调表达且属于DREB1/CBF亚家族的转录因子。同时，还筛选到17个MYB家族的转录因子。MYB转录因子家族是植物最大的转录因子家族之一，在逆境胁迫的调控方面具有重要作用[47-48]。Yang等[49]报道，超量表达OsMYB2基因能显著增强水稻植株对低温胁迫的抗性。本研究也发现了1个编码MYB转录因子的基因MYB2，其表达量在2个比较组中均上调(未发表数据)。此外，其他参与植物非生物胁迫过程的常见转录因子家族，如WRKY，NAC，bHLH等，也在我们的研究中被发现，这与在其他植物中的报道一致[11-12]。由此说明这些转录因子在高加索三叶草响应低温胁迫过程中可能具有重要的调控作用。

4 结论

高加索三叶草在面临不同降温模式的低温胁迫时具有不同的响应机制，应对缓慢降温的DEGs显著富集在光合作用、光合作用天线蛋白、淀粉和蔗糖代谢等通路。应对骤然降温的DEGs显著富集在淀粉和蔗糖代谢、倍半萜和三萜生物合成、类黄酮生物合成等通路。一些涉及淀粉和蔗糖代谢与植物激素信号转导通路的基因在两种降温模式下均上调表达，可作为潜在的抗寒基因。另外，响应低温胁迫的转录因子多分布在AP2/ERF，ZFP，MYB等家族。