RNA-seq揭示碱性盐(NaHCO3)对细叶百合鳞茎基因表达的影响

李冠嵘，何 好，朱国庆，陈诗雅，徐 阳，金淑梅

（东北盐碱植被恢复与重建教育部重点实验室，东北林业大学生命科学学院，哈尔滨150040）

0 引言

RNA-Seq转录组测序是分析基因表达水平、筛选差异表达基因、挖掘功能基因的重要手段，很多高等植物都进行了耐盐相关的转录组方面的研究。对耐盐棉花进行转录组测序，分析了与盐胁迫相关的基因[1]；通过转录组测序揭示了盐胁迫下柑橘根中基因表达的情况[2]；对在盐碱胁迫下高粱的转录组测序鉴定了与ABA相关的差异表达基因，推测了参与ABA信号传导的家族成员[3]；对烟草的转录组测序确定了分别响应NaCl和NaHCO3差异表达的miRNA，并确定了在NaCl和NaHCO3处理下的mRNA-miRNA相互作用[4]；通过转录组分析3种不同表型的陆地棉盐碱胁迫处理后反应的分子机制[5]。

百合的转录组学主要是在发育方面，如利用微阵列技术研究了麝香百合生殖细胞中基因的表达情况[6]，通过Illumina HiSeqTM2000测序平台研究了百合花粉萌发过程基因的表达[7]，在花色方面的研究主要有，为了研究与百合花颜色合成有关的调控信息，研究了东方百合和亚洲百合春花过程中的分子响应机制[8-10]；通过转录组测序，揭示了花色苷生物合成和运输所需表达的基因[11]；Li等[12]根据百合转录组测序结果研究了3种植物激素的信号传导途径。利用Roche454测序平台获得了4种基因型百合叶片的转录组测序信息[13]，在逆境方面也有研究，建立了卷丹冷胁迫盐分、激素压力下的转录谱[14-15]。

近几年，对细叶百合的转录组也有研究，对低温处理0天、30天、60天、90天后的细叶百合鳞茎进行转录组测序，探究鳞茎休眠解除的分子机制[16]。对参与细叶百合休眠解除过程中的WRKY差异表达基因进行筛选，详细解析了LpWRKY20基因在细叶百合休眠解除过程中的调控作用[17]。克隆得到MADS-box家族中的LpSVP阐明LpSVP基因的作用机制[18]。筛选得到21个 NAC基 因 ，并 得 到 含 LpNAC5、LpNAC13和LpNAC21基因的阳性植株进行研究[19]。

细叶百合(Lilium pumilum DC.)是百合科百合属的多年生草本植物，花美姿艳。不仅是重要的观赏性花卉，也具有很高的食用、药用等经济价值。大多数野生百合植株娇小细弱，其抗逆性不强，这极大地限制了百合在土壤盐碱化问题严重的地区的种植和应用[20]。细叶百合是在盐碱化土壤中生存的少数几种百合之一，是百合的优良种质资源之一，具有抗寒、耐旱、耐盐碱等优良的特性。细叶百合采自于盐碱土壤比较严重的黑龙江省大庆市地区，对于能在恶劣的盐碱条件下生长，其基因的表达有什么特殊之处呢？为了探究这个问题，对细叶百合在碳酸盐逆境下的转录组数据进行了详细的分析，希望揭示细叶百合在盐碱逆境中的基因表达方式，为合理利用盐碱地和大面积种植百合提供理论基础。

1 材料与方法

1.1 试验材料与处理

以培养在组培瓶中的一年生细叶百合为材料，将长势一致的细叶百合幼苗移植到1/2MS培养基(CK)和含有20 mmol/L NaHCO3（处理组）的MS培养基中，处理48 h后，将对照组(CK)和处理组(NaHCO3)各两组分别取样，液氮速冻后置于-80℃保存。试验在东北林业大学东北盐碱植被恢复与重建教育部重点实验室

1.2 细叶百合叶片总RNA的提取、文库构建、Illumina测序及数据质量控制

