无瓣海桑根响应盐胁迫的转录组分析

梁锐涛　韩维栋　杨少瑕等

关键词：无瓣海桑；盐胁迫；转录组；差异表达基因

中图分类号：S728.5 文献标志码：A 文章编号：1001-1498（2023）01-0068-11

土壤盐渍化现象几乎遍布所有的气候区域以及地理环境，是影响农业和林业生产的主要问题。土壤盐渍化使植物受到渗透胁迫、离子毒害、膜透性改变及生理代谢紊乱等危害，严重制约了植物的生长，是在植物生产与栽培上一个非常严峻的问题。因此，通过对植物抗盐的分子调控机制进行探析，挖掘耐盐相关基因，并应用于新种质创制和新品种选育，对于缓解土壤盐渍化问题具有重要的现实意义。

长期生长在滨海区域潮间带的红树林，由于其特殊的、有别于陆地和淡水植物的生长环境，演化出了一套特殊的盐逆境适应机制。植物的耐盐机制涉及渗透物质积累、蛋白积累、抗氧化酶的诱导、植物激素的诱导和光合作用路径的变化5个方面。对于红树植物，不同的红树物种具体的耐盐机制有所不同。无瓣海桑（Sonneratia apetala Buch.-Ham.）是海桑科海桑属红树乔木，具有生长快、耐淹浸、耐贫瘠、耐寒强等特点，是滨海滩涂的先锋造林树种，在维持海岸带生态平衡方面起着重要作用。先前有关其耐盐特性的研究主要集中在形态、生理生态层面，如廖岩等研究了不同盐度处理3个月下无瓣海桑根、茎、叶中的蛋白质、可溶性糖、丙二醛含量、抗氧化酶活性的变化情况，并探讨了它们的相互关系，为耐盐树种的选育提供了参考。然而，目前鲜见关于无瓣海桑耐盐分子机制方面的研究报道。

转录组测序是转录组学研究的核心技术，被广泛应用于新基因的挖掘等。目前，最新的测序技术已经发展到了第三代。三代测序技术具有单分子、超长读长、无需PCR扩增等优点，已在基因组学、转录组学、表观遗传学等领域上得到了应用。在转录组学上，三代全长转录组测序（lso-Seq）的应用，可弥补二代测序技术读段较短、无法覆盖整个转录本的缺点，从而做到从全长转录本水平上研究转录组。但三代测序技术存在碱基错误率高、成本高等缺点，目前应用并不广泛。二代测序有通量高、成本低的特点，且发展至今已取得了较好的工作基础，因此，二代和三代测序结合使用成为很多研究者的策略。如梅瑜等利用基于SMRT技术原理的PacBio平台三代测序技术，对甘葛藤（Pueraria thomsoruu Benth）进行全长转录组测序及分析，从转录数据中鉴定了黄酮类生物合成基因、转录因子、R基因和SSR标记，为甘葛藤的选育和利用提供了基础。然而，对于红树植物，当前研究仅限于对白骨壤（Avicennia marina （Forsk.） Vierh.）、秋茄（Kandeliacandel （Linn.） Druce／Kandelia obovata）、角果木（Ceriops tagal （perr.）C.B.Rob）、竹节树（Carallia brachiate （Lour.） Merr.）和老鼠簕（Acanthus ilicifolius L.）等个别物种开展了转录组学相关研究。特别是，至今对无瓣海桑仅开展了以三代全长转录组数据集为基础的盐胁迫应答转录因子挖掘工作。因此，本研究利用二代与三代结合的高通量测序技术，对盐胁迫处理前后的无瓣海桑根系进行转录组测序分析，从中鉴定并初步筛选出无瓣海桑根系响应盐胁迫的候选基因，可为深入研究无瓣海桑的耐盐基因组学、分子生物学和进一步揭示红树植物的耐盐机制提供理论基础。

1材料与方法

1.1材料与处理

