基于RNA-Seq挖掘桃低温应答相关基因的研究

李小兰 郝兰兰 张 帆 王 鸿,*

(1 甘肃省农业科学院林果花卉研究所，甘肃 兰州 730070；2 甘肃农业大学园艺学院，甘肃 兰州 730070)

桃(PrunuspersicaL.)原产于我国[1]，属蔷薇科(Rosaceae)李属(Prunus)桃亚属(Amygdalus)，是我国的主栽果树之一[2]。温度是影响桃树生长发育以及地理分布的一个重要环境因子。我国北方地区冬季气候严寒,如果发生冻害，会给果树生产造成严重的损失。此外，“倒春寒”的发生，也会使花芽受到损伤，严重影响果树产量。近年来，转录组测序RNA-Seq为植物响应低温胁迫的研究提供了一个新的视角，可以高效、快捷地反映生物在某一特定坏境和时间下基因的表达情况[3]，该技术已成为解析植物应答低温胁迫的分子机制以及挖掘功能基因的重要手段[4-5]。

植物不同的器官组织在低温诱导下表现出不同的表型和抗寒反应[6]。张丽[7]通过对低温储藏过程中两种桃果实进行转录组分析，发现ID为Prupe.1G220700的基因可能在桃的低温储藏过程中发挥重要作用。前人对已经开花的桃树枝条进行不同时间的低温处理，随后分离雌蕊进行高通量测序，发现与碳水化合物水解和低温响应相关的基因表达量显著上调[8]。Niu等[9]和Renaut等[10]分别对一年生枝条和树干皮层组织进行不同程度的低温诱导，并对获得的差异基因进行生物信息学分析，发现与碳水化合物代谢途径以及信号转导途径相关的基因在桃树抗寒过程中发挥重要作用。前人已经从不同层面揭示桃树的花[8]、果实[11-12]、枝条[9]和树干皮层组织[10]在低温诱导下的基因表达水平变化，但针对低温胁迫下桃叶片应激反应的研究较少。叶片作为重要的营养器官，在植物生长发育过程中发挥重要作用。因此，研究低温胁迫下桃叶片的抗寒反应，明确响应低温胁迫的分子调控机制，可为抗寒桃品种的培育提供理论基础。

甘肃河西走廊地区平均最低气温约为-27℃[13]，引进的桃品种无法适应这种环境。然而，属于李光桃抗寒品种群的丁家坝油桃品种能很好地适应这种环境，进化出优良的耐寒特性。本研究选取低温耐受型丁家坝[14]和敏感型加纳岩[14]桃叶片作为试验材料，比较其在低温胁迫下的生理生化水平变化，并利用生物信息学方法对丁家坝和加纳岩桃叶片在正常和低温条件下的转录组进行比较分析，鉴定其响应低温的相关基因及代谢途径，并对相关代谢途径进行深入研究，以期阐明桃树叶片抗寒的分子机制，为研究桃抗寒的分子机制提供新线索，也为桃抗寒育种提供新的候选基因。

1 材料与方法

1.1 试验材料

试验材料取自甘肃省农业科学院桃园中长势优良的冷耐受型丁家坝(D)和冷敏感型加纳岩(K)七年生桃树叶片。

1.2 试验处理

2021年7月下旬，于甘肃省农业科学院桃园各采集12枝粗细相近、长短相同、朝北直立生长的丁家坝和加纳岩一年生带叶枝条。将每个品种的枝条随机分为4组，用保鲜袋密封，放在4℃培养箱中进行低温处理。在处理0(D0/K0)、2(D2/K2)、6(D6/K6)和12 h(D12/K12)时，每个品种随机选取一组枝条采集叶片，每组3个枝条分别作为3个重复。选取从上往下数第6到第9片叶作为试验材料，立即用液氮冷冻(相对电导率用鲜叶测定)，在-80℃下保存备用。

1.3 叶片生理生化指标测定

