基于RNA-Seq技术的淹水胁迫下2个北美鹅掌楸无性系表达差异基因分析

郭雄昌，孙小艳，，伍祎翌，刘淑娟，李彦强，刘 晖

(1. 江西省科学院生物资源研究所，330096，南昌；

2. 江西省重金属污染生态修复工程技术研究中心，江西省科学院流域生态研究所，330096，南昌；
3. 江西省红壤及种质资源研究所，331717，南昌)

水淹胁迫是限制植物生长发育及影响植物分布的主要逆境环境因子之一。近年来，由于全球频繁出现如局部地区暴雨、洪涝灾害等极端气候，致使土壤长时间处于水分饱和状态，从而抑制绝大多数陆生植物生长，甚至导致死亡[1]。因此，全面解析植物的耐淹潜力，增强植物在淹水胁迫下的适应性，选育出耐淹植物品种，对淹水频发地区的生态恢复及环境改良具有重要意义。

植物对淹水胁迫的响应机制是一个复杂的过程，涉及多个信号和通路，参与的调控表达基因数量巨大。近年来，随着高通量测序技术的快速发展，植物逆境研究进入组学时代，其中，基于Illumina平台的转录组测序(RNA-Seq)技术以其高通量、低成本、组织检测差异性、系统性及重复性好等优势，被广泛应用于植物耐淹性研究中。目前，在水稻(Oryzasativa)[2-5]、玉米(Zeamays)[3,5-6]、拟南芥(ArabidopsisthalianaL.)[5]、黄瓜(Cucumissativus)[7]、落羽杉(Taxodiumdistichum)[8]、杨树(Populustomentosa)[9]、樱花(Cerasussachalinensis)[10]、猕猴桃[11]等多个物种中已利用高通量测序技术，从基因转录水平对植物淹水胁迫的分子响应机制进行了深入的研究，并成功分离、鉴定出一系列耐涝相关基因。所以，利用转录组测序技术分析植物耐水淹基因、揭示耐水淹分子机制、开展耐水淹品种选育具有重要意义。

北美鹅掌楸(Liriodendrontulipifera)树形挺直、花朵优美，是重要的园林绿化和用材树种，但在地下水位较高或排水不畅的地区栽植时，其根系易腐烂发黑，生长受到抑制。研究人员在鹅掌楸属植物耐淹机制和品种选育开展了大量工作，从根系形态发育、生理生化调控等不同方面阐述了鹅掌楸淹水适应机制[12-14]，并指出淹水胁迫下激活北美鹅掌楸抗氧化系统活性，诱发其形成通气组织等适应性机制均与抗逆响应的基因表达有关，所以，研究淹水条件下北美鹅掌楸基因表达特征对解释北美鹅掌楸耐淹机制具有重要意义。但截至目前，国内外尚未开展转录组测序技术在鹅掌楸属植物淹水胁迫响应方面的研究。因此，本研究以前期引种筛选出的1个耐淹北美鹅掌楸无性系和1个不耐淹北美鹅掌楸无性系为研究材料，利用RNAseq高通量测序技术，比较淹水胁迫4 d后不同北美鹅掌楸无性系的基因表达谱数据，通过生物信息学方法筛选淹水胁迫的差异表达基因，为进一步鉴定和淹水胁迫响应相关的关键基因，探讨北美鹅掌楸耐淹分子机制奠定理论基础。

1.1 研究材料

淹水实验在江西省科学院玻璃温室开展。基于前期研究，从北美鹅掌楸耐淹种源Louisiana筛选出1个耐淹无性系，不耐淹的Missouri种源中筛选出1个淹水敏感型无性系。取2个北美鹅掌楸无性系的扦插幼苗于春季2月未抽新叶时移植到规格25 cm×30 cm的盆栽中，培养基质为混合均匀的南方红壤、粗砂和蛭石(体积比3：2：1)。幼苗统一管理，待长出8～10片新叶后，选取高度50～60 cm、生长较一致的扦插幼苗开展淹水处理。

1.2 处理方式