2次重复的对照组和处理组的细叶百合鳞茎cDNA文库的构建以及测序实验由苏州金唯智生物科技有限公司完成。具体实验操作流程如下：（1）提取到总RNA，mRNA的富集、cDNA链的合成。（2）纯化、洗脱、末端修复后进行PCR扩增。（3）使用Illumina HiSeq高通量测序平台对cDNA文库进行测序，对测序结果原始图像数据利用软件bcl2fastq进行图像碱基识别，初步质量分析后得到原始测序数据。分析测序数据进行质量评估，使用二代测序数据质量统计软件Cutadapt对测序原始数据去除接头以及低质量序列，得到后续信息分析用的Clean Data。

1.3 测序数据组装及基因功能注释

采用组装软件Trinity对样品数据从头组装，组装结果通过序列聚类做进一步序列拼接和去冗余处理，得到长的非冗余的unigene序列并与Nr、Swiss-Prot、GO、KEGG、COG数据库进行比对，获得mRNA的注释信息。

1.4 序列比对以及差异表达基因(DEG)的鉴定

使用bowtie2软件进行短reads的比对，将Clean data比对Unigene上，进行分析，并RSEM软件使用FPKM方法计算基因的表达量。采用Bioconductor软件包的DESeq2进行样品组间的差异表达检测分析，以差异基因表达变化2倍以上且FDR≤0.05为标准进行筛选，得到最终的差异表达基因，选取9个差异表达基因进行QPCR验证。

1.5 实时荧光定量PCR验证

为了验证RNA-Seq中表达量数据的准确性，对细叶百合进行20 mmol/L NaHCO3处理24 h后提取细叶百合鳞茎RNA，将RNA反转录成cDNA作为模板，选取了9个（5个上调基因和4个下调基因）与盐碱胁迫密切相关的基因(BEL1、ERD6、GRP、MCM、FRO8、CPC、PPT、LBP、PEX7)设计荧光定量PCR引物（表1）进行qRT-PCR验证。UltraSYBR Mixture试剂进行实时荧光定量qRT-PCR检测，反应体系为20 μL（2×SYBR Green Mix 10 μL；10 μmol/L 上下游引物各0.5 μL；cDNA 模板 1 μL；ddH2O 8 μL）；反应条件为：95℃预变性10 min，95℃ 30 s，55℃ 30 s，72℃ 1 min，40个循环。得到细叶百合测序的转录组中上调和下调基因的表达量并进行分析。

表1 qRT-PCR所用引物序列

2 结果与分析

2.1 测序结果统计分析

对4个(CK1、CK2、Salt1、Salt2)cDNA文库进行的测序，通过质量控制后，共得到23.4 Gb Clean Data，各样品Q30碱基百分比均不小于89.35%。

Clean Data经组装后得到115061条unigene，N50长度为939 nt，115061条unigene中有65.18%的长度在200～500之间，具体信息见表3与图1。

图1 Unigene不同长度的占比情况

表3 经组装后的unigene情况

2.2 基因功能注释结果

将Unigene序列与Nr、COG、Swissprot、KEGG数据库进行Blast比对后，共获得56828个mRN的注释信息。55433个Unigene被注释到Nr数据库中，26973个Unigene被注释到COG数据库中，23610Unigene被注释到GO数据库中，37434个Unigene被注释到Swissprot数据库中，13142个Unigene被注释到KEGG数据库中。

表2 样品测序数据统计表

2.3 Nr注释结果

共55433个Unigene被注释到Nr数据库中，其中，注释到油棕(Elaeis guineensis)的Unigene最多，占总数的25%，其余注释到海枣(Phoenix dectylifera)、小果野蕉(Musa acuminatasubsp.malaccensis)、凤梨(Ananas comosus)、莲花(Nelumbo nucitera)、葡萄(Vitis vinifera)的Unigene分别占总数的19%、9%、5%、2%、2%。具体信息见图2。

图2 Nr注释分类图

2.4 KOG注释结果

共有26973条Unigene被注释到KOG数据库的26个分类中，其中数量最多的一类是一般功能类(R:General function prediction only，共3478条，12.89%)，其次是信号转导机制（T:Signal transduction mechanisms，共 3279条，12.16%），翻译后修饰（O:Posttranslational modificadification，共 3218 条 ，11.93%），具体信息见图3。