采用磺基水杨酸—酸性茚三酮法测定游离脯氨酸(Proline, Pro)含量[15]；采用蒽酮比色法测定可溶性糖含量[16]；考马斯亮蓝G-250染色法测定可溶性蛋白含量[17]；采用硫代巴比妥酸(thiobarbituric acid, TBA)反应测定丙二醛(malondiadehyde, MDA)含量[18]，并使用Yun等[19]描述的方法计算电解质渗透指数(electrolyte leakage index, ELI)，使用POD试剂盒(北京索莱宝科技有限公司)测定过氧化物酶含量(peroxidase, POD)。数据采用Excel 2010和SPSS 20.0进行整理分析。

1.4 RNA提取、文库构建及测序

用植物总RNA提取试剂盒(中国天根生物技术有限公司，北京)提取总RNA。转录组测序由南京派森诺基因科技有限公司完成。采用琼脂糖凝胶电泳法检测RNA的完整性及是否存在DNA污染。使用纳米分光光度计(IMPLEN，德国)检测RNA纯度。建库中使用的建库试剂盒为NEBNext® UltraTM RNA Library Prep Kit(Illumina，美国)。库检合格后，对不同文库进行Illumina HiSeq X Ten PE150高通量测序。

1.5 数据处理及拼装

对获得的原始读数(raw reads)进行过滤，去除带接头(adapter)的reads、含N的reads以及低质量的reads。同时，对干净读数(clean reads)进行Q20和Q30含量计算。使用HISAT2 v2.1.0[20]将配对末端clean reads与参考基因组比对。使用featureCounts(1.5.0-p3)计算映射到每个基因的读数，然后根据基因的长度计算每个基因的FPKM(每千个碱基转录每百万映射读取的fragments，Fragments Per Kilobase of exon model per Million mapped fragments)。

1.6 差异基因富集分析

采用DESeq[21]对基因表达进行差异分析，|log2FoldChange|>1、P<0.05作为差异基因筛选条件。使用topGO[22]进行GO富集分析，找出差异基因显著富集的GO term。使用clusterProfiler[23]软件进行京都基因与基因组百科全书(Kyoto Encyclopedia of Genes and Genomes，KEGG)通路富集分析，探究P<0.05的显著富集通路。

1.7 实时荧光定量PCR(quantitative real-time PCR，qRT-PCR)表达量验证

选取15个差异基因通过qRT-PCR对测序结果进行验证，使用与测序同批次的RNA样本，按照Prime Script TM RT reagent Kit反转录试剂盒(TaKaRa, 大连)说明书合成cDNA。使用Primer 3.0在线设计引物，由上海生物工程股份有限公司合成，Actin作为内参基因(表1)。使用Light Cycler®96 Real-Time PCR System PCR仪(Roche，瑞士)进行扩增，扩增体系为20 μL：上下游引物各1 μL，cDNA 2 μL，ddH2O 6 μL，TB Green-II 10 μL。反应程序：95℃变性5 min；95℃变性10 s，60℃退火10 s，72℃延伸10 s，40个循环。本试验进行3次生物学重复，目的基因表达量的统计方法采用相对定量法(2-△△CT)。

表1 本研究所需的引物信息

2 结果与分析

2.1 桃叶片在冷处理下的生理生化变化

冷处理下，POD可以阻止植物细胞氧化损伤，减轻胁迫危害。由图1-A可知，D和K桃叶片在冷处理12 h的POD活性与对照相比分别上升了58.94%和33.66%。

MDA含量和ELI与植物细胞膜受损伤程度及植物抗逆性强弱相关。D和K叶片的MDA含量(图1-B)和ELI(图1-C)均随着冷处理时间的延长呈上升趋势，且D的增长幅度小于K，表明D的细胞膜不易被破坏，抗寒性强。

注：不同小写字母表示同一品种在不同时间低温处理下差异显著(P<0.05)；*和**分别代表不同品种在相同时间低温处理下差异显著(P<0.05)和极显著(P<0.01)。

脯氨酸(Pro)、可溶性糖和可溶性蛋白含量的积累可以增强植物对逆境胁迫的耐受性，是植物在逆境下采取的自我保护措施。总体来看，D和K的脯氨酸、可溶性糖以及可溶性蛋白含量在冷处理12 h时与对照相比显著上升，且D的增长幅度大于K。

