盐胁迫下比拉底白刺差异表达基因的转录组分析

田 林，尹丹丹，成铁龙，夏新莉，尹伟伦

(1. 北京林业大学生物科学与技术学院，北京 100083；2. 中国林业科学研究院林木遗传育种国家重点实验室，北京 100091；3. 南京林业大学生物与环境学院，江苏 南京 210037)

土壤高盐是限制植物生长发育、地理分布的主要环境胁迫因子之一[1]。世界上有超过6%的土地被盐化，超过20%的灌溉土地遭受盐害[2-4]。盐生植物作为可以在盐化土壤中正常生存的一类植物，被广泛地用于盐碱地的生态改良和植被保护。因此，研究盐生植物的耐盐机制、发掘和选育耐盐植物对充分利用盐化土地、改善生态环境具有重要的生态和经济意义。

蒺藜科（Zygophyllaceae）白刺属（Nitraria）植物广泛生长在亚洲、非洲和欧洲盐碱地，因其较强的耐盐性，被应用于盐碱地的植被改善和生态环境改良。有关白刺耐盐的特性和生理生化机制研究已有一定的进展，研究报道白刺属可以耐受高达500 mmol·L-1的 NaCl胁迫条件[5-6]，NaCl胁迫可以显著增加 Na＋浓度，降低 K+浓度[7]。有研究发现，在盐胁迫下，唐古特白刺（N. tangutorum）可能通过将Na+积聚在液泡内进行隔离，并对其它矿质营养离子进行选择性吸收，进而控制Na+/K+比值，以提高白刺耐盐性[8]。唐古特白刺幼苗经高盐处理7 d时，叶绿素含量随盐胁迫浓度升高而呈现先升后降的变化，当NaCl浓度达 300 mmol·L-1时，叶绿素含量最高[9]；盐胁迫也提高了游离脯氨酸、可溶性蛋白和可溶性糖等渗透调节物质的含量[10-12]。低浓度NaCl处理可以提高叶片中超氧化物岐化酶（SOD）、过氧化物酶（POD）和过氧化氢酶（CAT）等抗氧化酶活性，进而提高植物对盐分胁迫的适应性[13]。

关于白刺耐盐的分子机制研究可以归结为两个方面：一方面，在500 mmol·L-1NaCl 高盐胁迫时，唐古特白刺幼苗体内涉及到的光合过程、氧化还原、逆境、防御、能量代谢、糖代谢、氨基酸代谢、信号转导、蛋白质合成、蛋白质折叠和装配、转录、膜运输、激素合成等途径的蛋白质表达发生了显著改变[4]；另一方面，200 mmol·L-1NaCl 中等盐胁迫诱导了比拉底白刺 (N. billardieri)幼苗叶片内多条代谢途径的蛋白质表达水平的改变，其中涉及到的有碳水化合物代谢、光合过程、氨基酸代谢、蛋白质降解、蛋白质折叠和组装、蛋白质合成、氧化还原、次级代谢、信号转导、细胞周期、细胞骨架、细胞防御、能量代谢、膜与运输、转录调控等[14]。上述2个白刺研究均是通过蛋白质组学的方法，不同强度盐胁迫导致了一些相同类别蛋白质的表达及变化趋势，同时也呈现出一些细微的差别。

上述研究工作为了解白刺耐盐的生理生化与分子机制奠定了一定的基础，而基因作为植物应答环境胁迫的早期响应分子，对调控植株应对胁迫生存条件并适应胁迫环境起到关键作用。为了寻找白刺应答盐胁迫的基因并解析其分子机制，本研究以比拉底白刺幼苗为材料，运用转录组测序技术挖掘200 mmol·L-1NaCl 胁迫下比拉底白刺幼苗叶片内的差异表达基因，并运用生物信息学方法分析差异表达基因的生物学功能及其参与的代谢途径。

1 材料与方法

1.1 实验材料

比拉底白刺幼苗的培养参照Tian等[14]方法，并稍作调整。将经温水浸种24 h后的比拉底白刺种子播种于装有纯净河沙的塑料盆（高14 cm，直径12 cm，底部有孔），5个塑料盆放在一个塑料桶里（高15 cm，直径80 cm）；然后放置在温室中培养，在（27± 2）℃ 下光照 (400～800 μmol·m-2·s-1)培养14 h和在（25±1）℃黑暗培养10 h，温室的相对湿度维持在60%～80%。待幼苗生长至2个月时，选取长势整齐的幼苗分为两组，第一组用200 mmol·L-1NaCl处理，第二组为对照，每组3个重复。7 d时分别取样，放置于-80℃冰箱中保存备用。