本实验原材料采自于广东省湛江市特呈岛红树林自然保护区（21°09′～21°10′N，110°25′～110°27′E）。2019年10月采集无瓣海桑的种子，播种于人工配制的营养土中。1a后，选取长势一致的无瓣海桑幼苗，分别设置对照组（RCK）和处理组（RT）。处理组用500 mmol·L^-1的NaCl溶液每3d进行1次盐胁迫处理，对照组使用等量清水。对照组（RCK_a、RCK_b、RCK_c）和处理组（RT_a、RT_b、RT_c）均进行了3次生物学重复。10 d后分别取对照组和处理组根部组织液氮速冻，于-80℃冰箱保存。

1.2方法

1.2.1测序文库构建及测序 提取6个样本的RNA，质检合格后构建三代全长转录组测序文库和二代RNA-seq测序文库，并进行测序及分析。二代测序文库构建方法如下：（1）用mRNA富集法或rRNA去除法对total RNA进行处理；（2）用打断buffer把获得的RNA片段化，随机的N6引物进行反转录，合成cDNA二链形成双链DNA，再将其末端补平并5'端磷酸化；（3）连接一个3'端有凸出“T”的鼓泡状的接头；（4）通过特异的引物进行PCR扩增得到PCR产物，将其热变性成单链，再用一段桥式引物将单链DNA环化得到单链环状DNA文库；（5）利用BGISEQ-500平台测序。

PacBio测序文库构建方法如下：取6个样本的RNA，等量混合，使用SMARTerTM PCRcDNA Synthesis试剂盒将混合的RNA样品反转录合成第-链cDNA，经PCR扩增合成第二链cDNA，然后对双链DNA经过2次PCR扩增后的片段进行损伤修复、末端修复，并连接SMRT适配体后构建全长转录组测序文库，利用PacBio平台进行测序。

文库构建及测序工作由华大基因公司完成。

1.2.2测序数据组装和基因定量 为保证数据质量，需对测序得到的原始数据进行过滤、参考基因比对等前期处理。本研究使用华大自主研发的过滤软件SOAPnuke对6个样本进行过滤，去除包含接头污染、未知碱基N含量大于10%、低质量的reads。使用Bowtie2将clean reads比对到三代全长转录组数据库；利用RSEM计算各个样品的基因表达水平。样本间的相关性分析通过皮尔逊相关系数（Pearson correlation coefficient）进行。

1.2.3差异基因的选择与功能注释分析 本研究基于负二项分布原理的Deseq2模型分析方法，对2组样本进行差异检测分析，并以差异倍数值（fold change， FC）≥2或≤0.5， P-value<0.01，Q-value<0.05作为筛选差异基因的标准，并将该筛选得到的差异基因用于后续的GO和KEGG Pathway富集分析，前者通过agriGO网站（http：//systemsbiology.cau.edu.cn/agriGOv2/）进行，后者借助基因数据库网站（https：//www.kegg.jp/）完成。

2结果与分析

2.1测序质量分析

通过SOAPnuke软件对原始数据进行过滤后，如表1所示，有效读序百分率均在98%以上，Q20（碱基量≥20%）均大于96%且Q30（碱基量≥30%）均>91%，比对率均在87%以上。6个样品转录本reads覆盖度的峰值在90～100；空白对照RCK组和盐处理RT组的三次生物学重复间的相关系数极高（>0.98，图1）。以上数据均说明6个样品的测序质量较高，测序数据可靠，可用于后续的差异表达分析。

2.2样本差异基因表达情况

无瓣海桑根系在500 mmol·L^-1的NaCl胁迫10 d后，共检测到20 289个基因表达量发生了改变（图2）。2组样本间共有14 401个基因差异表达，其中，7 153个差异基因上调表达，lOg2FC最大值为11.32（isoform_13 051）；7 248个差异基因下调表达，lOg2FC最低值为-12.69（isoform_11 580）。结果表明，下调表达的差异基因数量多于上调的数量。

2.3差异基因的GO分类和显著性富集分析

将差异基因比对到Gene Ontology （GO）数据库，共有11068个差异基因注释到47条GO条目，注释率为76.9%，涉及3大主要功能分类：生物过程、细胞组分和分子功能，占比分别为25.02%，43.76%，31.21%，由此可知，大多数差异基因的功能与细胞组分具有较大的相关性。根据GO功能注释分类结果（图3），这些差异基因在参与生物过程中，主要集中在细胞过程和代谢过程；在细胞组分中，主要集中在细胞、膜、膜部分；在分子功能中，主要集中在连接、催化活性。表明在高盐胁迫下，上述的亚类中涉及到的基因可能在无瓣海桑根系中响应强烈。

