锌胁迫下甘蓝型油菜发芽期下胚轴长的全基因组关联分析

魏丽娟 申树林 黄小虎 马国强 王曦彤 杨怡玲 李洹东 王书贤 朱美晨 唐章林 卢 坤 李加纳 曲存民

西南大学农学与生物科技学院 / 油菜工程研究中心 / 西南大学现代农业科学研究院, 重庆 400715

近年来, 土壤重金属污染问题愈发严重, 对人类身体健康和生态系统产生重大影响[1]。锌是植物生长发育所必需的微量元素之一, 适量锌可促进植物生长发育, 增强植物对不良环境的适应性[2-3]。然而过度使用不仅会造成严重的环境污染, 还会降低土壤微生物的活性, 减缓土壤有机质的分解, 进而影响植物根的生长和光合作用[4-6]。油菜被认为是修复重金属(如Cu、Cd、Zn)污染土壤的理想植物之一[7-10], 对锌具有高耐受性, 在锌胁迫下仍可获得较高的生物产量[7,11]。曹春信等[12]通过研究油菜地上部器官对锌富集的特征发现, 锌在茎秆中的含量大于籽粒。挖掘耐锌性关键位点和候选基因对耐锌性油菜品种的培育具有重要意义, 全基因组关联分析被广泛应用于作物重金属抗性位点的检测。Zhang等[13]利用全基因组关联分析方法, 在水稻中共检测到 31个耐锌QTL位点。Chen等[14]利用BrassicaIllumina 60K SNP芯片进行全基因组关联分析, 检测出25个甘蓝型油菜耐镉QTL位点, 并找到NRAMP6、IRT1、CAD1和PCS2四个候选基因。Zhang等[15]检测到4个甘蓝型油菜耐镉位点, 并鉴定出一些候选基因,如HIPP27、EXPB4、EMB1793和CDSP3。曲存民等[16]检测到 70个甘蓝型油菜耐砷性的显著关联位点, 并筛选出PHT3;3和PHT1;9等与甘蓝型油菜砷离子吸收转运相关的候选基因。超富集植物对重金属的积累和耐性机制研究也取得了一定进展, 许多蛋白参与锌的吸收、转运和贮存, 包括CDF (Cation Diffusion Facilitator)蛋白家族[17]、NRAMP蛋白家族[18]、ZIP蛋白家族[19]、HMA 蛋白家族(P-type ATPases)[20]、YSL (yellow stripe-like)蛋白家族[21]、ABC (ATP-binding cassette)转运蛋白家族[2,22]。但关于甘蓝型油菜耐锌性QTL定位和候选基因鉴定并未见相关报道。

为挖掘甘蓝型油菜与耐锌性状相关联的关键位点及候选基因, 本研究采用GWAS方法获得与下胚轴长显著关联的SNP标记, 结合转录组数据获得锌胁迫前后差异表达基因, 进一步确定影响下胚轴长的关键候选基因。本研究不仅为甘蓝型油菜耐锌性种质资源的鉴定提供分子标记, 而且可以为耐锌性油菜资源利用和新品种培育提供理论依据。

1 材料与方法

1.1 试验材料

供试材料为国内外收集的甘蓝型油菜自交系品种140份, 其中国内材料123份, 主要来自重庆、四川、湖北、湖南等地, 其余17份材料来自加拿大和德国等国家, 所有材料由西南大学重庆市油菜工程技术研究中心收集并提供。

1.2 试验处理与性状分析

从 140份研究材料中随机抽取10份材料(附表1), 分别用 0、5、10、15、20、30和 50 mg L-1浓度的 Zn(NO3)2进行处理, 比较分析发现, 处理浓度为30 mg L-1时材料间长势差异较大, 确定为最佳处理浓度(表1)。试验分为对照组和处理组, 重复3次, 挑选大小均匀、饱满的种子各 100粒, 分别播种于对照纯水和 30 mg L-1的 Zn(NO3)2溶液处理的培养盘(长 37.0 cm × 宽 24.5 cm × 高 6.5 cm, 体积为 4.5 L),覆盖保鲜膜。生长条件为昼夜温度为 25℃, 光照/黑暗时间为 16 h/8 h, 光照强度为 100 μmol m–2s–1, 相对湿度为60%[23]。播种7 d后, 选取长势一致的5株幼苗, 分别测量对照和锌胁迫下的下胚轴长(hypocotyl length, HL)。将处理和对照的相对值作为衡量锌胁迫对甘蓝型油菜生长的影响[24], 记为相对下胚轴长(relative hypocotyl length, RHL)。计算公式: 相对性状值(%)=处理组测定值/对照组测定值×100%[16]。利用 Microsoft Excel 2013绘制正态分布图, 利用 SPSS统计分析软件进行偏度和峰度分析。