1.2 叶片总 RNA 样品的提取

取比拉底白刺幼苗叶片放于已灭过菌的研钵中，添加适量液氮快速研磨，然后取0.3 g 磨好的样品全部装入1mL EP 管中；加入1 mL Trizol 试剂，摇匀后静置10 min；加入0.3 mL氯仿，摇匀后静置10 min，然后离心15 min (4℃, 12 000 g) ；吸取上清液。再加入等体积异丙醇，静置4 h后离心15 min (4℃, 12 000 g) ；去上清液，并用75%酒精洗涤，然后同样在高速冷冻离心机上离心5 min(4℃, 12 000 g) ；用 75% 乙醇重复洗涤后，在冰上晾干；用无 RNAase 的水溶解后，分别用 1% RNA琼脂糖凝胶电泳和 Nanodrop 仪检测 RNA 样品的质量和浓度。

1.3 RNA-seq 测序

总RNA经DNase I消化后，通过富集含有polyA尾巴的mRNA、mRNA打断成200 nt的小片段、反转录成第一条链cDNA、合成第二链cDNA等步骤，用 QIA quick PCR extraction kit 进行纯化，然后进行末端修复加 polyA。经回收纯化、PCR 富集获得 cDNA 文库。文库构建完成后，使用 Agilent 2100 Bioanalyzer和 ABI StepOnePlus Real-Time PCR System 分别对文库的插入片段长度和有效浓度进行检测，样品合格后在 Iilumina 测序平台 (Illumina HiSeqTM2000) 进行高通量测序。将得到的序列原始文件 Raw Reads 进行过滤、比对、de novo组装、Unigene 功能注释、基因定量、差异基因筛选。其中，差异表达基因的筛选条件是FDR (False Discovery Rate) ＜0.001 且基因表达差异倍数 (fold change) 大于2倍，选取这些基因作为应答盐胁迫的差异基因，用于后续分析。

1.4 差异基因的功能注释与富集分析

1.4.1 GO 功能注释 利用 Blast2GO 软件将比对到 Gene Ontology 数据库中的所有差异基因进行GO 功能注释和分类[15]，网站是 http://www.geneontology.org/，获得转录本的GO 功能注释信息。

1.4.2 GO 显著性富集分析 差异表达基因比对Gene Ontology数据库后，进行映射，得到目的GO条目（term）。以这些GO term为单位，通过超几何检验（phyper），得到检验p值，Bonferroni校正后，对于correctedp-value＜0.05的GO term定义为在差异表达基因中显著富集的GO term。

1.4.3 KEGG 显著性富集分析 将差异表达基因应用KEGG数据库预测其参与的代谢途径。以correctedp-value≤0.05为标准，筛选在差异表达基因中显著富集的代谢途径。

1.4.4 基因互作网络构建 利用分子相互作用分析软件(STRING)和相互作用蛋白数据库（DIP）进行蛋白质相互作用网络预测，运用Cytoscape程序构建蛋白质相互作用网络[16]。

1.5 实时荧光定量 PCR

以分别提取的对照组及盐胁迫处理7d时的幼苗叶片RNA为模板，使用Takara公司的PrimeScriptTMRT reagent Kit with gDNA Eraser试剂盒反转录合成cDNA第一链作为实时荧光定量PCR(qRT-PCR) 反应的模板。从转录组数据中随机选择7个差异表达基因，利用 Primer 5.0 设计 qRTPCR 引物，如表1。以18S rRNA 作为内参，引物为 5 ′-GCTGGATTTGCTGGTGGTAT-3′ 和 5 ′-TTCCTGGGTCTGTGCCTGT-3 ′ (表 1)。 使 用SYBR® Premix Ex TaqTM荧光定量试剂盒进行荧光定量检测。每个样品进行3次重复，采用 2-ΔΔCt计算相对表达量。

表 1 qRT-PCR引物序列Table 1 Primers for qRT-PCR

2 结果

2.1 盐胁迫下比拉底白刺叶片转录组分析