2.4差异基因的KEGG富集分析

在KEGG分析结果中，共有6189个差异基因富集到134条通路，以Q-value<0.05作为显著富集标准，共有14条显著富集的通路（图4），其中，富集差异基因数量最多的是植物激素信号转导（ko04075），涉及512个基因，最显著富集的植物一病原体相互作用（ko04626）涉及367个基因。富集分析结果表明，无瓣海桑根系基因主要通过14条显著通路参与盐胁迫的响应，并且其中涉及的基因可能对无瓣海桑根系耐盐性的调控发挥作用，关注这些显著通路，有助于对无瓣海桑根系耐盐候选基因的筛选。

2.5响应盐胁迫的差异表达基因功能分析

通过进一步对差异基因进行功能注释分析，本研究共筛选出89个无瓣海桑根系抗盐候选功能基因，涉及12条KEGG通路（表2，图5）。

2.5.1抗氧化基因的差异表达 正常情况下，植物体内存在着活性氧（reactive oxygen species，ROS）产生和清除的动态平衡，但外界的不良环境会打破这个平衡，导致ROS的积累，产生毒害。植物体内存在2种ROS清除机制——酶促清除和非酶促清除，前者主要涉及超氧化物歧化酶（SOD）、过氧化氢酶（CAT）等，后者主要涉及抗坏血酸、类黄酮物质等。本研究发现6类共24个活性氧清除基因（图6），其中，上调表达基因数量最多的是过氧化物酶基因，总数为7，表达量差异变化最大的为铜锌家族超氧化物歧化酶基因，差异倍数达到6.46。这些上调表达的基因有利于无瓣海桑在盐胁迫前期，合成多种抗氧化酶类物质，提高体内抗氧化和ROS清除的能力，增强耐受性。

2.5.2渗透调节物质基因的差异表达 植物在对盐胁迫的适应过程中，细胞会积累一定量的可溶性有机物质，作为渗透调节物质参与渗透调节，以适应外界的低水势。在本研究中，发现10种共22个编码生成有机渗透调节物质的基因（图6），涉及水苏糖、棉子糖、蔗糖、海藻糖、脯氨酸、甜菜碱、热激蛋白等渗透调节物质。在这些基因中，部分基因存在上调表达，表明在盐胁迫下，无瓣海桑可通过合成一些渗透调节物质，以调节盐渍环境下植物细胞内的渗透势。

2.5.3植物激素基因的差异表达 植物体内的多种激素在其生长发育、代谢等生理过程中发挥着重要的调控作用。本研究发现7种共19个涉及2种植物激素（IAA与ABA）的基因（图6），其中，IAA途径相关的有生长素应答因子（auxin response factor，ARF）、生长素应答GH3家族基因（auxin responsive GH3 gene family， GH3）、SAUR家族蛋白（SAUR family protein，SAUR）和生长素响应蛋白（auxin-responsive protein IAA，IAA）共11个基因；涉及ABA途径的有脱落酸应答元件结合因子（ABA responsive element bindingfactor，ABF）、脱落酸受体PYR/PYL家族（abscisic acid receptor PYR/PYL family，PYL）和蛋白磷酸酶2C（protein phosphatase 2C，PP2C）共8个基因。IAA相关基因除GH3基因均上调表达外，其他同时存在上下调表达，而ABA相关基因大部分下调表达。上述结果表明，植物激素能够在无瓣海桑的信号转导途径中发挥作用，不同激素共同调控盐胁迫响应过程。