图3 KOG注释结果

2.5 GO注释结果

共有23610条Unigene被注释到GO数据库的3个大类与53个小类中，在分子功能(molecular function)大类中，注释到催化活性(catalytic)与结合绑定(binding)分类的数量最多；在细胞组分(Cellular Component)大类中，注释到细胞部分(cell part)与细胞器(organelle)分类的数量最多；在生物学过程(Biological process)大类中，注释到代谢过程(metabolic process)及细胞生理过程(cellular process)中的Unigene数量最多。具体信息见图4。

图4 GO注释结果

2.6 KEGG注释结果

共13142条Unigene被注释到KEGG数据库的五大类中，其中注释到代谢(Metabolism)这一大类的Unigene数量最多，占总数的38%，其余注释到的大类中，有机体系(Organismal Systems)占20%，遗传信息处理(Genetic Information Processing)占17%，环境信息处理(Environmental Information Processing)占14%，细胞过程(Cellular Processes)占11%。具体信息见图5。

图5 KEGG注释结果

2.7 DEG的分析

经过NaHCO3胁迫处理细叶百合植株24 h后，处理组和对照组鳞茎中的基因表达对比，共鉴定出390个差异表达基因，其中有155个基因上调，235个基因下调，具体信息见图6。

图6 差异表达基因DEG统计

2.8 荧光定量PCR验证结果

qRT-PCR的结果显示见图7，BEL1基因与对照相比表达量降低了90%，ERD6基因与对照相比表达量降低了85%，GRP基因与对照相比表达量降低了43%，MCM基因与对照相比表达量降低了38%，PPT基因与对照相比表达量增加了47%，LBP基因与对照相比表达量增加了48%，FRO8基因与对照相比表达量增加了53%，CPC基因与对照相比表达量增加了176%，PEX7基因与对照相比表达量增加了191%。所选取的9个基因的表达结果与转录组的结果基本趋于一致，表明RNA-Seq的结果准确可信。

图7 荧光定量PCR验证转录组可靠性

3 讨论与结论

百合的转录组在发育[6-7]和花色[8-11]方面研究比较多，对激素信号途径的百合转录组也有研究[11-12]，在逆境方面研究了卷丹在冷胁迫下的基因表达，并且筛选出与盐碱逆境有关的基因[14-15]，在低温处理下研究了细叶百合鳞茎转录组[16]，本研究同样是以细叶百合的鳞茎做为植物材料进行的转录组方面的研究，但对于盐碱胁迫下百合尤其是碱性盐胁迫下的基因转录组研究未见报道，本文的创新之处是在碱性盐(NaHCO3)处理后转录组的研究，通过数据分析，力图找到处理组和对照组中差异表达的基因，观察差异基因的表达，为了探究细叶百合耐盐碱的分子机理提供理论基础。

在本研究转录组的差异数据中除了一些常见的与盐碱有关的基因（脯氨酸转运蛋白，过氧化物酶体生物合成酶）外，还出现几个与糖转运（蔗糖转运蛋白）及金属运输（铁还原氧化酶，蓝铜蛋白）等基因，选取了9个与基因进行qRT-PCR验证，BELL基因是参与调控植物的发育过程的转录因子，3% NaCl处理条件下一些基因的表达水平明显上调，尤其NtBELL6基因的表达水平上调了3倍[21]，结果显示BEL1基因与对照相比表达量降低了90%，BELL1的表达在NaHCO3处理后下调，是否在其他NaHCO3的浓度下，该基因的表达会有所不同，需要进一步的研究。ERD6蛋白属于蔗糖转运蛋白，发现AtERD6能够快速地响应低温诱导，在低温处理1 h后表达显著上调[22]，而茶树中的CsERD6-7冷驯化过程中的表达显著被抑制[23]。在转录组测序和qPCR验证中，ERD6基因与对照相比表达量降低了85%，在NaHCO3处理下，基因的表达量下调，其作用机理还有待进一步的研究。GRP富含甘氨酸的RNA结合蛋白，盐芥叶中的TsGRP表达量在300 mmol/L NaCl处理后表达下调，根中的表达量上调[24]，在转录组测序和qPCR验证中，GRP基因与对照相比表达量降低了43%。是否是因为所选的胁迫浓度或者选取的器官不同，基因的表达量有所不同。还有待进一步研究。MCM编码真核生物的转录因子家族，水分胁迫下大豆籽粒中发现被富集到DNA复制通路上的差异表达基因MCM[25]本转录组中MCM基因与对照相比表达量降低了38%，表达量降低的原因还有待进一步研究。

