芦竹茎对盐胁迫响应的转录组分析

孙源长，罗 琳，林 辉，林冬梅，林占熺

(福建农林大学国家菌草工程技术研究中心，福州 350002)

【研究意义】随着全球人口的急剧增长以及科学技术的快速发展，人们对能源的需求量逐渐增大。石油和煤炭等不可再生能源日益减少，因此迫切需要开发利用可再生能源，而生物质能是其重要组成部分。目前,一些作物和生物量较大的植物作为生物能源被广泛应用，但这两者竞争生长，不适合一起种植在耕地中，将适应能力强、生物量大的植物种植到盐渍地等生长环境恶劣的区域是很好的选择。最近，研究表明全球超过10×108hm2的土地受盐分影响，约占地球陆地表面的7%，并且还在持续蔓延。近期研究结果表明，伊朗西北部水土流失和盐渍化问题严重[1]，而亚洲(中国、印度、巴基斯坦)、美洲(美国)、非洲(南非)等地的土壤盐含量也严重超标[2]。土壤盐分是影响植物生长和发育最严重的非生物胁迫之一，盐度主要以两种方式影响植物的生长发育。渗透胁迫是在高浓度盐环境下根表面的水势降低，导致植物吸水量减少，扰乱植物的正常生命过程；离子胁迫是植物细胞内 Na+/ K+的失衡引起的，许多生化反应将不能进行[3-4]。芦竹又称“荻芦竹”，为禾本科芦竹亚科芦竹属的多年生草本植物，已在亚洲、南欧和北非等地种植上千年，能够适应各种不同环境，在温带地区可以长期在寒冷、盐碱、干旱和涝渍等条件下生存而不会造成重大损失[5-7]。这种植物能够耐受盐胁迫，在高盐度下仍保持较高的光合作用和生物量，不会对光化学反应以及叶片结构和完整性造成较大损伤[8]。因此，探究芦竹茎的盐响应方式，鉴定出参与盐响应的重要基因，对未来分子育种，提高植物耐盐性，减轻能源压力具有重大意义。【前人研究进展】动物可以移动迁徙来应对即将和已经发生的胁迫，而植物却只能通过调节自身状态来应对环境的变化[9]。自身状态的调节主要是通过基因的表达调控来实现，并且这种表达调控是动态的。为应对各种生物和非生物胁迫，植物在长期进化过程中产生了多种自我保护和防御机制。例如当环境发生剧烈变化时，植物会调整一系列信号通路，来提高自身生存率[10-11]；通过关闭气孔来调节叶片的光合过程[10]；通过调节体内抗氧化系统消除多余的活性氧，维持体内氧化还原平衡[10-11]；通过细胞渗透调节物质改变细胞内水势，维持一定的渗透压，从而保护植物在盐胁迫下的持续生长[10]；植物激素通过协同或拮抗作用来调节自身[12]；转录因子(TF)与启动子区域中顺式作用元件发生特异性结合，进而调控靶基因的表达[13]。同样，芦竹主要也是通过这些方式来减轻外界不利环境对自身的影响。【本研究切入点】目前关于芦竹盐响应的相关研究并不多，所以迫切需要对其进行深入研究以丰富芦竹的耐盐基因资源，为未来的分子育种奠定基础。RNA测序 (RNA-Seq) 一直是研究应激耐受分子机制的有效方法，被广泛应用于各个领域[14-18]。【拟解决的关键问题】本研究以不同浓度NaCl盐溶液处理的“绿洲1号”芦竹的茎为材料，利用RNA-seq手段分析，探究芦竹茎的盐响应方式，鉴定出参与盐响应的重要基因，进一步加深对芦竹耐盐分子机制的理解。

1 材料与方法

1.1 供试材料

