甘蓝型油菜盐胁迫响应基因表达谱分析及共表达网络的构建

杨 闯 王 玲 全成滔 余良倩 戴 成 郭 亮 傅廷栋 马朝芝

华中农业大学作物遗传改良全国重点实验室 / 国家油菜工程技术研究中心 / 洪山实验室, 湖北武汉 430070

土壤盐渍化是威胁全球农业生产的主要环境挑战之一[1]。世界上大约有20%的灌溉农田会受到土壤盐渍化的不利影响[2]。自然环境恶化、灌溉方式不当以及气候变化等因素也正在进一步加剧土壤盐渍化问题[3]。预计到2050年, 盐渍化土壤面积将达到全球可利用耕地的50%以上[4]。近些年来, 我国的农业耕地也在不断受到盐渍化侵蚀, 盐渍地总面积已达到3.69×107hm2, 给作物生产和农业可持续发展带来了极大的挑战[5-6]。油菜是世界四大油料作物之一, 具有重要的油用、饲用、菜用以及旅游观光价值[7]。作为食用植物油和植物蛋白的主要来源, 油菜在我国农产品中占有重要地位, 对于保障我国油料供给安全尤为重要[8]。然而, 盐胁迫目前是影响油菜生长发育的主要环境条件之一, 可能会造成油菜减产、品质下降甚至死亡。因此, 解析甘蓝型油菜响应盐胁迫的机制, 并从中挖掘油菜响应盐胁迫的核心基因, 对于油菜的耐盐性改良具有重要意义。

盐胁迫对植物的毒害作用主要包括渗透胁迫和离子胁迫2种[9]。盐胁迫的渗透作用会导致土壤溶液的水势低于根细胞水势, 从而抑制根系的吸水作用; 此外, 渗透胁迫还会导致叶片气孔关闭从而使光合效率下降, 最终影响植物的正常生长发育。离子胁迫主要是由细胞中的Na+和Cl-过量积累导致的。Na+的毒性主要是通过抑制新陈代谢途径中关键酶活性从而影响植物的生长发育; 同时,过量的Na+也会影响植物对氮、磷、钾、钙等其他矿质元素的吸收利用。由于植物对NO3-、SO42-以及Cl-的吸收都是通过相同的非选择性阴离子通道, 因此过量的Cl-可能会导致氮素以及硫素的缺乏。

植物为了响应盐胁迫, 也进化出了互相关联的调节通路并引起细胞内的广泛变化, 其中许多变化是适应性反应, 可导致抗逆性增加, 这也是植物在长期进化过程中保留下来的适应性机制[10]。盐胁迫早期响应阶段, 植物可以通过高渗感受器OSCA1以及低渗感受器MSL8感受渗透压的变化[11-12], 同时植物还可以通过钠离子感受器MOCA1感受盐浓度变化[13]。植物感受到盐胁迫后, 可通过渗透调节、组织特异性离子转运、激素信号传导、以及活性氧清除等机制来应对盐胁迫, 其中也涉及到一系列基因的表达, 包括相关功能基因和调节基因, 如转录因子、蛋白激酶等。SOS信号通路是近20年来发现的典型的盐胁迫响应通路,SOS1、SOS2以及SOS3编码构成SOS信号通路的蛋白质, 在多种植物的盐胁迫响应过程中发挥着重要功能[14]。AtNHX1和AtNHX2编码液泡膜上的Na+/H+逆向转运蛋白, 主要在根和叶中表达, 并选择性地将Na+运输到液泡, 其突变体对盐胁迫更加敏感[15]。其中,转录因子作为关键的调控因子, 其介导的基因调控网络在植物响应盐胁迫时起着不可或缺的作用。拟南芥AtNIG1是植物中首先被鉴定出来参与耐盐性的BHLH转录因子家族成员, 过表达AtNIG1表现出更高的耐盐性[16]。玉米转录因子ZmWRKY114作为负调控因子, 通过ABA依赖性途径参与盐胁迫响应[17]。水稻转录因子基因OsDOF15通过介导乙烯生物合成途径, 在盐胁迫抑制初生根伸长的过程中起着关键作用[18]。棉花转录因子基因GhABF3可以通过降低细胞氧化程度显著提高陆地棉的耐盐性[19]。大豆转录因子基因GmNAC06可以通过控制毛状根Na+/K+比值维持离子稳态, 从而增强大豆的耐盐性[20]。这些研究结果表明, 相关转录因子的鉴定及功能解析对于理解植物响应盐胁迫的分子机制至关重要。