淹水处理使用温室盆栽法：在45 cm×80 cm的塑料箱中放置盆栽容器(2盆/箱)，设置4个处理组，分别为：淹水敏感无性系对照(CS)，淹水敏感无性系淹水处理(TS)，耐淹无性系对照(CR)，耐淹无性系淹水处理(TR)。每个处理4个重复。

处理条件：淹水处理时保持水面高于盆栽容器3 cm左右，玻璃温室中日最高温32 ℃，夜最低温25 ℃，平均光照强度范围为1 000～2 000 μmol/s·m2，昼夜湿度60%±5%。

采样：淹水4 d时，取距顶部第6片叶作为构建转录组文库的样本，每个处理设置4个生物学重复。

1.3 高通量测序方法

采用Trizol试剂提取叶片总RNA后，用Nanodrop分光光度计和琼脂糖凝胶电泳分别检测总RNA浓度和质量；
mRNA分离与纯化依据Invitrogen公司的FastTrack®2.0 Kit说明书，再加入Fragmentation Buffer随机打断mRNA，随后添加随机引物和反转录酶合成cDNA第1链；
随后加入DNA聚合酶、RNase H和dNTPs合成第2链。纯化后的cDNA末端修饰后连接接头序列，在Illumina HiSeqTM PE100平台进行双末端测序(paired-end sequencing)，测序工作由北京诺禾致源科技股份有限公司完成。

1.4 转录数据分析

10 Gb测序数据量完成后过滤raw reads，去除接头序列及低质量reads，获得高质量测序序列，用OAP软件将clean reads从头组装成Unigenes。采用Gene Inde Assembly程序，将Unigenes与NCBI数据库中北美鹅掌楸EST序列进行合并，并用Corset程序对测序转录本进行层次聚类，以聚类后的序列作为参考基因组进行后续分析。使用BlastX将参考基因组与七大数据库(Nt，Nr，Pfam，KEGG，GO，KOG/COG，Swiss-prot)进行功能注释、功能分类、CDS(coding sequence)预测等。

以拼接得到的转录组作为参考基因组序列，比对每个样品中去除低质量后的clean reads；
利用RPKM(reads per kilobase per million mapped reads)值反映基因表达丰度，量化基因的表达水平。因为RPKM可消除基因长度和测序量差异导致的对基因表达计算的误差，所以以FDR﹤0.01和log2(RPKMT / RPKMC)>2.00为标准，可以直接鉴定淹水处理前后的差异表达基因。通过K-means聚类分析、SOM聚类分析等，将淹水上调表达基因和下调基因进行层次聚类，并通过GO(Gene Ontology)功能显著性富集分析，以不同家系淹水前后差异表达基因的分布情况，探究差异表达基因参与的生物学调控过程；
并通过KEGG Pathway富集分析，分析差异表达基因参与调控代谢途径、信号转导途径等。

2.1 高通量测序数据质量评估及其de novo组装

从北美鹅掌楸cDNA文库共获得210 235 018条raw data，去除带接头、低质量、N的比例大于10 % 的序列后，共获得201 817 794条clean reads；
碱基正确识别率达99 %的Q20设置的高质量序列占 96.16%，碱基正确识别率达99.9 %的Q30设置的高质量序列占90.75%，GC含量占总碱基48.20%，碱基错误率为0.02%，说明高通量测序平台获得的北美鹅掌楸转录数据质量较高，满足后续生物信息学分析的条件。

表1 北美鹅掌楸转录本拼接组装统计

采用Trinity软件对clean reads进行拼接获得参考转录本序列，获得395 733个转录本，214 161个Unigene，平均碱基长度984 nt，最长序列达到19 144 nt，最短序列为201 nt，而且随着unigenes长度增加，unigenes数量逐渐减少(表1)。淹水敏感北美鹅掌楸无性系在对照和淹水胁迫下分别有74.41%、73.26%的clean reads与参考转录本序列比对上，耐淹无性系分别有74.86%、73.93%的序列比对上参考转录本序列。

2.2 北美鹅掌楸转录组unigenes的功能注释