铁还原氧化酶(ferric reduction oxidase，FRO)基因家族参与植物中广泛存在的各种生物过程，AtFRO8为线粒体FRO家族基因成员，在不同发育阶段参与植株体内Fe3+螯合物的还原[26]，FRO8基因与对照相比表达量增加了13%；CPC是一种蓝铜蛋白，影响生物的生理发育活动。拟南芥蓝铜蛋白能催化离子的氧化，参与了在细胞膜区域的电子转移反应来提高植物的逆境胁迫能力[27]本研究中CPC基因与对照相比表达量增加了22%；脯氨酸转运蛋白（ProT或PPT）是典型的Na+依赖型亚氨基酸转运蛋白类，在NaCl的胁迫下，马铃薯脯氨酸转运蛋白2个基因在2 h就开始急剧响应，其表达量均是对照的19倍以上[28]，本转录组所得到的数据与文献中的结果相一致，PPT基因与对照相比表达量增加了47%；BIP为热激蛋白家族成员，是内质网的一个主要的分子伴侣，大豆BIP能增强抗旱能力，减弱由干旱诱导的叶片衰老[29]本研究中BIP基因与对照相比表达量增加了48%，过氧化物酶体生物合成酶(PEX)是过氧化物酶体在生物体内合成的关键酶，胡杨PePEX11能够提高拟南芥在盐胁迫下的抗氧化能力，提升耐盐能力[30]PEX7基因与对照相比表达量增加了59%。这些在转录组中表达量上调的与文献中基因在逆境处理后，表达量上调的结果是一致的，说明细叶百合，BELL、ERD、GRP、FRO、CPC、PPT、LBP、PEX7基因与盐有一定的应答关系，为以后进一步深入研究这些基因的功能奠定一定的基础。

通过对细叶百合鳞茎进行转录组测序，共获得23.4 Gb的数据，通过组装得到115061条unigene。将Unigene序列与Nr、COG、Swissprot、KEGG数据库进行Blast比对后，共获得56828个mRNA的注释信息。其中，55433个Unigene被注释到Nr数据库中，其中，注释到油棕的Unigene最多，占总数的25%，其余注释到海枣、小果野蕉、凤梨、莲花、葡萄的Unigene分别占总数的19%、9%、5%、2%、2%。26973条Unigene被注释到KOG数据库的26个分类中，其中数量最多的一类是一般功能类，其次是信号转导机制。有23610条Unigene被注释到GO数据库的3个大类与53个小类中，在分子功能(molecular function)大类中，注释到催化活性与结合绑定分类的数量最多；在细胞组分(Cellular Component)大类中，注释到细胞部分与细胞器分类的数量最多；在生物学过程(Biological process)大类中，注释到代谢过程及细胞生理过程中的Unigene数量最多。13142条Unigene被注释到KEGG数据库的五大类中，其中注释到代谢(Metabolism)这一大类的Unigene数量最多，占总数的38%，其余注释到的大类中，有机体系(Organismal Systems)占20%，遗传信息处理(Genetic Information Processing)占17%，环境信息处理(Environmental Information Processing)占14%，细胞过程(Cellular Processes)占11%。在两组转录组数据分析中，共鉴定390个差异表达基因，其中有155个基因上调，235个基因下调。qRT-PCR结果显示所选取的9个基因的表达结果与转录组的结果基本趋于一致，结果说明了RNA-Seq数据的可靠性。细叶百合在碳酸盐(NaHCO3)逆境下的转录组数据对比，提供了许多的基因表达通路和表达量的差异，希望为细叶百合在盐碱逆境中的基因表达方式研究提供理论基础。