从福建农林大学国家菌草工程技术研究中心选取生长6个月并且长势一致的“绿洲1号”芦竹(原产地：莱索托)组培幼苗，之后移栽到温室(在30 ℃/25 ℃以14 h光照/10 h黑暗的光周期生长，湿度为60%，光照强度为350 μmol/(m2·s)并用1/5霍格兰营养液水培1个月进行缓苗，之后分别用0、50、100、150和200 mmol/L NaCl盐浓度溶液处理，处理时间为36 h。

1.2 RNA提取与转录组测序

分别收集不同NaCl盐浓度溶液处理后的“绿洲1号”茎，每个重复由2～4株芦竹自下往上数第1～8节混合获取，共3次重复。样品用液氮研磨成粉末状混匀。

称取样品粉末0.2 g，用Trizol试剂提取茎的总RNA。取2 μL RNA在1%琼脂糖凝胶上初步评估RNA降解程度及污染情况。通过Nano Drop 2000®分光光度计测定RNA的纯度与浓度。使用Agilent2100(Agilent Technologies)检测RNA完整性。对合格的样品委托诺禾致源公司进行转录组测序。简要流程如下：用带有Oligo(dT)的磁珠富集mRNA，随后加入fragmentation buffer将mRNA打断成短片段，以mRNA为模板，合成双链cDNA。之后进行末端修复，加A尾并连接测序接头。最后进行PCR扩增，并用AMPure XP beads纯化PCR产物，得到最终文库。对合格的样本在Illumina Novaseq 6000平台上进行双端测序，读长为150 bp。

表1 用于qPCR反应的引物

1.3 转录组数据的质量控制与组装

为了得到更加完整的转录本，此次组装额外使用了本实验室同一品种的15个叶片转录组数据。原始数据(Raw reads)由FastQC(https://www.bioinformatics.babraham.ac.uk/projects/fastqc/)软件进行质量检测。任何低质量的原始数据和接头序列均通过Trimmomatic[19]软件过滤，最终得到高质量数据(Clean data)。使用Trinity-v2.13.2[20]软件对转录组数据进行从头组装，利用CD-HIT[21]软件去冗余(相似性设置为0.95)得到最终单基因簇(Unigenes)序列。

1.4 基因表达定量和差异表达分析

使用Bowtie2和RSEM对基因进行表达定量[22-23]。在R中使用DESeq2[24]进行差异表达基因(Differentially expressed genes，DEGs)的检测。DEGs的过滤标准是|log2FoldChange|≥1并且P≤ 0.05。

1.5 GO和KEGG富集分析

使用clusterProfiler[25]R包对DEGs进行GO与KEGG富集分析。参数设置为校正后的P<0.05, Unigene的注释来源为eggNOG5.0[26]数据库。

1.6 转录因子鉴定

所有差异表达基因均通过iTAK(v1.6)[27](http://itak.feilab.net/cgi-bin/itak/online_itak.cgi)在线程序进行转录因子的预测，比对类型为核酸序列。

1.7 基因的qRT-PCR验证

随机选取5个差异表达基因验证转录组数据是否可靠。根据Michele等[28]的实验结果，选取较为稳定的RPN6基因作为内参基因，所有基因通过Primer 5.0软件设计引物，qPCR引物见表1。使用TransScript®Uni All-in-One First-Strand cDNA Synthesis SuperMix for qPCR(One-Step gDNA Removal)试剂盒进行反转录。qPCR反应根据PerfectStart®Green qPCR SuperMIX试剂盒说明书进行。采用2-△△Ct法计算基因的相对表达量[29]，用GraphPad Prism 8 软件进行图像绘制。

2 结果与分析

2.1 测序数据的组装结果

使用Trinity[20]软件对Clean data进行De novo组装，共产生971 309个转录本，GC含量为47.14%。对组装后的结果文件用CD-HIT[21]软件去冗余，得到最终的Unigenes。Transcript contigs的N30达到2053 bp, N50达到1228 bp，contig 平均长度为823.47 bp；Unigene contigs的N30达到1915 bp, N50达到1126 bp，contig平均长度为769.52 bp。转录本和基因的长度数量分布统计情况见图1。以上结果表明此次转录本组装效果优良，可以进行后续分析。

图1 转录本和单基因簇长度数量分布统计Fig.1 Length and number distribution statistics of transcripts and unigenes

2.2 不同盐浓度胁迫下的基因表达量化和差异表达分析

在R中通过 DESeq2包分析“绿洲1号”响应盐胁迫的差异表达基因。在叶片中检测到的差异表达基因为7305个，其数量分布情况见图2。随着盐浓度的增加，上调差异表达基因、下调差异表达基因与总的差异表达基因的数量都逐渐增加。上调的差异表达基因数目始终大于下调的差异表达基因数目。以上结果表明“绿洲1号”茎在响应盐胁迫时，随着盐浓度的增加其调控过程越来越复杂，并且主要通过基因的上调起作用。

图2 不同盐浓度下茎差异表达基因的数目Fig.2 Number of differentially expressed genes in stems under different salt concentrations

由图3可知，在茎中有104个基因是不同浓度盐处理下的差异表达基因所共有的，将它们作为芦竹盐响应核心差异表达基因。随着盐浓度的增加，各个组别中特有的差异表达基因数目也逐渐增加。stem_CK_vs_50 mmol/L到stem_CK_vs_100 mmol/L、stem_CK_vs_150 mmol/L、stem_CK_vs_200 mmol/L的特异性差异表达基因数目分别从111个上升到292、1141、3347个。这可能是由于随着胁迫程度的加深，越来越多特定的信号传导机制被激活。以上结果表明茎存在一些较为保守的基因来响应盐胁迫，并且随着盐浓度的增加越来越多的特异性差异表达基因参与到表达调控中来。

图3 不同盐浓度下茎差异表达基因的Venn图Fig.3 Venn plot of differentially expressed genes in stems under different salt concentrations

2.3 差异表达基因的GO和KEGG富集分析

由图4可知，在芦竹茎响应盐胁迫过程中，50 mmol/L 盐浓度下差异表达基因主要富集到钾离子稳态、硝酸盐响应、细胞壁修饰、水杨酸合成相关条目。100 mmol/L 盐浓度下，臭氧、泛素化、渗透压调节相关条目富集到的基因数目最多。150 mmol/L盐浓度下富集到对硝酸盐、无机阴离子、光形态、端粒酶活性等条目。与150 mmol/L盐浓度下富集结果类似，200 mmol/L盐浓度下除了富集到无机阴离子转运条目还富集到与水响应、光反应调节相关的条目。

图4 茎中差异表达基因的GO富集Fig.4 GO enrichment of differentially expressed genes in stems

由图5可知，在茎中，50 mmol/L盐浓度只富集到了单萜生物合成途径；而100、150和200 mmol/L盐浓度富集到氮代谢、昼夜节律、萜类合成、玉米素合成和维生素B6代谢等通路。

图5 茎中差异表达基因的KEGG富集Fig.5 KEGG enrichment of differentially expressed genes in stems

2.4 盐响应核心差异表达基因分析

由图6可知，在芦竹茎中大部分盐响应核心差异表达基因上调，只有少数基因下调；并且与下调基因相比上调基因的差异倍数更大。说明芦竹响应盐胁迫主要是通过基因上调起作用,并且上调基因的反应更加强烈。在茎中88315_c0_g1_i1、1200_c2_g1_i1、13330_c0_g1_i8、46571_c0_g1_i7、3434_c0_g1_i11和713_c0_g1_i32的差异倍数最大。这些基因在“绿洲1号”茎基础盐响应过程中发挥极其重要的作用，注释的具体情况见表2。

图6 茎中核心差异表达基因的差异倍数分析Fig.6 Fold analysis of core genes expression differences in stems based on log2FoldChange

表2 茎中重要的盐响应核心差异表达基因

2.5 转录因子家族鉴定

对茎参与盐响应表达调控的差异表达基因进行转录因子预测，一共预测到41个转录因子基因家族，家族成员数目最多的是AP2/ERF-ERF家族，共94个；其次是MYB-related家族和WRKY家族，家族成员分别为65、57个(图7)。

图7 茎中差异表达基因转录因子预测数目Fig.7 Predicted numbers of transcription factors of differentially expressed genes in stems

为探究AP2/ERF-ERF、MYB-related和WRKY转录因子家族的表达模式，对其基因表达差异倍数进行分析，AP2/ERF-ERF的log2FoldChange在-1.85～10.88；MYB-related的log2FoldChange在-3.00～7.85；WRKY的log2FoldChange在-1.58～6.13。并且这3类转录因子在响应盐胁迫时大部分呈上调状态，说明其主要通过增加基因的表达起作用(图8)。

A：AP2/ERF-ERF；B：MYB-related；C：WRKY图8 转录因子家族的差异倍数分析Fig.8 Fold analysis of transcription factor families in stems based on log2FoldChange

对前8个差异表达倍数较大的转录因子进行分析(表3)，AP2/ERF-ERF、MYB-related和WRKY转录因子家族差异倍数最大的转录因子分别为：18947_c0_g1_i3、7562_c1_g1_i2和24870_c0_g1_i2，平均log2FoldChange达到7.88、5.81和5.28,它们在芦竹响应盐胁迫时可能发挥着重要作用。

表3 差异表达倍数较大的转录因子

2.6 实时荧光定量PCR验证

随机挑选5个差异表达基因对转录组数据进行qRT-PCR验证，可以看出这些基因的相对表达量和转录数据的基因表达模式基本一致(图9)，也说明转录组数据较为可靠。

左侧纵坐标为qPCR相对表达量，右侧纵坐标为TPM值。A：13361_c0_g1_i15；B: 14873_c0_g1_i6；C: 70740_c0_g1_i2；D: 129260_c0_g2_i1；E：122017_c0_g1_i1The left ordinate is the relative expression level of qPCR, and the right ordinate is the TPM value. A: 13361_c0_g1_i15; B: 14873_c0_g1_i6; C: 70740_c0_g1_i2; D: 129260_c0_g2_i1; E: 122017_c0_g1_i1图9 差异表达基因qRT-PCR 验证Fig.9 qRT-PCR analysis of differentially expressed genes

3 讨 论

如果为植物提供如光照、温度和水分等适宜的生长环境，植物就会相对高产，反之有可能导致植物发育不良甚至死亡。而产量与植物发育过程中的非生物胁迫条件呈负相关[30]。在全球范围内，由于非生物胁迫(干旱、盐渍、洪水、霜冻、营养缺乏等因素)导致的损失估计为每年1.2×1011美元[31](http://www.fao.org/docrep/008/y5800e/Y5800E06.htm)。盐胁迫是降低作物产量最重要的非生物胁迫之一，因此，对植物进行盐胁迫的相关研究将有利于筛选出耐盐的关键基因，对提高植物的耐盐性和产量具有重要意义。芦竹具有较强的抗逆性，能够在高盐环境下生长。本实验以“绿洲1号”芦竹的茎为材料，使用不同浓度的NaCl盐溶液进行处理，利用RNA-seq手段分析芦竹茎的盐响应方式，鉴定出重要的盐响应相关基因，为以后的遗传改良提供指导。

3.1 芦竹盐响应不同盐环境的转录差异

本研究共检测到7305个盐胁迫差异表达基因，随着盐溶液浓度的增加,差异表达基因的数目也相继增加，并且越来越多的特异性差异表达基因参与到表达调控中来。以上结果表明不同浓度的盐环境均可诱导大量的基因来响应盐胁迫，随着浓度的增加表达调控的复杂程度也可能随之上升。

对差异表达基因进行GO与KEGG富集分析，与代谢产物相关的条目主要富集到与萜类合成、玉米素合成、亚油酸代谢和水杨酸生物合成途径，这些代谢物质作为调节剂在植物应激过程中发挥着重要作用[32-34]，对提高芦竹的盐耐受性有重要作用。无论从低盐还是高盐环境都富集到与萜类合成相关的条目，可以推测芦竹能通过快速调节萜类化合物合成来应对盐环境的变化。富集到钾离子稳态、与细胞壁相关[35]、调节对渗透压的反应、对盐胁迫反应的正调节、应对缺水的正调控、无机阴离子跨膜转运和磷酸根离子跨膜转运等条目，将对芦竹感应盐胁迫、维持细胞的渗透压稳态和调节相关的信号通路有所帮助。在盐浓度达到200 mmol/L时，即较高盐环境下，还富集到了维生素B6代谢途径。维生素B6近些年来被广泛研究，它是多种代谢酶的必需辅酶，参与植物的生长发育、生物胁迫和非生物胁迫等一系列过程。它还是一种有效的抗氧化剂，研究表明维生素B6在植物体内抗氧化防御中具有作用，在强光下限制1O2积累并防止1O2介导的氧化损伤[36]。

3.2 芦竹盐响应过程中的关键基因

共筛选到104个盐响应的核心差异表达基因，其中大部分基因上调，只有少数基因下调，说明这些基因主要是通过增加表达量来响应盐胁迫。对其中差异倍数最大的基因进行注释，其中GATL1编码糖基转移酶，有助于木聚糖的生物合成。此外，植物糖基转移酶能够将糖添加到小分子中，例如植物激素[37]、水杨酸[38]和类黄酮[39]。而木聚糖与细胞壁合成相关，细胞壁参与了盐分胁迫的感知和响应，对于细胞的压力保护和反应都至关重要；CRF4是ERF/AP2 转录因子家族，之前的研究表明它在响应非生物胁迫中起着重要作用[40]；SF3B5A编码剪接因子3B亚基5/RDS3复合亚基，在选择性剪接中起作用；SYR1为系统素受体，主要位于质膜上，有助于植物防御过程中系统素介导信号传导的调节[41]。

3.3 芦竹盐响应过程中起重要作用的转录因子

植物编码大量转录因子，它们主要按照 DNA 结合结构域分类[42]。越来越多的证据表明，转录因子可调节多种生物过程，例如各种生物及非生物胁迫、光和压力信号、种子成熟和花发育等[43]。目前已经鉴定到了一系列的转录因子家族在盐胁迫过程中发挥着重要作用，例如bZIP、AP2/ERF-ERF、MYB、WRKY、NAC和bHLH等[44]。

为筛选出在芦竹盐响应过程中起重要作用的转录因子，本实验对差异基因进行研究，共注释到41个转录因子基因家族，推测它们可能在芦竹对盐胁迫的反应中起重要作用，其中成员数目最多的包括AP2/ERF-ERF、MYB-related和WRKY家族。AP2/ERF-ERF家族是植物所特有的一类基因家族，该家族中的每一个成员都具有单个或2个AP2结构域，主要参与植物的生长发育、细胞周期和生物及非生物胁迫过程，广泛参与脱落酸、赤霉素和细胞分裂素等信号通路。MYB-related家族是植物最大的转录因子家族之一，根据 DNA 结合结构域相邻重复的数量可分为四大类，在植物的发育和胁迫反应中起重要作用[45]。WRKY家族是调节植物生物过程信号网络的重要组成部分，通常作为抑制因子和激活因子发挥作用，在调节植物的许多胁迫反应中起关键作用。例如，在水稻中，热休克诱导型HSP101启动子驱动的OsWRKY11过表达将会导致耐热和耐旱性增强[46]。同样，OsWRKY45的过表达除了提高抗病性外，还导致耐盐性和耐旱性增强[47]。此外，NAC、bZIP、bHLH等转录因子家族也在本研究中被发现，这与前人的研究一致，并且这些转录因子绝大部分通过上调表达量来响应盐胁迫。以上结果表明AP2/ERF-ERF、MYB-related和WRKY转录因子家族在“绿洲1号”芦竹响应盐胁迫的过程中可能具有重要的调控功能。

4 结 论

对0、50、100、150和200 mmol/L NaCl盐溶液进行处理的“绿洲1号”芦竹的茎进行转录组测序与分析，共筛选到10个差异表达倍数较大的基因和一些差异表达倍数较大的转录因子，它们可能在芦竹响应盐胁迫的过程中发挥着重要作用。本研究丰富了芦竹的耐盐基因资源，为分子育种提供了理想的改造靶点。