使用BLAST将总unigenes分别与Nr、Nt、Pfam、KOG/COG、Swiss-prot、KEGG、GO等7个数据库进行注释。结果显示，97 725条unigenes在Nr数据库比对成功，比对效率45.63%；
68 729条unigenes(32.09%)在Nt数据库比对成功，15 556条序列(7.26%)同时在所有数据库中注释。以NR数据库注释的序列为基础进行物种比对，结果显示与北美鹅掌楸大部分基因序列匹配最多的物种是荷花(Nelumbonucifera)和葡萄(Vitisvinifera)，分别占注释序列总数的30.08%和12.37%，与其他油棕(Elaeisguineensis)、海枣(Phoenixdactylifera)等物种匹配的同源序列均小于10%(图1)。

Nelumbo nucifera荷花；
Vitis vinifera葡萄；
Elaeis guineensis油棕；
Phoenix dactylifera海枣；
Amborella trichopoda无油樟；
Theobroma cacao可可；
Musa acuminata小果野蕉；
Hordeum vulgare大麦；
Gossypium raimondii二倍体棉花；
Jatropha curcas麻风树

2.3 差异表达基因的筛选

对比淹水胁迫处理后淹水耐性无性系TR和淹水敏感无性系TS与未胁迫前的对照组CR、CS的基因表达量后发现，当差异倍数设置为2时，结果显示在未进行淹水处理的条件下，淹水耐性无性系CR相对淹水敏感无性系CS的差异基因163个，其中80个基因上调表达，84个基因下调表达。淹水4 d后，TR相对CR有3 808个差异基因，其中2 015个上调表达基因，1 793个下调表达基因；
TS相对CS检测到10 389个差异基因，其中4 488个上调表达基因，5 901个下调表达基因；
TR相对TS的差异基因639个，其中325个基因上调表达，314个基因下调表达(图2)。差异基因统计结果表明，淹水胁迫下植物参与应激与防御涉及的调控途径和基因变化非常复杂。

图2 淹水胁迫前后差异表达基因数量统计

2.4 差异表达基因GO分类和KEGG通路分析

利用转录组比对获得的DGEs按照基因的生物反应、细胞组分、分子功能三大类中的具体功能进行GO分类。以校正P≤0.05为显著富集标准，3 808个DGEs至少含有一种GO功能注释(P<0.05)，其中注释到生物反应、细胞组分、分子功能中富集数量分别为2 415个、419个、891个。差异基因主要富集于生物反应中氧化还原过程、代谢过程、磷酸化、碳水化合物代谢过程、蛋白磷酸化，分子功能中的催化活性、激酶活性、蛋白激酶活性、转移酶活性，细胞组分中涉及的差异基因相对较少。1 256个淹水差异表达基因被KEGG数据库注释到113条通路(pathway)中。图3显示差异表达基因中显著性富集的前20条通路，结果显示光合作用、淀粉和蔗糖代谢、半乳糖代谢、果糖和甘露糖代谢、植物激素信号转导、氨基糖和核苷酸糖代谢、苯丙烷类生物合成等途径基因显著性富集。

图3 北美鹅掌楸耐淹无性系淹水后差异表达基因的KEGG分类

淹水敏感无性系淹水胁迫条件下差异基因高达10 389个差异基因，下调表达基因的数量显著大于上调表达的基因数目。GO注释结果显示差异表达基因在GO数据库中生物反应、细胞组分、分子功能中富集数量分别为5 873个、1 533个、3 088个。差异基因主要体现在分子功能和生物学过程，其中生物反应中的差异基因主要富集于氧化还原过程、代谢过程、单一生物代谢运输、跨膜转运、单一生物过程、单一生物分解代谢过程等反应中，分子功能分类中的氧化还原酶活性、催化活性、铁离子结合、血红素结合、四吡咯结合、辅助因子结合、黄素腺嘌呤二核苷酸结合、辅酶结合等功能相关的差异基因显著富集。在淹水敏感无性系，3 156个差异基因被KEGG数据库注释到123条通路，基因显著富集于苯丙烷类生物合成、植物激素信号转导、脂肪酸降解、精氨酸和脯氨酸代谢、光合作用、苯丙氨酸代谢、半乳糖代谢等途径(图4)。