2.2 转录组测序数据统计与质量评估

分别对D和K两种桃叶片经4℃处理0、2、6和12 h并进行取样，每个样本进行3次生物学重复，共建立24个cDNA文库(表2)。经过高通量测序后，样本产生的原始序列(raw reads)数为4 300万～5 700万，过滤后的序列(clean reads)为4 000万～5 300万。比对到参考基因组(GCF_000 346 465.2)的序列总数占clean reads的87.59%～97.00%。对每个样本的clean reads进行统计后发现Q20平均达到97.61%，Q30值平均达到 93.27%。以上结果说明测序所得的cDNA文库质量较高，可以进一步进行生物信息学分析。

表2 测序数据的统计与质量评估

2.3 低温胁迫下差异表达基因的分析

为了筛选与低温响应相关的基因，分别对两种桃叶片不同时间段低温处理产生的差异基因进行比较分析(图2)。D在低温胁迫2、6和12 h时与对照相比分别有659、1 632和2 458个基因上调，534、1 064和1 951个基因下调(图2-A)，3个处理时间段共同表达的DEGs有569个(图2-B)。与对照相比，K在低温胁迫2、6和12 h时分别有407、1 742和2 016个基因上调，506、946和1 815个基因下调(图2-A)，共同表达的差异基因有505个(图2-B)。

图2 两种桃叶片在不同时间低温处理下的差异基因柱形图(A)和韦恩图(B)

2.4 差异基因GO富集分析

在D对比组中，569个共同表达的差异基因富集在1 097个GO项(基因功能国际标准分类体系，Gene Ontology)，主要包括分子功能(288个)、生物过程(714个)和细胞组分(95个)。根据错误发现率(false discovery rate，FDR)值选择富集最显著的20个GO项作图(图3-A)。显著富集中的前5个GO项为：属于分子功能的转录调节活性(GO:0140110)、DNA结合转录因子活性(GO:0003700)和分子功能(GO:0003674)以及属于细胞组分的膜的固有成分(GO:0031224)和属于生物过程的生物学过程(GO:0008150)。

注: 1: 转录调节活性; 2: 转录，DNA模板; 3: 序列特异性DNA结合; 4: 转录调控，DNA模板; 5: RNA代谢过程的调控; 6: RNA生物合成过程的调控; 7: 初级代谢过程的调节; 8: 含碱基的化合物代谢过程的调控; 9: 核酸模板转录的调控; 10: 分子功能; 11: 膜; 12: 膜的固有成分; 13: 膜的组成部分; 14: 杂环化合物结合; 15: DNA结合转录因子活性; 16: DNA结合; 17: 细胞组成; 18: 细胞解剖实体; 19: 生物过程; 20: 结合; 21: 转移己糖基的转移酶活性; 22: 东莨菪碱β-葡萄糖苷酶活性; 23: 细胞过程的调节; 24: 细胞生物合 成过程的调节; 25: 生物合成过程的调控; 26: 葡萄糖苷酶活性; 27: β-葡萄糖苷酶活性。

在K对比组中，505个共同表达的差异基因富集到1 021个GO项，其中分子功能(268个)、细胞组分(79个)以及生物过程(674个)。挑选FDR值最小即富集最显著的20个GO项进行展示(图3-B)。前5个显著富集的GO项为：MF的DNA结合转录因子活性(GO:0003700)、转录调节活性(GO:0140110)、序列特异性DNA结合蛋白(GO:0043565)和β-葡萄糖苷酶活性(GO:0008422)以及BP的生物学过程(GO:0008150)。

2.5 差异基因KEGG富集分析

为了进一步分析DEGs的生物学功能，采用KEGG数据库对DEGs进行通路富集分析。在D对比组中，569个差异基因涉及57个KEGG途径，包括细胞过程、遗传信息处理、环境信息处理、新陈代谢、有机体系统等五个方面，大部分途径属于新陈代谢类别。显著富集的代谢通路有植物与病原菌互作(plant-pathogen interaction)、MAPK信号通路-植物(MAPK signaling pathway-plant)和倍半萜和三萜生物合成(Sesquiterpenoid and triterpenoid biosynthesis)(图4-A)，分别涉及11、8和3个相关基因。