随着高通量测序技术的快速发展, RNA-seq技术已经被广泛地用于多种植物的研究中, 同时也为转录组学研究提供了大量的分析方法。权重基因共表达网络分析(WGCNA)是最常用的分析RNA-seq数据的方法之一。相比于其他分析方法, WGCNA在产生具有生物学意义的结果方面也更加稳健[21]。该方法可以通过构建基因共表达网络将上万个基因划分成十几个具有相似表达模式的基因模块, 并通过分析与性状高度相关的关键模块从中挖掘核心基因。Zheng等[22]利用324个RNA-Seq数据构建了玉米表皮蜡质合成相关基因的共表达网络, 发现几乎所有已知的有关蜡质合成的基因都在同一个模块, 这为后续深入挖掘玉米表皮蜡质合成相关基因提供了重要的参考依据。李旭凯等[23]利用WGCNA预测出与多种胁迫高度相关的25个抗逆关键基因, 为水稻的多抗逆基因挖掘以及利用提供了参考依据。Ye等[24]利用盐梯度处理的RNA-seq数据, 阐明了互花米草高耐盐性的转录组动力学, 并通过WGCNA分析发现SaOST1等蛋白激酶编码基因是耐盐调控网络中的核心基因。Zhao等[25]以花生为材料, 利用WGCNA筛选出一个与干旱胁迫生理数据密切相关的重要模块, 并从中预测到9个抗旱候选基因。Yang等[26]基于剪接异构体的WGCNA分析, 筛选出与不同非生物胁迫相关的23个共表达模块, 并发现许多重要的核心基因在单一应激或多种应激反应中起关键作用, 其中主要包括转录因子、RNA结合蛋白等。

利用WGCNA构建基因共表达网络, 已经被广泛地应用于多种植物的转录组研究以及核心基因筛选。本研究通过测得的90份RNA-seq数据, 获得了甘蓝型油菜盐胁迫处理前后的高分辨率时间动态转录表达谱, 进一步通过差异基因分析和WGCNA分别构建油菜根系以及叶片组织响应盐胁迫的基因共表达网络, 并从中挖掘油菜响应盐胁迫的核心基因。研究结果可为甘蓝型油菜的盐胁迫共表达网络研究以及油菜的耐盐性分子育种提供理论支持。

1 材料与方法

1.1 植物材料和盐胁迫处理