图4 北美鹅掌楸淹水敏感无性系淹水后差异表达基因的KEGG分类

淹水敏感无性系和耐淹无性系淹水后的163个显著差异表达基因主要富集与糖代谢、苯丙氨酸、酪氨酸和色氨酸生物合成，甘氨酸、丝氨酸、苏氨酸、丙氨酸、天冬氨酸和谷氨酸代谢，以及植物激素转导等途径，说明淹水导致的厌氧胁迫显著影响了北美鹅掌楸的基础代谢，氨基酸代谢途径与鹅掌楸耐淹能力具有密切相关性。

植物耐涝性是一个复杂的综合性状，淹水4 d后，北美鹅掌楸耐淹无性系的生长状态明显比淹水敏感无性系表现良好，植株叶片萎蔫黄化情况较少，顶梢出现短暂萎蔫后恢复一定活力，而淹水敏感无性系根系出现明显腐烂黑化的现象，叶片黄化脱落，顶梢萎蔫失去活力，从淹水胁迫表型响应说明不同无性系淹水耐性能力存在较大差异[13]。

转录组是功能基因组研究的重要技术，受淹水胁迫后不同阶段和不同组织部位的植物基因表达水平及表型差异较大，具有特定的时间和空间特异性。本研究首次利用Illumina高通量测序技术对北美鹅掌楸耐淹无性系和淹水敏感无性系淹水前后的叶片转录组进行测序，获得395 733个转录本，214 161个Unigene，为进一步挖掘鹅掌楸属耐淹相关的功能基因提供了丰富的数据信息。

利用基因差异表达量阈值分析，发现不同无性系在淹水胁迫后，表达差异基因类型存在较大差异。在耐淹无性系，淹水后上调的基因数量大于下调的基因数量；
淹水敏感无性系中，下调的基因数量显著大于上调的基因数量。Arismendi等[15]发现樱桃缺氧胁迫6 h后只有764 DEGs差异表达，其中65%的DEGs上调表达，淹水胁迫72 h时差异表达基因数量最多。而玄参科植物地黄淹水处理后筛选出7 249个DEGs，2 505个(约35%)基因上调表达，4 744个(约65%)基因下调表达[16]。Thirnavukkarasu等[4]对玉米淹水胁迫转录组研究发现，淹水胁迫后，HKI 1105(耐淹型)中1 188个基因上调表达，891个下调表达，而V372(淹水敏感型)1 337个基因上调表达，1 040个下调表达。Qi等[7]利用高通量标签测序发现黄瓜淹水胁迫2 h检测到1 180个DGEs，其中551个上调表达，629个下调表达，且随着淹水时间延长至4 h、8 h、24 h时下调基因比例持续增加，分别为900个(55.3%)、2 116个(67%)、2 153个(66%)，说明淹水胁迫抑制基因表达，而上调表达的基因主要参与碳循环、光合作用、激素介导的信号转导，以及活性氧产生和清除等通路。所以，植物逆境胁迫下的转录基因差异变化随胁迫时间和品种存在较大差异，而且敏感材料在淹水胁迫后表现出更强烈的差异基因表达水平[17]。

光合作用变化是植物受淹水、干旱等逆境胁迫的最直接反映，叶片光合作用途径相关的主要基因在淹水胁迫下普遍下调表达[10,18]。位于类囊体膜上的光合作用光反应包括光合作用色素(2个反应中心和2个光系统，PSII和PSI)、光合作用电子传递链、ATP合成元件等。在环境变化过程中，植物通过动态地调节叶绿体色素的分配和耗散光能保障其光合系统的正常功能，耐性较强的植物品种能够维持较高的光合速率[19]。在淹水北美鹅掌楸叶片中，淹水敏感无性系检测到编码光反应(Light Harvesting Chlorophy Protein Complex，LHC)相关蛋白的35个差异基因均显著下调表达，而耐淹无性系检测到8个编码LHC I相关蛋白的差异基因显著下调表达，编码LHC Ⅱ相关蛋白的差异基因15个下调表达，而4个Lhcb1基因显著上调表达，这表明淹水后北美鹅掌楸光合速率普遍下降，但耐淹无性系可能通过调节叶绿体色素的分配和效率保持相对较高的光合速率以维持其生存。