2.5.4植物蛋白激酶基因的差异表达 蛋白激酶可调控植物的生长，在逆境中起到调控感受胁迫信号、启动各种非生物逆境响应的作用。本研究发现6种共10个蛋白激酶相关基因（图6），其中，2种促分裂原活化蛋白激酶基因（mitogen-activated protein kinases，MAPK）MPK4、MPK6和2种丝氨酸，苏氨酸蛋白激酶基因PBS1、CTR1均上调表达，1种促分裂原活化蛋白激酶基因（mitogen activated protein kinase，MEKK）MEKK1和1种钙依赖性蛋白激酶基因（calcium-dependent protein kinase， CDPK）中均有2个上调表达，1个下调表达。可见这些蛋白激酶在无瓣海桑根系响应高盐胁迫时，发挥着尤为重要的调控作用。

2.6响应盐胁迫的转录因子分析

转录因子是一类功能蛋白质，具有调控基因表达的作用。当植物受到逆境胁迫时，转录因子通过与相应的顺式作用元件相互作用，调控并减轻逆境胁迫带来的伤害，以适应各种逆境。本研究发现8种共14个差异表达的转录因子（图6），涉及5个家族，分别为WRKY家族、MYB家族、GRAS家族、bHLH家族和EREBP家族，其中，WRKY家族发现的成员数量最多，暗示WRKY家族在无瓣海桑响应盐胁迫中起着重要的调控作用。

3讨论

植物对盐胁迫的响应是一个多基因参与的复杂过程，从微观上的基因响应到宏观上生理的变化，经历了一系列的信号传递与调节过程，最终适应盐逆境。通过转录组测序技术对植物盐胁迫下的转录调控机制进行研究分析，能够快速准确地初步筛选出其中的耐盐相关基因，为进一步利用分子生物学手段深入探究其耐盐机制奠定数据基础。然而，对于红树植物，当前研究对无瓣海桑仅限于利用转录组学技术开展了盐胁迫相关转录因子的挖掘，而对转录因子以外的其他可能参与调控盐逆境适应过程的基因并未关注。本研究联合二代和三代转录组测序技术，获得了无瓣海桑对照组及500mmol·L^-1盐胁迫处理10d后的根系转录组数据，在此基础上进行生物信息学分析，最终筛选出89个无瓣海桑根抗盐候选基因，其中一些基因对盐胁迫的响应起着十分重要的作用。

在植物的ROS酶促清除机制中，SOD、过氧化物酶（POD）、CAT是常见的酶。研究表明，NaCl胁迫后耐盐性强的A34株系葡萄根系中的SOD、CAT、POD活性表现出快速升高并在较高水平后小幅度降低的趋势，这有利于维持活性氧代谢的平衡，缓解细胞膜脂过氧化。鲁琳等通过高通量测序技术对盐胁迫下花烟草（Nicotiana alata Link et Otto）的活性氧清除基因进行挖掘，发现编码SOD、POD、CAT的基因均发生了显著的上调表达，同时存在下调表达的基因。本研究也发现了类似的结果，在SOD中，SOD1 2个上调，1个下调，SOD2 1个上调表达；CAT5个基因全部上调表达；在POD中，存在较多下调表达的基因。上述结果推测可能是由于植物体内的抗氧化机制存在着正负2种调控方式，且不同抗氧化酶在不同植物中的调控方式可能有所不同；对于无瓣海桑来说，CAT可能是以正向调控为主导，而POD可能主要是以负向调控为主。

在逆境下，植物可通过调节海藻糖的含量，稳定蛋白质的结构与功能，提高在逆境中的生存能力。李辉等克隆了海藻糖合成关键酶基因海藻糖-β-磷酸磷酸酶基因（trehalose-β-phosphate phosphatase，TPP），并研究其在非生物胁迫下的作用，结果表明该基因表达显著上调。研究表明，高盐胁迫主要促进了龙须菜（Asparagus schoberioides Kunth）海藻糖-6-磷酸合成酶/磷酸酶1基因（trehalose 6-phosphatesynthase/phosphatase 1，TPS1）、海藻糖-6-磷酸合成酶，磷酸酶2基因（trehalose 6-phosphatesynthase/phosphatase 2，TPS2）和海藻糖-6-磷酸合成酶/磷酸酶4基因（trehalose 6-phosphatesynthase/phosphatase 4，TPS4）的表达，而渗透胁迫抑制了TPS1、TPS2和海藻糖-6-磷酸合成酶，磷酸酶3基因（trehalose 6-phosphate synthase/phosphatase 3，TPS3）的表达。本研究发现2种共7个合成海藻糖的关键基因：海藻糖-6-磷酸磷酸酶基因（otsB）和海藻糖-6-磷酸合成酶基因（TPS）。在这些基因中，3个表现出下调表达，可能是由于盐胁迫产生的渗透胁迫抑制了这些基因的表达。