K中的505个DEGs被分配到47个KEGG途径中，包括细胞过程、遗传信息处理、环境信息处理、新陈代谢、有机体系统等五个方面。最显著富集的途径是氰基氨基酸代谢(cyanoamino acid metabolism)和内质网蛋白加工(protein processing in endoplasmic reticulum)(图4-B)。以上结果表明，低温胁迫对两种抗寒性不同的基因型的影响不同。比较D和K富集到的KEGG通路，重点研究D中显著富集的代谢通路，如植物与病原菌互作(plant-pathogen interaction)和MAPK信号通路-植物(MAPK signaling pathway-plant)。

注: 1: 胞吞作用; 2: MAPK信号通路-植物; 3: 植物激素信号转导; 4: 磷脂酰肌醇信号系统; 5: 内质网中的蛋白质加工; 6: 倍半萜和三萜生物合成; 7: 萜类化合物生物合成支柱; 8: 甜菜色素生物合成; 9: α-亚麻酸代谢; 10: 油菜素类固醇生物合成; 11: 酮体的合成与降解; 12: 二萜生物合成; 13: 类胡萝卜素生物合成; 14: 苯丙氨酸、酪氨酸和色氨酸生物合成; 15: 卟啉与叶绿素代谢; 16: 核黄素代谢; 17: 半乳糖代谢; 18: 硫胺素代谢; 19: 亚油酸; 20: 昼夜节律-植物; 21: 光合作用-天线蛋白; 22: 生物素代谢; 23: 类黄酮生物合成; 24: 丁酸代谢; 25: 叶酸生物合成; 26: 莨菪碱、哌啶和吡啶生物碱的生物合成; 27: 淀粉和蔗糖代谢; 28: 角质、木栓质和蜡的生物合成; 29: 植物与病原菌互作; 30: 苯丙烷的生物合成; 31: 蛋白质输出; 32: 氰基氨基酸代谢; 33: 芥子油苷生物合成; 34: N-聚糖生物合成; 35: 精氨酸和脯氨酸代谢; 36: 赖氨酸生物合成; 37: 不饱和脂肪酸生物合成; 38: 烟酸盐和烟酰胺代谢; 39: 色氨酸代谢; 40: 苯丙烷代 谢; 41: 脂肪酸生物合成; 42: 脂肪酸降解。

2.6 丁家坝抗寒相关基因

对丁家坝显著富集的代谢通路进行分析表明，在植物与病原菌互作途径，由真菌几丁质、Ca2+、细菌鞭毛等不同PAMPs激发的反应中，CaM、CDPK、CML和WRKY等相关基因转录本显著上调，促进ROS和一氧化氮(NO)的产生，参与丁家坝抗寒反应。此外，MAPK信号通路-植物代谢通路中FLS2受体识别来自细菌鞭毛蛋白的flg22肽，激发植物的MAPK级联反应，引起细胞内转录及代谢水平变化，进一步响应低温胁迫。对植物与病原菌互作和MAPK信号通路-植物途径涉及的基因进行热图分析(图5)，结果表明CPK26、CALM7、CML45、CML31、CML24、CML27、CML38、WRKY22、WRKY24和WRKY33等相关基因在丁家坝抗寒反应中发挥重要作用。

图5 丁家坝抗寒相关基因热图

2.7 qRT-PCR验证

为验证转录组测序数据的可靠性，从转录组鉴定到的差异基因中选取了15个可能参与低温响应的基因进行实时荧光定量分析。如图6所示，qRT-PCR分析的基因表达模式与RNA-Seq方法分析的基因表达模式相似，证实了RNA-Seq数据的可靠性。

注：qRT-PCR：直方图；FPKM：折线图。

3 讨论