以盐胁迫处理及对照组白刺叶片为材料，分别提取RNA后建立cDNA文库，用Illumina HiseqTM2000测序。通过de novo拼接、组装后共获得168 463条unigenes。对测序获得的所有unigenes分别与COG、GO、KEGG、KOG、Pfam、Swiss-Prot、NR七大功能数据库进行比对和功能注释，注释率分别为 21.5%、41.2%、14.2%、37.2%、42.4%、40.4%、59.7%（表2）。结果共注释到101 016条unigenes，占所有unigenes的60.0%，剩余67 447条（40.0%）unigenes未得到注释。在注释到的所有unigenes里，长度在300～1 000 bp的有24 329条，占24.1%；长度≥1 000 bp的有76 687条，占75.9%（表 2) 。

表 2 Unigene在多个数据库中注释结果统计Table 2 Statistics of annotation results in the databases

2.2 盐胁迫下比拉底白刺叶片差异表达基因

对盐胁迫处理及未处理对照的比拉底白刺幼苗叶片RNA-seq测序后得到的所有基因进行差异表达筛选，以log2fold change (基因差异表达倍数)的绝对值＞1且错配率FDR＜0.001为差异表达基因(DEG）的筛选标准， 盐胁迫7 d时，共筛选得到差异表达基因196个，其中，79条DEGs在盐胁迫响应中上调表达，下调表达DEGs有117条。基因表达差异变化绝大部分集中在3倍左右，其中，comp73009_c0_seq1差异变化为-7.10，下调变化最大，该unigene注释为5′-甲基硫腺苷/s-腺苷同型半胱氨酸核苷酶（5′-methylthioadenosine/S-adenosyl-homocysteine nucleosidase）；comp79737_c1_seq2和comp82550_c1_seq4差异变化分别为5.29和5.21，上调变化最高，然而前者尚未得到注释，功能未知，后者注释为赤霉素调节蛋白（Gibberellin-regulated protein）。推测以上显著差异unigenes可能参与比拉底白刺的耐盐过程，需要进一步研究。

2.3 盐胁迫下白刺差异表达基因的 GO 分析

Gene Ontology（GO）分析能够比较全面地描述差异表达基因的主要生物学功能，是一个国际标准化的基因功能分类系统。对盐胁迫处理7 d的比拉底白刺叶片DEGs进行GO分析，划分为64个功能组（图1），GO数据库中的结构根据功能组的不同又进一步分为三类描述：生物学过程、分子功能和细胞组分图1显示了比拉底白刺盐胁迫后差异基因表达谱的总体情况。生物学过程主要集中在代谢过程、单体过程、细胞过程和刺激反应。在分子功能中，注释为结合和催化活性的差异基因最多，其次为转运活性和抗氧化活性。在细胞组分大类中，差异基因被注释到最多的亚类依次为细胞、细胞组分、细胞器、膜、细胞器组分和膜部分（图1）。以上注释较多的过程涉及到的基因可能参与了比拉底白刺的盐胁迫响应过程，进一步挖掘这些差异表达基因有助于研究其抗逆机制和品质改良。

图 1 盐胁迫下比拉底白刺叶片差异表达基因的GO功能聚类注释图Fig. 1 Gene Ontology classification annotation of differentially expressed genes in the leaves of N. billardieri under salt stress

通过进一步对比拉底白刺叶片DEGs的GO富集分析（校正后p-value＜0.05）发现：在生物学过程中DEGs富集最显著的GO term分别是细胞壁高分子分解代谢过程、几丁质分解代谢过程、多糖分解代谢过程、氨基糖代谢过程、防御反应。在细胞组分中，富集最显著的GO term分别是液泡、叶绿体类囊体膜、细胞外区域、质外体、细胞壁。在分子功能中的DEGs富集最显著的GO term分别是几丁质结合、几丁质酶活性、氧化还原酶活性、醌结合、血红素结合。

2.4 盐胁迫下白刺差异表达基因的 KEGG 分析

在生物体中，基于通路的分析有助于进一步了解基因的生物学功能和上下游基因的相互作用，深入理解基因与功能的关系。本文使用KEGG功能富集分析盐胁迫处理7 d时比拉底白刺幼苗叶片的差异表达基因可能涉及的代谢通路，发现比拉底白刺响应盐胁迫的DEGs涉及到25条KEGG代谢通路（图2），其中，碳水化合物代谢过程的差异基因最多，占14%，其次为运输和分解代谢，占9%；蛋白质折叠、分类、降解和其它次级代谢产物的生物合成均占7%，免疫应答系统占6%，氨基酸代谢和信号转导分别占5%和4%。与盐胁迫应答相关的环境适应占2%（图2）。KEGG分析可以从功能的角度聚焦到通路及基因，也可以从基因的角度锁定功能和互作关系，直观地显示了比拉底白刺在盐胁迫下差异表达基因的代谢过程和信号通路，有助于更好地研究比拉底白刺响应盐胁迫的分子机制。

图 2 盐胁迫下比拉底白刺叶片差异表达基因的KEGG通路分析Fig. 2 KEGG pathway analysis of differentially expressed genes in the leaves of N. billardieri under salt stress

以校正后p-value＜0.05为筛选阈值，盐胁迫7 d比拉底白刺叶片差异表达基因的KEGG富集分析发现，有9条代谢途径产生显著变化（图3），分别是苯丙烷类生物合成、氨基糖和核苷酸糖代谢、托烷，哌啶和吡啶生物碱的生物合成、抗原加工和呈现、内质网中的蛋白质加工、核苷酸切除修复、戊糖和葡萄糖醛酸的相互转化、内吞作用、嘌呤代谢（图3）。以上9条代谢通路有助于获得比拉底白刺在抵御盐胁迫时的代谢信息，以便更清楚地理解其在盐胁迫应答中的过程。

2.5 盐胁迫下白刺差异表达基因的互作网络

为了研究比拉底白刺幼苗是如何通过基因之间的相互作用来传递盐胁迫信号，通过String数据库和搜索软件构建了盐应答基因在细胞内的互作网络（图4、表3）。在基因作用网络中检测到一组相互作用的基因簇，共包含32个基因（图4）。组成相互作用网络的这些基因主要与信号转导、氧化还原平衡、蛋白质合成、蛋白质折叠和组装、蛋白质降解、转录调控、次级代谢、细胞拯救/防御等途径相关。转录调控和氧化还原平衡在基因相互作用网络中占据着重要位置，它们可能在比拉底白刺幼苗响应盐胁迫反应中起重要作用。此外，基因 comp71223_c0_seq1、基因 comp83599_c0_seq4和基因comp71691_c0_seq2，作为节点，占据了网络的中心位置。

2.6 差异表达基因的 qRT-PCR 验证

为了验证RNA-Seq数据的准确性，随机选择7个差异表达基因，取对照和盐胁迫7 d的比拉底白刺叶片进行qRT-PCR实验。结果显示：这7个基因在qRT-PCR的结果与转录组分析的趋势一致（图5），而在表达变幅倍数上有些许差异，可能是由于两类实验的检测灵敏度及数据分析方法不同造成的。盐胁迫处理前后，7个差异基因的表达特征与测序结果呈现相同的变化趋势，验证了转录组测序结果的可靠性。

图 3 盐胁迫下比拉底白刺叶片差异表达基因的KEGG富集分析Fig. 3 KEGG Enrichment Analysis of differentially expressed genes in the leaves of N. billardieri under salt stress

图 4 盐胁迫诱导的差异表达基因相互作用网络.Fig. 4 The interaction network of differentially expressed genes induced by salt stress.

3 讨论

本研究进行了多次预实验，通过观察盐胁迫下植株的表型发现，盐生植物比拉底白刺在7 d时开始出现叶片轻微卷曲的现象，之后开始逐渐萎蔫，到9 d时萎蔫症状较明显。相对于非盐生植物，盐生植物能够耐受一定强度的盐胁迫环境一段时间，在一定盐浓度范围的土壤环境下正常生长[17-18]。为了研究比拉底白刺应答盐胁迫的分子机制，本实验选择在其表型性状刚开始出现时进行取样和转录组测序分析，以便更好地理解盐生植物特有的耐盐机制。

表 3 盐胁迫诱导的基因相互作用网络中差异表达基因的名称及注释Table 3 The name and annotation of differentially expressed genes responding to salt stress in gene interaction network analysis

图 5 差异表达基因的qRT-PCR验证Fig. 5 Verification of differently expressed genes using qRT-PCR

植物盐胁迫是一个由多个基因参与表达调控的作用网络，转录组提供了基因的多样性、表达水平差异以及时空变化，利用转录组测序技术的实用性和高效性可以促进植物响应非生物逆境胁迫的分子机制研究[19-20]。本研究采用Illumina HiseqTM2000高通量测序技术对盐胁迫7 d与未做处理的比拉底白刺幼苗叶片进行了转录组分析，共获得了168 463条unigenes，与七大数据库的注释率达60.0%，注释结果较好，大部分基因功能得到注释，与泡泡刺（N. sphaerocarpa）的注释率（62.1%）接近[21]。剩余67 447条（40.0%）unigenes未得到注释，可能是比拉底白刺中的新基因。

植物中各类转录因子对于调控各种诱导型基因的表达以及植物生长发育、适应环境等起主要作用[22]。该研究筛选出的差异基因中涉及到了WRKY转录因子（comp85183_c0_seq4、comp85183_c0_seq10、comp74597_c0_seq1）、乙烯响应转录因子（comp82253_c0_seq3）、NAC转录因子（comp-69459_c1_seq1）、GATA转录因子（comp74878_c0_seq4）、发病机制相关转录因子（comp83152_c0_seq5），这些unigenes在盐胁迫下均有显著响应，其中，WRKY转录因子具有高度保守的WRKY结构域，研究表明能够参与损伤、衰老、发育、抗病等抗逆反应[23]。乙烯响应转录因子（ERF）家族属于AP2/ERF超转录因子家族，具有保守的一个AP2/ERF结构域，在植物生长发育、环境应答反应（干旱、低温、高盐、病虫害等胁迫）中具有重要作用[24]。本研究发现的几个转录因子需要进一步挖掘其功能，以期为探究比拉底白刺的耐盐机理作出贡献。

通过盐胁迫诱导的差异表达基因GO分析结果可见，盐胁迫主要引起了比拉底白刺叶片内与代谢过程、环境刺激应答、防御反应等有关的生物学过程变化，以及与催化活性、结合活性、转运活性和抗氧化活性有关的分子功能的变化，这与比拉底白刺在蛋白质水平上对盐胁迫的应答模式有些类似[14]。初步表明，比拉底白刺在受到低盐胁迫时，可能通过合成或者降解一些蛋白酶类，改变蛋白质结构和含量，以抵抗低盐逆境。

从比拉底白刺幼苗叶片差异表达基因参与的代谢途径可见，盐胁迫主要引起了白刺幼苗叶片内与碳水化合物、氨基酸代谢相关的基因表达发生变化，同时也诱导了与蛋白质折叠、分类和降解相关的基因表达发生改变，这可能与盐胁迫导致了幼苗叶片内一些蛋白质代谢变化相关[4,14]。比拉底白刺幼苗叶片在受到低盐逆境胁迫时，氨基酸和核苷酸糖代谢途径有显著富集，可能其光合作用也受到一定影响，碳水化合物代谢活跃，推测其可能通过减少物质合成，控制能量代谢来抵抗低盐环境[25]。由盐胁迫引起的植株叶片表型的轻微变化可能是由于差异表达基因所参与的代谢途径的改变引起的，正是由于植株体内代谢途径的改变引起了诸如光合作用、质膜渗透力等生理现象的变化，进而导致表型的变化[4,14]。

将比拉底白刺幼苗叶片响应盐应答差异基因在细胞内构建互作网络，作为网络中的“hubs”（连接到许多其他的蛋白质）和“bottlenecks”（在一个网络的子网络的关键连接器），如基因comp71223_c0_seq1注释为热激同源蛋白（Heat shock cognate 70 kDa protein），基因 comp83599_c0_seq4注释为L型凝集素域受体激酶（L-type lectin-domain containing receptor kinase IV.1），和基因comp71691_c0_seq2注释为Win类蛋白（Win-like protein），它们占据了网络的中心位置，可能在白刺幼苗叶片对盐胁迫处理的响应过程中发挥重要作用。因此，通过基因相互作用网络分析可以初步了解这些基因在比拉底白刺响应盐胁迫过程中起的作用，互作关系有待进一步研究。

4 结论

本研究通过对盐胁迫 7 d 与未处理的比拉底白刺叶片进行转录组测序及相关生物信息学分析，共获得168 463条 unigenes ，筛选到196条差异表达基因。通过GO 功能和 KEGG 通路分析，理清了差异表达基因富集的分子功能与代谢通路。此外，进一步构建了差异基因互作网络，得到了重要的节点基因。综上所述，本研究为挖掘和找寻重要耐盐候选基因提供了参考，同时为研究比拉底白刺响应盐胁迫的分子机制，以及为下一步培育比拉底白刺耐盐新品种奠定了基础。