选取甘蓝型油菜ZS11作为试验材料, 采用水培法进行培养, 培养条件为: 光照16 h、黑暗8 h, 温度24℃左右, 相对湿度为75%。培养27 d后, 选取长势一致的幼苗,使用200 mmol L–1氯化钠模拟盐胁迫处理, 分别在盐胁迫处理后0、0.25、0.5、1、3、6、12、24 h进行取样, 同时各个时间点均用正常培养条件下的幼苗作为对照组, 取样时分离根部组织和叶片组织, 立即将其包入锡箔纸并投入–80℃的液氮中储存, 每个处理3个生物学重复。取样完毕后送公司进行RNA提取以及转录组测序。盐胁迫处理条件如图1所示, 其中0 h的样本组未显示。本研究所使用的原始数据尚未公开发表, 但单个基因受到盐胁迫处理后的转录表达量已上传并可通过BnIR[27]网站(http://yanglab.hzau.edu.cn/BnIR/eFP_single_gene)查询。

图1 盐胁迫样品处理条件Fig. 1 Treatment conditions of salt stress samples

1.2 RNA-seq数据处理与分析

1.2.1 RNA-seq数据的质控与定量 使用FastQC和MultiQC软件[28]对RNA-seq数据进行质控, 使用Trim Galore过滤掉低质量的reads片段, 获得clean reads用于后续分析。使用Hisat2[29]将其批量比对到ZS11参考基因组, 统计所有样本的比对率。使用Samtools[30]将sam文件转换成bam文件, 使用FeatureCounts[31]对回帖到参考基因上的CDS序列进行计数, 获得原始counts值。根据counts值以及CDS序列长度, 使用TBtools[32]计算出基因的TPM值。同时使用原始counts值, 利用R包DESeq2[33]进行标准化, 对盐胁迫根系以及叶片组织在各个时间点的处理组样本与其对照组样本进行两两差异基因分析,最终获得组织在各个处理时间点的差异基因以及上下调基因数目。筛选差异基因的阈值设置为|log2(Fold Change)|>1 & FDR<0.05。

1.2.2 利用WGCNA构建基因共表达网络 WGCNA(weighted gene co-expression network analysis), 即权重基因共表达网络分析, 可以通过构建共表达网络寻找协同表达的基因模块, 探究基因模块与样本性状之间的相关性, 从中挖掘核心基因, 经常适用于非生物逆境胁迫不同时间点应答、病原菌侵染后不同时间点应答等研究方向[34]。本研究考虑到低表达基因一般不具备生物学意义, 将Deseq2分析得到的盐胁迫组织响应差异基因赋予TPM值,进一步筛选得到至少在3个样本中TPM值大于1的组织响应差异表达基因用于WGCNA的输入值。最小模块的基因数不低于100个, 其余参数按照默认设置。使用R包Cluster Profiler[35]对模块进行GO富集分析, 阈值设置为FDR<0.05。本研究定义模块MEs与样本时间点之间的相关性绝对值大于0.6,P值小于0.05为显著相关; 以权重值(weight)大于0.2作为基因间的有效连接, 连通度前10%的基因为核心基因。使用Cytoscape软件对关键模块的共表达网络进行可视化分析[36]。

2 结果与分析

2.1 盐胁迫转录组数据的产生及初步验证

本研究于盐胁迫处理0、0.25、0.5、1、3、6、12、24 h对油菜叶片组织和根系组织进行转录组测序, 获得90份RNA-seq数据。通过质检和过滤, 共获得52.78亿条clean reads。将其比对到ZS11参考基因组, 发现所有样本的比对率都超过95%, 表明本研究所获得的RNA-seq数据具备较高质量。为了验证盐胁迫处理的有效性, 利用拟南芥已报道参与盐胁迫响应的marker基因, 观察其对应的油菜同源基因在盐胁迫处理前后的表达模式。SOS1、WRKY33通常在盐胁迫早期响应阶段被激活[37-38], 油菜根系组织中的BnaC09.SOS1以及BnaC04.WRKY33在盐胁迫1 h内表达量显著上升(图2-A, B)。HKT1作为Na+运载体, 可以装载叶片中的钠离子到韧皮部再循环到根部[39],油菜叶片组织中的BnaA09.HKT1在盐胁迫处理3 h时出现显著差异表达(图2-C)。ABI1可以负调控ABA促进气孔关闭从而影响植物耐盐性[40], 油菜叶片以及根系组织中的BnaA03.ABI1在盐胁迫处理3 h时表达量显著上升(图2-D)。以上结果均与拟南芥中的报道结果类似, 这表明本研究对甘蓝型油菜幼苗进行的盐胁迫处理是有效的,也进一步验证了RNA-seq数据的可靠性。

图2 盐胁迫响应marker基因的表达模式Fig. 2 Relative expression pattern of marker genes under salt stress

2.2 盐胁迫的相关性分析

为了探究样本处理之间的相关性, 本研究进行了主坐标分析以及层次聚类分析。从盐胁迫的主坐标分析结果来看, 第1主坐标主要解释的是组织表达的差异性, 其贡献度较大, 约占75.08%, 表明后续在利用WGCNA构建基因表达网络时, 应将胁迫的根系和叶片组织分开分析,避免组织差异性覆盖胁迫处理差异性; 第2主坐标主要解释的是对照组以及盐胁迫处理之间的差异性, 盐胁迫处理3 h、6 h、12 h、24 h的样本表达与对照组形成明显差异,这也表明胁迫处理的有效性和数据的良好性(图3-A)。从层次聚类结果来看, 根和叶由于组织表达的差异性较大聚为2类, 这与主坐标分析结果相对一致。另外, 无论是根系组织还是叶片组织, 1 h处理前后的组织样本均形成了明显的聚类差异。根据组织响应盐胁迫的样本聚类结果进行初步划分, 可以将盐胁迫处理1 h前后划分为早期响应阶段以及后期响应阶段(图3-B)。

图3 盐胁迫的相关性分析Fig. 3 Correlation analysis of salt stress

2.3 差异基因分析

本研究利用DESeq2包进行差异基因分析, 获得了油菜根系以及叶片组织在盐胁迫各个处理时间点的差异基因以及上下调差异基因数目, 并分析了差异基因数目随盐胁迫处理时间的变化趋势。盐胁迫处理1 h内, 叶片组织在0.25 h、0.5 h和1 h的差异基因数目始终小于根系组织, 表明根系组织在早期阶段的响应程度大于叶片组织; 盐胁迫处理1 h后, 叶片组织在任意一个处理时间点的差异基因数目都大于根系组织, 表明叶片组织在后期阶段的响应程度大于根系组织(图4-A)。其中, 根系组织在盐胁迫处理3 h时达到最高的响应程度, 此时差异基因为10,883个; 而叶片组织在处理24 h时达到最高的响应程度, 此时差异基因数目为18,414个。进一步分析盐胁迫上下调差异基因数目的变化趋势, 发现盐胁迫处理24 h内, 根系组织中的上调基因数目始终大于下调基因, 而叶片组织仅在处理3 h内的上调基因数目大于下调基因, 而随着盐胁迫时间的增长,叶片组织中的上调基因数开始降低并逐渐少于下调基因(图4-B)。表明油菜幼苗对盐胁迫的响应程度存在组织差异性。盐胁迫处理12 h内, 根系组织和叶片组织的响应程度均呈现出先升高后下降的变化趋势, 变化拐点均出现在处理3 h附近, 表明油菜幼苗对盐胁迫的响应程度也存在组织相似性。本研究定义至少在任意一个时间点存在显著差异的基因作为盐胁迫组织响应差异基因集, 通过去重取并集, 得到根系组织以及叶片组织响应差异基因分别是20,462个和29,334个, 表明全程盐胁迫处理24 h内, 叶片组织整体上的响应程度比根系更剧烈。

图4 盐胁迫处理24 h内的差异基因数目及变化趋势Fig. 4 Number and trend of differential genes under salt stress treatment for 24 hours

2.4 权重基因共表达网络的构建

为了构建油菜幼苗响应盐胁迫的基因共表达网络,并从中挖掘核心基因, 基于盐胁迫处理24 h的转录组表达矩阵用于WGCNA分析。同时, 为了防止组织差异性过大影响分析结果, 在保证不低于WGCNA所要求的样本数的情况下, 分别构建了根系组织以及叶片组织响应盐胁迫的基因共表达网络。过滤低表达基因, 筛选26,839个叶片组织差异表达基因, 使用动态树切割方法, 选取最佳软阈值β=14构建叶片组织响应盐胁迫的基因共表达网络(图5-A), 共识别出14个不同的共表达模块, 最大和最小模块分别是turquoise和tan模块, 分别含有8651个和260个基因(附表1)。过滤低表达基因, 筛选18,952个根系组织响应差异表达基因, 选取最佳软阈值β=18构建根系组织响应盐胁迫的基因共表达网络(图5-B), 共识别出11个不同的共表达模块, 最大和最小模块分别是turquoise以及greenyellow模块, 分别含有4058个和234个基因(附表2)。

图5 最佳软阈值的选择Fig. 5 Selection of optimal soft threshold (power)

2.5 早期及后期响应模块的鉴定

在叶片组织响应盐胁迫的共表达网络中, 通过计算模块与样本处理时间点之间的相关性发现, tan模块、salmon模块分别与盐胁迫早期响应阶段0.5 h、1 h显著正相关; red模块、cyan模块、turquiose模块以及green模块分别与盐胁迫后期响应阶段3 h、6 h、12 h、24 h显著正相关(图6-A)。在根系组织响应盐胁迫的共表达网络中,通过计算模块与样本处理时间点之间的相关性发现, yellow模块与盐胁迫早期响应阶段0.25 h以及0.5 h显著正相关; blue模块和turquiose模块分别与盐胁迫处理1 h、3 h显著正相关; black模块与盐胁迫处理12 h显著正相关;purple模块以及red模块均与盐胁迫处理24 h显著正相关(图6-B)。综合考虑模块与样本处理时间点的相关性、模块内的转录因子所占比例以及模块的GO富集结果, 选取叶片组织中的tan模块和根系组织中的yellow模块作为组织早期响应盐胁迫的关键模块, 叶片组织中的green模块和根系组织的red模块作为组织后期响应盐胁迫的关键模块。叶片组织中的早期响应tan模块包含260个基因, 其中有75个转录因子, 该模块内的绝大多数基因在盐胁迫处理0.5 h的表达量显著高于其他时间点(图7-A)。根系组织中的早期响应yellow模块包含2322个基因, 其中含有276个转录因子, 该模块内的绝大多数基因在盐胁迫处理0.25 h以及0.5 h显著表达(图7-B)。叶片中的后期响应green模块包含1654个基因, 其中含有138个转录因子;根系中的后期响应red模块包含1238个基因, 含有99个转录因子, 2个模块内的绝大多数基因均在盐胁迫处理24 h时的表达量显著高于其他时间点(图7-C, D)。

图6 基因模块与样本时间点的相关性热图Fig. 6 Correlation heat map between gene modules and sample time point

图7 早期及后期响应模块的基因表达谱Fig. 7 Gene expression profiles of early and late response modules

2.6 早期及后期响应模块的GO富集分析

为了分析组织响应盐胁迫的共表达模块特征, 本研究分别对早期响应阶段的tan模块和yellow模块, 后期响应阶段的green和red模块进行GO富集分析。早期响应阶段, 根系最先受到盐胁迫的刺激, 根系组织中的yellow模块富集到了对生物刺激反应的调控、对外部刺激反应的调控、磷酸化信号转导、对盐的响应等生物学过程(图8-B)。叶片中的tan模块富集到与叶片衰老、盐胁迫响应的正向调控、自噬调节等生物学过程(图8-A)。通过对比分析发现, 早期响应模块均富集到了脱落酸、乙烯、茉莉酸等激素信号通路, 这可能表明油菜幼苗早期响应盐胁迫就会有多种激素参与全身器官的系统性反应。早期响应阶段,叶片中还富集到了对细菌的防御反应和对昆虫的防御反应, 根系也富集到了对细菌防御反应的调控以及对细菌来源分子的响应, 这可能表明油菜幼苗早期响应盐胁迫可能与早期响应细菌等生物胁迫有相似之处。除此之外,叶片中也富集到了抗毒素合成通路, 根系中富集到了类黄酮生物合成过程的调控。这些早期响应阶段的生物学通路富集, 可能表明油菜幼苗在早期响应盐胁迫时也会通过合成抗毒素等物质来响应盐胁迫。盐胁迫后期响应阶段,green和red模块显著富集到了对过氧化氢的响应、细胞对盐胁迫的反应、对脱水反应的调控等盐胁迫响应相关通路(图8-C, D)。其中, 叶片中的green模块也富集到了与早期响应相似的生物学通路, 如脱落酸激活信号通路;green模块也特异富集到了有机羟基化合物合成、肌醇生物合成等通路(图8-C)。根系中的red模块特异富集到了有机酸转运、植物型次生细胞壁的发生、苯丙烷生物合成过程、硫代葡萄糖苷代谢过程、木质素代谢过程等通路(图8-D)。总体来看, 无论是早期响应模块, 还是后期响应模块, 都富集到了盐胁迫响应的相关通路, 表明利用WGCNA可以构建具有生物学意义的共表达模块, 这些模块可成为本研究的重点关注对象。

(图8)

图8 早期及后期响应模块的GO富集分析Fig. 8 GO enrichment of early and late response modules

2.7 核心基因的鉴定与分析

为了识别关键模块中的核心基因,利用cytoscape进行共表达网络可视化。以权重值大于0.2作为阈值, 筛选连通度最大的前10%基因作为核心基因, 早期响应tan模块筛选到23个核心基因, 包含15个核心转录因子(图9-A)。其中连通度最高的3个核心基因分别为BnaA04G0009600ZS(AGP20)、BnaC04G0304800ZS(CAMBP25)、BnaA10G 0150300ZS(RING1)。BnaA04G0009600ZS(AGP20)编码调节细胞壁扩张的阿拉伯半乳糖蛋白, 其表达主要由盐胁迫诱导[41]。早期响应yellow模块筛选到221个核心基因,包含26个核心转录因子 (图9-B)。其中连通度最高的3个核心基因分别为BnaC09G0584800ZS(NHL3)、BnaC08G 0479500ZS(WAKL8)、BnaA09G0091000ZS(SUPA)。BnaA09G 0091000ZS对应的拟南芥基因SUPA定位于过氧化物酶体中, 盐胁迫处理15 min内迅速上调表达, 过表达可导致拟南芥中活性氧类物质积累增加以及抗逆性改变[42]。后期响应green模块筛选到核心基因137个, 包含14个核心转录因子(图9-C)。其中连通度最高的3个核心基因分别为BnaA03G0371900ZS(AOX1A)、BnaC07G0256700ZS(PBL19)、BnaA02G0027300ZS(DGK1)。BnaA03G 0371900ZS对应的拟南芥基因AOX1A编码线粒体电子传递链的替代氧化酶, 会限制脯氨酸诱导的氧化应激, 从而利于拟南芥盐度恢复[43]。后期响应red模块筛选到核心基因115个, 包含12个核心转录因(图9-D)。其中连通度最高的3个核心基因分别为BnaC05G0203500ZS(PME6)、BnaA02G0180600ZS(WRKY57)、BnaC06G0284600ZS(KTI1)。BnaC05G0203500ZS对应的拟南芥PME6编码果胶甲酯酶,在保卫细胞中高度表达, 并且是气孔功能所必需的, 基因恢复可以挽救保卫细胞壁果胶甲酯化状态、气孔功能和植物生长[44], 该基因也在油菜盐胁迫后期显著差异表达,可能通过气孔调节应对盐胁迫。

图9 早期及后期响应模块的基因共表达网络可视化Fig. 9 Co-expression network visualization of early and late response modules

转录因子作为关键调控因子, 其介导的基因调控网络经常在植物响应盐胁迫时发挥着重要作用。本研究着重关注核心转录因子, 通过双向blast将核心基因比对至拟南芥, 并借助TAIR数据库注释核心转录因子的功能(附表3), 发现油菜同源基因WRKY33、DDF1、AZF2、ARR1、ZFHD1、DREB2B等是拟南芥中已知参与编码盐胁迫响应的核心转录因子。其中, 拟南芥中的同源基因WRKY33对应油菜中的多个拷贝BnaA04G0247200ZS、BnaC03G 0217300ZS、BnaA03G0185200ZS、BnaC04G0562300ZS, 4个拷贝均在早期响应tan模块中占据较高连通度, 并参与多种非生物胁迫响应, 尤其是对盐胁迫的响应[37]。Yellow模块中的早期响应基因BnaA06G0081500ZS、BnaC05G 0100700ZS在拟南芥中的同源基因DDF1编码ERF/AP2转录因子家族DREB亚家族成员, 过表达可增强对高水平盐的耐受性[45]。Yellow模块中的早期响应基因BnaC05G 0398500ZS在拟南芥中的同源基因AZF2编码锌指蛋白转录因子, 响应ABA、高盐和轻度脱水反应[46]。Green模块中的后期响应基因BnaA01G0191600ZS在拟南芥中的同源基因ARR1属于细胞分裂素信号成分拟南芥响应调节因子1, ARR1蛋白可以通过MPK3/6负调控作用保持稳定从而增强耐盐性[47]。Red模块中的后期响应基因BnaA07G 0311300ZS在拟南芥中的同源基因ZFHD1编码锌指转录因子家族成员, 可由干旱、高盐和脱落酸诱导表达[48]。无论是早期响应模块还是后期响应模块, 均有已知的核心转录因子参与拟南芥的盐胁迫响应, 这表明本研究所筛选的关键模块以及核心转录因子具有较高的可靠性。据此初步推测, 这些核心基因可能也在甘蓝型油菜的耐盐性过程中发挥着重要作用。

2.8 核心转录因子在505份盐胁迫油菜群体变异数据中的SNPs及单倍型分析

进一步利用已发表的505份盐胁迫处理油菜群体变异数据库, 分析核心转录因子在甘蓝型油菜群体变异中的SNPs及单倍型情况(http://yanglab.hzau.edu.cn/BnIR/single_locus)[49]。结果表明, 多个核心转录因子在505份油菜群体中存在较为丰富的SNPs变异以及单倍型类型。其中,BnaC04G0015700ZS(BnWRKY46)是根系组织早期响应核心基因, 其在盐胁迫处理0.5 h时可上调表达约6倍。BnaA07G0271900ZS(BnWRKY57)是根系组织后期响应核心基因, 其在盐胁迫后期处理24 h可达到最大上调倍数约7倍。BnWRKY46和BnEWRKY57均由2个内含子和3个外显子组成, 经在群体内序列变异分析发现, 这2个核心基因均呈现出较为丰富的SNPs变异(图10-A, C)。单倍型分析发现,BnWRKY46可初步划分为3种单倍型, 不同单倍型在盐胁迫处理后的根长变短, 其中hap.B和hap.C型在盐胁迫处理前后的根长相对值显著小于hap.A型(图10-B)。BnWRKY57可划分为3种单倍型, 不同单倍型的根长在盐胁迫处理后显著缩短, 其中hap.C型在盐胁迫处理前后的根长相对值显著大于hap.A和hap.B型(图10-D)。表明本研究所挖掘到的核心转录因子中很有可能存在参与盐胁迫响应的关键候选基因, 深入鉴定盐胁迫相关的极端单倍型也将有利于为耐盐型油菜提供优质的种质资源。

图10 核心基因BnWRKY46和BnWRKY57在505份群体数据中的SNPs及单倍型分析Fig. 10 SNPs and haplotype analysis of hub genes BnWRKY46 and BnWRKY57 in 505 population data

3 讨论

3.1 甘蓝型油菜响应盐胁迫的时间动态转录表达谱分析

油菜作为世界上重要的油料作物, 在生长过程中经常受到盐胁迫的影响, 可能会导致减产、品质下降甚至死亡。而甘蓝型油菜被认为是中等耐盐作物, 因此解析甘蓝型油菜响应盐胁迫的调控网络, 对于油菜的耐盐性改良以及盐碱地的开发和利用都具有重要意义。本研究利用转录组测序获得了甘蓝型油菜幼苗在盐胁迫处理前后的90份RNA-seq数据, 将其比对到ZS11参考基因组, 发现所有样本的比对率都超过95%, 盐胁迫响应marker基因的表达模式也与拟南芥报道结果类似, 表明盐胁迫处理是有效的, 数据质量是可靠的。

由于差异基因数目可以在一定程度上反映胁迫响应程度, 本研究发现油菜幼苗对盐胁迫的响应程度存在组织差异性和相似性。全程盐胁迫处理24 h内, 叶片组织整体上的响应程度比根系更剧烈。根系组织24 h内响应盐胁迫的差异基因动态变化趋势与2007年发表的拟南芥盐胁迫24 h内的转录动态极为类似, 均呈现出先上升后下降的趋势。并且, 根系组织在24 h内的任意处理时间点的上调基因数目始终大于下调基因数目, 这一点也与拟南芥中的研究结果相对一致[50]。然而, 不同之处在于, 油菜中的组织响应差异基因数目远高于拟南芥, 这可能是由于异源四倍体的油菜经历过自然加倍等事件具有更多的同源基因。综上所述, 本研究获得了甘蓝型油菜响应盐胁迫的高分辨率时间动态转录表达谱, 可为后续研究人员提供可靠的数据参考。

研究人员曾在拟南芥中发现非生物逆境胁迫早期反应与后期反应的转录差异, 并认为早期反应可能是非特异性的, 主要是在胁迫处理1 h内, 而基因转录的特异性应激反应是在1～3 h后被检测出来[51]。根据盐胁迫样本层次聚类结果, 本研究也发现油菜响应盐胁迫存在明显的早期响应和后期响应的聚类差异。对于全程盐胁迫处理24 h的样本, 处理1 h前后的根系以及叶片组织均形成了明显的聚类差异。

3.2 利用WGCNA构建盐胁迫响应基因共表达网络

WGCNA分析经常被用于各种植物的转录组数据分析, 但多数限制于样本数目或处理时间点有限而采用多胁迫或多组织混合分析的方式。本研究利用盐胁迫处理24 h内的多时间点样本数据, 可以更好地反映基因响应盐胁迫的表达模式, 从而最大程度地排除其他干扰因素,将不同表达模式的基因归类到相应模块。同时, 为了避免组织差异性覆盖胁迫处理的差异性, 本研究分别利用根系以及叶片组织响应差异表达基因进行WGCNA分析,构建组织响应盐胁迫的基因共表达网络。通过分析模块与样本时间点的相关性可以获得组织在不同处理时间点响应盐胁迫的显著模块, 并着重关注早期响应阶段的tan模块和yellow模块以及后期响应阶段的green和red模块。

已有研究表明, 拟南芥中的自噬水平在盐胁迫处理30 min内可达到峰值, 并且这些盐诱导的PCD事件是由Na+特异引起的[52]。本研究在叶片组织中的tan模块显著富集到了自噬反应的调控通路, 该结果与拟南芥中的报道相对一致。油菜在盐胁迫处理0.5 h左右, 体内自噬调节通路就会显著富集, 可能表明此时油菜幼苗的叶片不单受到盐处理的渗透作用, 盐胁迫的离子作用已经可以通过根系传导到叶片, 并引起叶片中盐胁迫特异的离子胁迫反应。植物可以通过多种串扰途径协调各种激素的合成、信号传导以及新陈代谢来应对盐胁迫, 从而建立有效的防御系统[53]。本研究在早期响应模块中显著富集到了脱落酸、乙烯、茉莉酸、水杨酸等激素信号通路, 这可能表明油菜在盐胁迫1 h内就会有多种激素信号参与早期胁迫反应的调控。除激活多种激素信号之外, 盐胁迫早期响应阶段的根系和叶片组织中也显著富集到了植物抗毒素合成、植保素代谢、类黄酮合成以及对细菌的防御反应等通路。已有研究表明, 植物抗毒素主要包括倍半萜、类黄酮和植保素等物质, 其中植保素是拟南芥中一种主要的植物抗毒素, 可以由多种植物病原体诱导产生[54]。这可能表明盐胁迫早期响应阶段可能和病原菌等生物胁迫具有一定的相似性, 同时也会引起机体合成抗毒素等物质进行初级的防御反应。

盐胁迫造成的生理干旱会引起脱落酸积累从而调控气孔关闭。本研究在早期响应模块和后期响应模块均富集到由不同共表达基因参与的脱落酸信号激活通路。这可能表明盐胁迫全程响应过程中, 会有不同的共表达基因在不同的响应阶段参与调控相似的生物学通路。与早期响应阶段相比, 后期响应阶段也具有特异性。后期响应阶段的green模块特异富集到了肌醇合成相关通路。已有研究表明, 原本是受到盐胁迫影响的植物, 在肌醇生产保持不间断时会表现出正常的生长情况, 推测肌醇代谢途径的选择性控制可能是提高耐盐性的方法之一[55]。后期响应阶段的red模块特异富集到了植物型次生细胞壁的发生以及木质素代谢过程的通路。已有研究表明, 植物可以通过半纤维素和木质素沉积来加强次生细胞壁来增加细胞壁厚度, 从而维持盐胁迫的稳态[56]。Red模块也特异富集到了有机酸转运以及苯丙烷合成通路。苯丙烷代谢作为植物重要的次级代谢途径之一, 其代谢产物如木质素以及有机酸等, 在调控植物适应性生长的过程中发挥着重要功能[57]。以上结果表明本研究所挑选的关键模块具有较高的可靠性, 同时也说明利用WGCNA可以构建出具有生物学意义的共表达模块和通路, 这为下一步挖掘核心基因奠定了基础。

在筛选核心基因时, 本研究发现BnaA04G0247200ZS由于具有较高的连通度而被筛选出来, 其对应拟南芥同源基因WRKY33。已有研究表明, 过表达WRKY33可增强拟南芥的耐盐性, 该基因可通过控制根细胞凋亡屏障的形成, 赋予根组织的耐盐能力[58]。油菜中WRKY33可通过增强植保素合成基因的表达增强油菜对菌核病的抗性[59]。本研究分析发现,WRKY33在油菜ZS11中有7个同源拷贝,其中4个拷贝均在盐胁迫早期响应阶段的同一个tan模块中被筛选到, 这说明WGCNA对于挖掘具有多拷贝的核心基因, 也具有较高的适用性和可靠性, 同时推测这些核心转录因子可能也在油菜响应盐胁迫的过程中发挥着重要功能。进一步结合505份盐胁迫处理的油菜群体变异数据库, 发现一些核心转录因子在甘蓝型油菜群体中存在较为丰富的SNPs变异以及单倍型类型。其中, 早期响应基因BnaA08G0015900ZS在盐胁迫处理0.5 h时可上调表达约18倍。其在拟南芥中的同源基因ERF8参与ABA以及免疫信号传导, 并介导植物对病原体的防御反应[60]。单倍型分析发现, 该基因可初步划分为5种单倍型, 且不同单倍型在盐胁迫处理前后的根长方面存在极显著差异。后期响应模块中的核心转录因子基因BnaA05G0156300ZS(NFYA5)在3 h前与对照组基本无差异, 3 h表达量显著上升并在24 h达到最大差异倍数。该基因有3种单倍型, 不同单倍型在幼苗期的植株总干重比(盐胁迫/对照组)方面存在极显著差异。如BnaA07G0311300ZS(ZFHD1)也是后期响应基因, 有2种单倍型, 不同单倍型在幼苗期的植株总干重比(盐胁迫/对照组)方面存在极显著差异。如BnaC01G0510700ZS(NAC047), 有2种单倍型, 不同单倍型在幼苗期的根长比(盐胁迫/对照组)方面存在极显著差异。这些基因在模式植物拟南芥中有所报道, 但在油菜中报道还相对较少, 可为后续甘蓝型油菜耐盐性改良提供重要的候选基因资源。

附表请见网络版: 1) 本刊网站http://zwxb.chinacrops.org/; 2) 中国知网http://www.cnki.net/; 3) 万方数据http://c.wanfangdata.com.cn/Periodical-zuowxb.aspx。