低温是植物生长发育过程中常见的环境胁迫因子，植物通过改变其形态和生理生化特性来适应低温胁迫[24]。前人研究表明，植物对低温胁迫的响应与许多基因的表达有关[25]。转录组测序技术的快速发展为解析植物响应低温胁迫的分子机制和挖掘抗寒关键基因提供了新的途径[26-27]。本研究对冷处理下的耐冷型丁家坝和敏感型加纳岩桃叶片的生理生化水平变化进行研究，并利用转录组技术对其冷处理下的转录组进行对比和分析，以鉴定出与丁家坝抗寒性相关的调控机制及关键基因。结果表明，丁家坝和加纳岩桃叶片的电解质外渗率和MDA含量整体呈上升趋势，在冷处理6 h时出现激增。与冷处理2 h相比，6 h丁家坝和加纳岩的ELI分别增加了45%和53%，MDA含量分别增加21%和24%。这与Megha等[28]的研究中关于橄榄型油菜ELI随着冷处理时间的延长逐渐增加的结果相同；加纳岩在冷处理各时期的MDA含量均高于丁家坝，这与低温胁迫下甜瓜不同品种间的耐低温性与MDA含量成显著负相关的研究结果一致[29]。MDA含量的积累是ROS过度积累导致的，抗氧化物酶在植物活性氧解毒过程中发挥重要作用[30]。本研究中，冷处理12时，MDA增长率在与6 h相比有所下降，推测可能与POD活性在冷处理12 h时与6 h相比显著上升有关。前人研究表明，植物体内可溶性糖、可溶性蛋白和脯氨酸等重要渗透调节物质的含量与植物自身的抗寒性成显著正相关[31-32]。本研究中，丁家坝和加纳岩的可溶性蛋白、可溶性糖以及脯氨酸的含量随着冷处理时间的延长逐渐升高，且丁家坝的增长幅度高于加纳岩，与前人研究结果一致。

转录组分析结果表明，丁家坝在冷处理下的差异基因显著富集在植物与病原菌互作、MAPK信号通路-植物和倍半萜和三萜生物合成等代谢通路。信号转导途径是冷感机制和遗传反应之间的联系[33]。前人已经研究并阐明了可能与植物冷响应相关的信号通路[34]。Niu等[9]对丁家坝休眠枝在不同程度低温胁迫下的转录组研究结果表明，差异基因最显著富集的途径是植物与病原菌互作，并认为Ca2+信号转导途径在丁家坝抗寒过程中发挥重要作用。本研究中，由真菌几丁质、Ca2+、细菌鞭毛等不同的PAMPs激发的反应中，CDPK和CaM/CML的转录本显著上调。CDPK基因上调可促进ROS的积累，是植物早期防卫反应的标志之一[35]；CaM/CML的上调加速NO的产生，促发过敏性坏死反应和抗病反应[36]。此外，FLS2受体识别来自细菌鞭毛蛋白的flg22肽，激发丁家坝的早期防卫反应后进一步激发MAPK级联反应，将细胞表面的信号放大并传递到细胞内，引起细胞内的转录及代谢水平变化，从而响应低温胁迫。这与田玉珍等[37]对天汪一号苹果花芽响应不同时间低温胁迫的研究结果相同，该研究发现差异基因最显著富集的途径为植物与病原菌互作，并表明Ca2+在诱导CaM/CML上调表达参与低温响应的同时也激发MAPK级联通路中的WRKY22、WRKY33和WRKY25等基因参与低温响应。本研究中，MAPK级联反应的功能酶基因MPK3/6以及WRKY22、WRKY24和WRKY33等相关基因在冷处理下的丁家坝叶片中具有较高表达，参与丁家坝对低温的响应过程。WRKY22、WRKY24和WRKY33分别在菜心和笃斯越橘的研究中被证实响应低温胁迫[38-39]。类似的低温响应机制在Niu等[9]研究休眠枝响应不同时间低温胁迫的试验中也有印证。

4 结论

本研究对低温耐受型丁家坝和敏感型加纳岩在冷处理下的生理生化水平变化和转录组进行比较分析，发现植物与病原菌互作和MAPK信号通路以及CDPK、CaM、CML、MPK3/6、WRKY22、WRKY24和WRKY33等相关基因在丁家坝抗寒过程中发挥重要作用。研究结果有助于了解丁家坝响应冷胁迫的分子机理，为桃耐冷品种的培育提供理论依据。