1.3 基因型分析及SNP标记定位

利用芸薹属60K SNP芯片, 参照Qu等[25-26]方法, 对 140份材料进行基因型分析, 去掉最小基因型频率(minor allele frequency, MAF)低于 0.05和SNP得率(call frequency)小于 80%的标记。为确定SNP在染色体上的位置, 利用探针序列与基因组进行BLASTN比对, E-value设定为1e−5。去掉Bit-Score相同的SNP, 最终获得30,288个SNP标记用于群体性状的关联分析。

1.4 群体结构分析与亲缘关系

利用STRUCTURE软件进行甘蓝型油菜群体结构分析[27], 组群数目(K)设定为1～10,K值运行5次(number of iterations), 将MCMC (Markov Chain Monte Carlo)不作数迭代(length of burn-in period)和不作数迭代后的MCMC全部设为100,000次, 在混合模型下进行运算。最后用运算得到的后验概率值和 2个连续的后验概率值的变化速率(ΔK)来确定群体中存在的亚群数目[28]。利用SPAGeDi v1.4软件进行亲缘关系(relative kinship)分析, 并计算亲缘关系的矩阵[29]。

1.5 全基因组关联分析及候选基因鉴定

利用 TASSEL 5.0软件[30]中的一般线性模型(General Linear Model, GLM)和混合线性模型(Mixed Linear Model, MLM), 结合 SNP芯片基因型数据,以Q值和亲缘关系K值作为协变量, 进行全基因组关联分析, 确定关联位点和标记对表型变异的解释率(R2)。本研究采用了6种模型进行关联分析, 包括GLM中的naive、Q和PCA模型, MLM中的K、Q+K和 PCA+K模型。根据每个 SNP位点的-log10(P)观察值和期望值, 利用 R软件(https://www.r-project.org/)绘制Quantile-Quantile散点图(QQ plot)。曼哈顿图(Manhattan plot)利用 R语言程序包的 mrMLM(Multi-Locus Random-SNP-Effect Mixed Linear Model)进行绘制[31]。与性状显著关联的SNP位点的阈值为P< 3.30×10-5, 计算方法为: 1/N(N=30,288)[16]。利用 TASSEL 5.0软件[30]计算 linkage disequilibrium (LD)在油菜不同染色上的分布情况,确定不同染色体的LD衰减距离, 根据LD确定候选基因区间。

1.6 甘蓝型油菜锌胁迫转录组测序

取甘蓝型油菜自交系07191锌胁迫处理前和处理7 d后的幼苗, 于液氮中保存, 参照曲存民等[16]方法提取总RNA, 利用Illumina Hiseq 2500进行转录组测序, 以甘蓝型油菜为参考基因组[32], 利用Cuffdiff筛选差异表达基因, 筛选标准为FDR < 0.01,log2|(fold change)|>1。上调和下调的差异表达基因GO 富集分析采用基迪奥云平台进行分析(https://www.omicshare.com/)。

2 结果与分析

2.1 甘蓝型油菜发芽期锌胁迫最佳浓度确定

为了确定油菜发芽期锌胁迫的最佳处理浓度,本研究随机挑选 10个材料, 分别利用不同浓度的Zn(NO3)2溶液和纯水为对照进行处理, 发现处理浓度为30 mg L-1时, 这些材料下胚轴长度明显受到抑制, 且材料间变异较大(表 1)。因此, 在发芽期筛选耐锌油菜的最佳浓度为30 mg L-1。

2.2 甘蓝型油菜锌胁迫性状表型变异分析

对140份油菜正常和锌胁迫条件下的下胚轴长度进行统计发现, 在正常条件下, 甘蓝型油菜下胚轴长度为 3.66～7.99 cm, 变异系数为 15.8%; 锌胁迫下, 下胚轴长度为 2.28～6.47 cm, 变异系数为18.8%; 相对下胚轴长的变异范围为 0.38～1.08, 均值为 0.73, 变异系数为 17.6%, 表现出广泛的表型变异(表2)。锌胁迫下, 大部分材料的下胚轴长度都较正常生长状态的下胚轴短, 只有5个半冬性材料例外(甘油 5号、SWU94、中油 821、川油 18和wx10213), 表明30 mg L-1锌胁迫严重抑制了甘蓝型油菜下胚轴伸长。相对下胚轴长呈连续性正态分布, 说明其受多基因控制, 表型为典型的数量性状遗传特点(图1)。

表2 锌胁迫下甘蓝型油菜发芽期性状统计分析Table 2 Statistical analysis of the traits at germination stage under zinc stress

2.3 群体结构和亲缘关系分析

利用30,288个SNP标记进行关联群体的群体结构和亲缘关系分析, 当K=2时, ΔK值最大(图2-A),因此, 该群体的亚群数为2, 生成的最优群体结构的Q矩阵见图2-B。140份甘蓝型材料中, 89%材料间亲缘关系小于0.1 (图2-C), 表明该群体内材料间的亲缘关系较远。

2.4 耐锌胁迫下相对下胚轴长度的全基因组关联分析

对相对下胚轴长度的QQ图进行分析(图3), 结果表明, K+Q模型检测到的标记的P值较其他模型更接近期望值, 其降低了假阳性, 关联分析结果更可靠。全基因组关联分析共检测到 8个显著关联的SNP位点(图 3和表 3), 分别位于 A02、A07、C04和 C06染色体上, 单个 SNP解释的表型变异为22.0%～33.2%。为了确定显著关联 SNP位点的候选基因区间, 本研究对这些关联位点的染色体进行LD分析, 利用R2值估计LD, A02和A07染色体的衰减距离为300 kb, C04染色体衰减距离为700 kb, C06染色体的衰减距离为1 Mb (图4-A)。另外, A07染色体上22.2 Mb～23.9 Mb与C06染色体34.9 Mb～37.2 Mb位于共线性区间内(图4-B), 说明找到的显著关联的SNP位点位于这2个共线性区间内。

表3 锌胁迫下甘蓝型油菜发芽期相对下胚轴长显著关联的SNP位点Table 3 Significant SNP locus associated with the relative hypocotyl length at germination stage under zinc stress in B. napus

2.5 锌胁迫下油菜幼苗差异表达基因

转录组分析结果表明, 锌处理油菜幼苗 7 d后,与对照相比, 上调基因 2733个, 下调基因 3652个,其中差异转录因子/基因 616个, 主要包括 bHLH(77)、ERF (73)、NAC (65)、MYB (60)、WRKY (55)和 C2H2 (33)家族(图5)。ERF、NAC和WRKY家族基因上调数目较多, 而bHLH、MYB和C2H2家族基因下调基因数目较多。本研究对上调和下调基因分别进行 GO功能富集分析, 上调基因主要富集在细胞外区域和细胞膜, 分子功能是硫酸盐跨膜转运功能, 主要参与氧化还原反应、离子转运、胁迫反应、防御反应和硫化合物转运(图6-A)。下调基因主要富集在细胞外区域、质外体和细胞膜区域, 主要的分子功能是抗氧化活性和过氧化物酶活性, 主要参与抗氧化反应、氧化还原反应和细胞壁生物合成(图 6-B)。

2.6 GWAS与转录组结合鉴定共有基因

根据各染色体的 LD确定候选区间, 参照油菜基因组注释信息[32], 共找到 19个与锌胁迫相关的候选基因(表 4), 这些基因主要包括编码锌指蛋白家族成员(B-box型和 ZFP1)、谷胱甘肽转移酶GSTU21、过氧化物酶家族蛋白、ABC和 MFS转运蛋白及细胞壁相关激酶蛋白, 并找到一些重要的转录因子, 如MYB、bHLH、NAC、ERF和WRKY,大部分基因(除GSTU21、ZFP1)。

3 讨论

表4 甘蓝型油菜锌胁迫下相关性状候选基因Table 4 Candidate genes of the traits correlated with zinc stress treatment in B. napus

重金属污染可归纳为 2类, 一类是一些具有显著生物毒性的重金属物, 如镉(Cd)、砷(As)等, 另一类为超过一定量后具有生物毒性的一些重金属, 如锌(Zn)、铜(Cu)等[33]。近年来, 土壤重金属污染越来越引起人们关注[34]。在土壤中过量的锌富积会抑制植物生长, 导致营养失衡、萎黄叶片, 以及光合作用受损[5-6]。研究表明, 甘蓝型油菜对重金属锌具有较强的耐受性, 在较高浓度下仍可获得较高的生物产量[2-5]。然而关于油菜响应锌胁迫的分子机制尚未见报道。在甘蓝型油菜中, 通过GWAS分析与基因组信息相结合已成为挖掘油菜数量性状基因座的常规手段, 并获得了一些与非生物胁迫相关的候选基因[14-16,35-37]。因此, 本研究通过对锌胁迫下140份甘蓝型油菜发芽期的相对下胚轴长进行了全基因组关联分析, 共检测到8个显著关联的SNP位点, A02、C04和C06染色体上显著关联的SNP位点分别位于22.0、4.4和36.8 Mb, A07染色体上显著关联的SNP位点位于22.0 Mb和23.5 Mb。曲存民等[16]将砷胁迫下甘蓝型油菜的RHL定位在A02: 19.2 Mb～22.0 Mb、A07: 19.2 Mb～22.5 Mb、A07: 23.1 Mb～23.5 Mb、C04:4.3 Mb～4.5 Mb 和 C06: 31.6 Mb～36.8 Mb, 与本研究定位结果一致, 这些位点具有一因多效性, 可能在不同重金属胁迫过程中都起着重要作用。

在植物中, 重金属可以引起许多代谢过程发生变化, 这些变化主要是通过氧化还原反应引起的, 同时造成抗氧化酶的活化、细胞膜通透性改变和离子的置换等反应[38]。本研究发现锌胁迫后, 上调的差异表达基因主要富集在细胞膜, 参与氧化还原反应、离子转运、防御反应和硫化合物转运。本研究结合全基因组关联分析和转录组测序, 找到 19个与锌胁迫相关的候选基因, 包括转录因子/基因, 如WRKY57、ORA47和BnaA07g28000D(NAC),WRKY57可以通过清除活性氧提高植物抗旱和抗盐能力[39-40]; 植物受到伤害和水分胁迫时,ORA47可通过调控脱落酸和茉莉酸信号途径来提高抗性[41]。本研究中WRKY57和ORA47在锌胁迫后显著上调表达, 可能在锌胁迫过程中起着重要作用, 但具体的调控机制还需进一步研究。

本研究还找到了 B-box型和锌指蛋白 ZFP1,锌指蛋白中的锌指结构可以结合锌离子形成稳定的指状结构, 编码锌指蛋白的众多基因已经被证明在生物和非生物胁迫过程中起着重要作用[42-43]。Sun等[44]研究表明, 菜豆 MTF通过激活色氨酸生物合成调控烟草耐镉性, 锌指基序是其重要的结合金属反应的元件。拟南芥锌指蛋白ZAT6通过谷胱甘肽途径提高镉抗性[45], 而锌指蛋白 ZAT11负调控拟南芥镍抗性[46]。本研究找到的B-box型锌指蛋白基因BnaA02g30040D在锌胁迫后表达量上调, 而锌指蛋白ZFP1基因表达量下调, 它们可能通过不同的机制来响应胁迫。另外, 本研究找到 2个编码Na+/H+逆向转运蛋白的基因NHX6(BnaA07g35030D和BnaC06g39970D), 锌胁迫后其表达量上调, Yang等[47]研究表明, 浮萍NHX1基因通过减少Cd2+流入和增加Cd2+外排提高重金属镉抗性。本研究中NHX6表达上调, 可能增加了Zn2+的外排, 从而提高锌抗性。

重金属处理下植物细胞产生活性氧造成细胞氧化损伤, 为了抵抗活性氧对细胞的损伤, 植物会产生抗氧化酶及非酶抗氧化剂清除自由基[48]。锌处理下小麦过氧化氢含量、过氧化物酶及过氧化氢酶活性均显著提高, 增强了小麦对锌的耐受性[49]。本研究找到一些与氧化还原反应相关的基因, 如编码过氧化物酶家族蛋白和谷胱甘肽转移酶基因GSTU21。GST可使GSH与内源电子化合物结合而产生脱毒[50],重金属(Cu和 Al)胁迫下, 拟南芥GST基因表达上调, 超表达水稻GSTU30基因提高了拟南芥重金属耐性[51-52]。但是本研究中,GSTU21基因在锌胁迫后表达量下调, 可能是由于锌处理植株后, 重金属造成的氧化胁迫反应超过了油菜自身的协调能力, 抗氧化能力受到明显抑制。

另外, 本研究还定位到一些编码 ABC和 MSF转运蛋白的基因, 其功能是将重金属转运到液泡内,在重金属转运过程中起着重要作用[53-54], 这些转运蛋白基因表达量上调, 可能通过相似的途径增强植株对锌的抗性。拟南芥铝胁迫后, 可以诱导WAK1基因的表达[55];WALK4基因启动子的突变导致水稻对钾、钠、铜和锌超敏感,WALK4突变导致水稻茎中锌含量下降[56];WAKL10基因表达量上调, 其可能在锌胁迫过程中起着重要作用。本研究为深入研究甘蓝型油菜响应锌胁迫的分子机制奠定了基础。

4 结论

全基因组关联分析共鉴定出8个与锌胁迫显著关联的SNP位点, 结合转录组测序结果, 共找到19个候选基因, 可能在锌胁迫过程中起着重要作用。

附表1 140份甘蓝型油菜材料信息Table S1 Detailed information of the 140 B. napus accessions

(续附表 1)

(续附表 1)