淹水导致的低氧或缺氧胁迫抑制植物的有氧呼吸过程，ATP生成受到严重抑制，从而导致一系列依赖ATP的生理代谢过程受到抑制，植物迅速激活糖酵解、乙醇发酵、乳酸发酵以及丙氨酸发酵途径等以维持植物根系ATP的产生，所以在淹水胁迫下植物能量代谢途径受到巨大影响[8,20-21]。本研究中，在淹水胁迫下，耐淹水和淹水敏感无性系之间存在大量共差异表达基因，这些基因显示在淹水胁迫下北美鹅掌楸启动了很多共同的基因应答途径。大量研究表明，淹水胁迫下植物的差异表达基因主要集中在能量代谢、光合作用、碳水化合物合成、糖酵解、氨基酸代谢、植物激素信号转导等代谢通路，以保证植物在水淹胁迫下维持生长所需的能量和物质供应[10,16]。富集基因的KEGG结果表明，在淹水胁迫下，耐淹无性系叶片中与糖代谢、激酶活性相关的代谢途径，如淀粉和蔗糖代谢、半乳糖代谢、果糖和甘露糖代谢、氨基糖和核苷酸糖代谢等差异基因的数量较多、基因富集程度较显著。缺氧环境下，碳水化合物作为植物首要能量来源，植物的耐淹型主要取决于可溶性糖、淀粉等碳水化合物的利用效率和储备[22-23]。与蔗糖合成的反向途径相比，蔗糖合成途径消耗能量较少，蔗糖合成反向酶降解途径中，每降解1分子蔗糖消耗2分子ATP，而蔗糖合成酶途径每合成1分子蔗糖消耗1分子ATP和1分子UDP[24]。所以，北美鹅掌楸耐淹无性系通过激活淀粉、蔗糖等碳水化合物合成途径的相关基因，激活糖酵解、乙醇发酵等途径，来维持其在缺氧环境下生存所需的能量消耗[25]。

总之，本研究通过高通量测序技术，获得大量北美鹅掌楸淹水诱导基因，而且敏感无性系在淹水胁迫后表现出更丰富的差异基因表达水平，通过Pathway显著性富集分析发现淹水胁迫后光合作用、植物激素信号转导等途径基因显著性富集，但耐淹无性系中淀粉和蔗糖代谢、果糖和甘露糖代谢，以及糖酵解等途径基因显著性富集。因此，利用RNA-seq技术获得差异表达基因为进一步深入研究北美鹅掌楸耐淹胁迫调控网络和分子调控机制提供了理论依据，也为林木的耐淹研究提供了基因资源。

猜你喜欢差异基因高通量北美北美灰熊被杀案小哥白尼(野生动物)(2022年3期)2022-06-16高通量血液透析临床研究进展透析与人工器官(2020年1期)2020-11-16基于RNA 测序研究人参二醇对大鼠心血管内皮细胞基因表达的影响 （正文见第26 页）心电与循环(2020年1期)2020-02-27Ka频段高通量卫星在铁路通信中的应用探讨铁道通信信号(2019年8期)2019-10-10北美纪行中国三峡(2017年9期)2017-12-19中国通信卫星开启高通量时代中国发展观察(2017年8期)2017-04-26紫檀芪处理对酿酒酵母基因组表达变化的影响江苏农业科学(2017年5期)2017-04-15北美星鸦知道松子藏在哪儿学苑创造·A版(2016年8期)2016-07-06护理干预在高通量血液透析患者中的应用效果中国当代医药(2015年33期)2015-03-01只身闯北美海峡姐妹(2015年7期)2015-02-27