脯氨酸是植物中常见的有机渗透调节物质，而脯氨酸的积累是植物应对盐胁迫的普遍反应。在大多数植物物种中，δ-1-吡咯啉-5-羧酸合成酶（P5CS）催化谷氨酸生物合成脯氨酸的第一步。在Funck等的研究中，发现在拟南芥（Arabidopsis tha/iana （L.） Heynh）NaCl处理的p5cs2突变体中，几乎没有Na+的积累，植株既没有褪绿，也没有光合作用减弱，可见P5CS2介导的脯氨酸合成在调节叶片Na+积累和抗盐胁迫方面具有重要调控功能。本研究中发现2个P5CS基因，1个上调，1个下调，表明P5CS基因在渗透调节中发挥了作用。

脱落酸在盐、干旱、低温等非生物胁迫下会在植物体内被诱导产生，同时启动相关的信号途径，以响应非生物胁迫。ABA信号途径由ABA受体、PP2C、SnRK2和ABF转录因子四部分组成。研究表明，盐处理抑制了拟南芥中PP2C31基因的表达，且在不同时期盐处理下PP2C31基因出现了不同程度的下调表达，说明PP2C31蛋白在拟南芥响应盐胁迫中起负调控作用。本研究发现4个PP2C基因均下调表达，推测PP2C蛋白负调控无瓣海桑的盐胁迫响应。

丝裂原活化蛋白激酶（ MAPK）途径在盐、高温、病原菌等逆境的应答中起着重要作用，它通过MAPKKK→MAPKK→MAPK逐级传递放大胁迫信号，由下游的MAPK激活调控相应基因的表达，引起植物对逆境的应答。在野生大麦（Hordeum spontaneumL.）中，MEKK1-MKK2-MPK4/6联级系统参与了盐胁迫的响应。拟南芥的MAPK联级中多个MAPK基因在多种非生物胁迫的响应中起着重要作用。汪芳珍等通过对沙芥根进行6h50 mmol·L^-1的盐处理实验，结果表明，共有11个MAPK/MAPKK/MAPKKK蛋白激酶基因的表达显著上调。本研究中发现MEKK1、MPK4、MPK6基因大都上调表达，推测这些基因在无瓣海桑的MAPK联级途径响应盐胁迫中起着重要的作用，但无瓣海桑具体响应盐胁迫的MAPK联级系统中的基因及机制，还有待进一步研究。

植物体转录因子在盐胁迫调控下游机制中发挥着重要作用。参与植物盐胁迫响应的转录因子家族有WRKY、MYB、bHLH、NAC等。在盐胁迫的响应中，WRKY家族转录因子存在着正负2种调控机制。在拟南芥中，RtWRKY23可提高自身过氧化物酶水平，以增强盐胁迫耐受性。转基因烟草中过表达的CaWRKY27抑制了抗氧化酶合成途径相关基因的表达，在响应盐胁迫中起负调控作用。本研究发现WRKY22、WRK Y25均呈上调表达，说明这些转录因子可能正向调控无瓣海桑的耐盐性。此外，WRKY2、WRKY33呈现上下调表达，这表明了WRKY家族转录因子参与无瓣海桑耐盐机制方式的复杂性。

4结论

本研究通过二代RNA-seq和三代全长转录组测序技术的联合使用对无瓣海桑响应盐逆境的转录调控机制进行了解析，获得了盐逆境下无瓣海桑根系的差异表达基因及功能注释信息，并鉴定出一批与活性氧清除、渗透调节、植物激素、蛋白激酶、转录因子等相关的潜在耐盐候选基因，研究结果对揭示红树植物的抗盐机制、耐盐型红树新型品种选育及沿海生态环境保护具有重要意义。