大豆根部水压胁迫耐逆指数遗传体系解析

王吴彬 童 飞 KHAN Mueen Alam 张雅轩 贺建波郝晓帅 邢光南 赵团结 盖钧镒,*

1南京农业大学大豆研究所 / 国家大豆改良中心 / 农业农村部大豆生物学与遗传育种重点实验室(综合) / 作物遗传与种质创新国家重点实验室 / 江苏作物生产协同创新中心, 江苏南京 210095; 2江苏徐淮地区徐州农业科学研究所, 江苏徐州 221131

大豆[Glycine max(L.) Merr.]是重要的植物蛋白质和油料作物, 近年来在世界上迅速发展。地球上41%的陆地面积属于干旱地区[1], 即使湿润与半湿润地区, 由于气候变化问题, 季节性干旱也经常发生[2]。干旱或持续缺水是限制大豆生长发育及产量实现的最重要环境因素。温室及田间研究表明, 水分胁迫能使大豆产量降低 24%～50%[3-4], 严重情况下甚至绝收。因而, 解析大豆耐旱性的遗传机制, 具有重要的育种价值。

干旱可发生在大豆生长发育的不同阶段, 包括苗期、开花期和鼓粒期等。苗期耐旱性鉴定具有试验周期短的特点, 能有效加快耐旱性遗传研究及育种利用进程。耐旱性评价方法有多种, 包括田间鉴定法[5]、干旱棚盆栽法[6]、实验室法[7](水培法和高渗溶液法等)等。耐旱性评价指标包括株高、干重等生长性状和叶片相对含水量和脯氨酸含量等生理生化性状等。

在水分胁迫下, 植物适应逆境, 主要通过积累大量有机或无机物质来调节细胞质的渗透势, 保护酶、蛋白质和生物膜系统[8]。有机物质以脯氨酸、可溶性糖和甜菜碱甘露醇为主, 其中, 脯氨酸是一种溶解度高的小分子渗透物质, 具有良好的水合性, 在植物受到逆境环境时能够大量积累, 可有效减少植物体内水分散失, 实现细胞和组织的保水或吸水, 抑制逆境对植物细胞膜透性的破坏, 同时提高抗氧化酶的活性, 从而减轻逆境造成的伤害[9]。目前, 脯氨酸含量是植物抗逆生理生化研究中通用的一项重要指标[10]。

由于大豆耐旱性易受环境影响, 表型鉴定困难, 导致了其遗传研究滞后于其他性状。截止目前, SoyBase网址(https://www.soybase.org/)共收录88个耐旱性QTL, 其中62个QTL来自于9个不同杂交组合衍生的重组自交系群体, 另外26个QTL是通过不同群体QTL元分析获得, 这些QTL分布在除15号染色体以外的所有染色体上, 其中12号和17号染色体上分布的QTL最多, 可能携带主效耐旱性位点, 具体信息见表1。

为解析大豆苗期耐旱性的遗传基础, 前期研究构建了一套来自蒙8206 (Meng 8202, M)、通山薄皮黄豆甲(Tongshanbopihuangdoujia, T)和正阳白毛平顶(Zhengyangbaimaopingding, Z)的巢式关联作图群体(population from Meng 8202, Tongshanbopihuangdoujia and Zhengyangbaimaopingding, MTZ), 并以根长、茎长、植株长度(根+茎)和干重等生长性状为指标, 检测到 157个耐旱性 QTL[19-20]。在此基础上, 本研究将耐旱指标拓展到生理性状, 在 15% PEG-6000模拟干旱条件下, 以脯氨酸含量为基础的根部水压胁迫耐逆指数为指标, 分别在春季和夏季 2个环境下评价 MTZ群体苗期耐旱性; 然后, 联合 MTZ群体表型与基因型, 利用限制性二阶段多位点全基因组关联分析方法(restricted two-stage multilocus genomewide GWAS, RTM-GWAS), 解析了MTZ群体苗期根部水压胁迫耐逆指数的遗传体系, 以期发掘控制大豆苗期耐旱性的主要位点及其候选基因, 为大豆耐旱性分子育种提供依据。

1 材料与方法

1.1 试验材料

为解析大豆耐旱性的遗传基础, 课题组前期构建了一套以敏旱材料蒙 8206 (M, 熟期组 V)为共同母本, 耐旱材料通山薄皮黄豆甲(T, 熟期组 VI)和正阳白毛平顶(Z, 熟期组 V)为父本的巢式关联作图群体MTZ。MTZ群体共包含409个重组自交系, 其中285份家系来自蒙 8206×通山薄皮黄豆甲(M8206 ×Tongshanbopihuangdoujia, MT), 124份来自蒙8206×正阳白毛平顶(M8206 × Zhengyangbaimaopingding,MZ)。MT和 MZ群体构建方法如下: 以蒙 8206为共同母本, 分别与通山薄皮黄豆甲和正阳白毛平顶杂交, 获得F1, 后经单籽传法自交至7代, 得到重组自交系群体。上述材料均由南京农业大学国家大豆改良中心种质资源库提供。

表1 文献报道的部分大豆耐旱性位点信息Table 1 Information of drought tolerance QTLs reported in the literatures

1.2 试验设计及群体表型评价

试验在南京农业大学国家大豆改良中心温室中进行, 采用完全随机区组设计, 3次重复, 分别在春季(2017年4月, 温度大约6～24℃)和夏季(2017年6月, 温度大约21～32℃) 2个环境下进行。鉴定方法如下: 首先, 剪裁边长约20 cm的正方形滤纸, 选取MTZ群体及其3个亲本各7粒种子, 均匀排列在滤纸上方三分之一处, 匀速卷动滤纸, 保证种子不移位, 卷好后分上中下 3处将纸卷用胶带固定, 然后垂直插入500 mL (直径9 cm, 高11 cm)的塑料杯中,杯中注水; 选择长势一致的大豆植株, 每个处理留2株, 为保证不伤及试验所用幼苗, 将其余幼苗用剪刀剪除; 待植株生长至V1期(第1片三出复叶展开)时, 对幼苗进行含有 15% PEG-6000 (聚乙二醇)的Hoagland营养液胁迫处理, 营养液每周更换 1次;待植株长到 V2期(第 2片三出复叶展开)时, 采取第 2片三出复叶中间小叶进行叶片脯氨酸含量测定。

采用酸性茚三酮显色法测定叶片脯氨酸含量[21-22], 根据耐旱亲本通山薄皮黄豆甲和正阳白毛平顶叶片脯氨酸含量均值计算每份材料的根部水压胁迫耐逆指数(RHSTI), 以增加不同重复和环境间表型数据的可比性, 计算公式为 RHSTIijk=LPCijk/(LPCijT+ LPCijZ)/2, RHSTIijk代表第i个环境、第j次重复中第k个家系的根部水压胁迫耐逆指数,LPCijk表示第i个环境、第j次重复中第k个家系的叶片脯氨酸含量; LPCijT表示第i个环境、第j次重复中耐旱亲本通山薄皮黄豆甲的叶片脯氨酸含量;LPCijT表示第i个环境、第j次重复中耐旱亲本正阳白毛平顶的叶片脯氨酸含量。RHSTI值越大, 表示材料越耐旱。

1.3 MTZ群体基因型鉴定

采用 CTAB法从大豆叶片中提取大豆基因组DNA, 经TaqI酶切后, 选取 400～700 bp之间的DNA片段, 利用Illumina HiSeq 2000测序仪的多元鸟枪法基因分型策略[23], 测定 DNA片段两端碱基序列, 可读取长度90 bp。采用SOAP2软件[24]将测序获得的原始测序reads与Williams 82参考基因组(Glyma.Wm82.a1.v 1.1)进行比对。利用基于贝叶斯模型的realSFS软件[25]进行SNP基因型的鉴定。用以下标准对MTZ群体的SNP进行筛选, 缺失和杂合率≤30%, 最小等位基因频率≥1% (如果存在第 3个或更多个等位基因时, 则用缺失等位基因替换)。通过RAD-seq测序, MZ和MT群体最终分别获得66,677和55,958个SNP, 用于后续划分SNP连锁不平衡区段(SNP linkage disequilibrium block,SNPLDB)标记。此部分工作由深圳华大基因科技有限公司完成。

对MTZ群体SNP数据进行严格过滤, 仅保留基因分型率 80%以上的 SNP标记, 通过并集方式, 把MTZ群体的 SNP数据进行合并, 获得 55,936个SNP。使用fastPHASE软件[26]对缺失SNP进行填补,利用Haploview软件[27]根据MTZ群体的55,936个SNP标记划分单倍型区段, 使用默认设置的置信区间法定义区段, 最大距离和最小 MAF分别设置为200 kb和0.01。位于1个 LD (linkage disequilibrium)区段内紧密连锁的多个 SNP, 整合在一起作为该基因组区域的一段序列, 划分为1个SNPLDB标记。基于以上方法, MTZ群体共检出55,936个SNP组成的6137个SNPLDB标记用于后续分析。该群体基因型已应用于以根长、茎长、植株长度(根+茎)和干重等生长性状为指标的耐旱性遗传研究[19-20]。

1.4 MTZ群体全基因组关联分析

采用贺建波等[28]提出的限制性二阶段多位点全基因组关联分析方法(RTM-GWAS)进行大豆根部水压胁迫耐逆指数遗传解析。所用分子标记是MTZ群体的6137 SNPLDB, 表型数据是2个环境下大豆根部水压胁迫耐逆指数(含每个环境下的 3次重复数据)。RTM-GWAS关联分析分2个阶段进行, 在第1阶段, 基于简单线性模型(simple linear model)进行单位点关联检验(single-locus association test), 对SNPLDB标记进行初步筛选(从6137个SNPLDB中选出 2708个); 在第 2阶段, 对第 1阶段筛选到的2708个显著标记, 利用多位点逐步回归模型进行筛选, 并估计等位变异的效应。在这2个阶段中, 将基于SNPLDB标记的个体间遗传相似系数矩阵的前10个特征向量作为协变量进行群体结构控制。为检测到更多的候选关联位点, 初步筛选标记和多位点逐步回归关联分析时, 显著水平分别为P=0.05和P=0.01, 遗传力控制在37.9%以内。

1.5 MTZ群体QTL-等位变异矩阵构建

将RTM-GWAS关联分析检测到的QTL及等位变异效应构建 QTL-等位变异(QTL-allele)矩阵, 构建方法如下: 以 409个家系为横坐标, 以关联SNPLDB位点为纵坐标, 每个单元格表示某一家系在特定关联标记上的等位变异效应。为简化表示该矩阵, 用不同颜色梯度来表示等位变异效应值的变化趋势, 红色表示增效等位变异, 逐步过渡到蓝色来表示减效等位变异。

1.6 大豆苗期耐旱性候选基因分析

以课题组前期鉴定的耐旱材料光泽黄荚豆为材料, 用Hoagland营养液培养至V1期, 后用含有15%PEG-6000的Hoagland营养液进行胁迫处理12 h, 分别取处理前后的叶片进行转录组测序, 测序深度为5倍。该测序工作委托上海伯豪生物技术有限公司进行。联合关联分析和转录组测序结果, 将位于QTL±500 kb区间内(衰减距离)的差异表达基因认为本研究检测到的耐旱性候选基因。

2 结果与分析

2.1 大豆MTZ群体及其亲本根部水压胁迫耐逆指数表现

在春季和夏季2个环境中, 3个亲本蒙8206、正阳白毛平顶和通山薄皮黄豆甲的根部水压胁迫耐逆指数RHSTI分别为0.35与0.17、1.05与1.04和0.95与 0.96, 共同亲本蒙 8206较其他 2个亲本 RHSTI存在显著差异, 且较低, 说明其耐旱性较弱(表2和图1)。在春、夏2个环境下, MTZ群体RHSTI均表现连续分布, 且具有较高遗传力, 分别为 97.2%和97.9%, 表明PEG模拟干旱条件下大豆RHSTI属于多基因控制的数量性状。在春季环境下, MTZ群体RHSTI平均值为 0.69, 变幅为 0.11～2.94; 在夏季环境下, RHSTI均值为 0.65, 变幅为 0.03～1.93, 表明MTZ群体RHSTI存在广泛变异, 适合大豆耐旱性遗传解析。在 2个不同环境下, 误差变异系数分别为13.2%和 12.1%, 均低于 15.0%, 说明本研究误差控制较好。

表2 MTZ群体及其亲本根部水压胁迫耐逆指数RHSTI的描述性统计Table 2 Descriptive statistics of root hydraulic stress tolerance index (RHSTI) in the MTZ population

对2个不同环境下MTZ群体的RHSTI进行联合分析发现, 群体在 RHSTI上表现连续分布, 均值为 0.68, 变幅为 0.09～2.80, 同样显示了广泛的表型变异。联合方差分析结果显示, 家系以及家系×环境间均存在极显著差异。在RHSTI上, 家系遗传力为37.9%, 家系与环境互作遗传力为 60.1%, 这与叶片脯氨酸含量在春季与夏季间差异较大相符, 可能与不同环境下气温差异较大有关(表3)。

2.2 MTZ群体根部水压胁迫耐逆指数 QTL-等位变异(QTL-allele)体系解析

联合MTZ群体2个环境下的RHSTI表型数据和 6137个 SNPLDB标记基因型数据, 通过RTM-GWAS方法, 在遗传力控制在37.9%以内的情况下, 共检测到了45个与根部水压胁迫耐逆指数显著相关的 QTL, 显著水平(-lgP)变幅为 11.7～73.0,12号染色体上的位点RHSTI12.2显著性最高(图2-A,B); 这些位点中, 有34个位点与环境存在显著互作(QEI, QTL × Environment interaction), 具体结果见表4。检测到的 QTL分布在大豆 Gm01、Gm02、Gm03、Gm04、Gm06、Gm07、Gm08、Gm09、Gm10、Gm11、Gm12、Gm13、Gm14、Gm15、Gm17、Gm18、Gm19和Gm20等18条染色体上, 每条染色上携带有 1个(Gm01、Gm04、Gm08和 Gm20)至 4个(Gm02、Gm09、Gm10和 Gm19) QTL。这些位点总计解释37.58%的表型变异, 单个位点的贡献率变幅为0.46%～3.09%; 其中有 7个 QTL的贡献率≥1%, 归为大贡献 QTL, 应为控制大豆耐旱性的主要位点,累计解释12.00%的表型变异; 另外38个QTL的表型贡献率＜1%, 归为小贡献 QTL, 累计解释 25.58%的表型变异。这些QTL中, 有33个位点与环境存在显著互作(QEI位点), 单个QEI的表型贡献率变幅为0.06～1.70, 互作累计可以解释 12.50%的表型变异;其中3个为大贡献位点, 累积贡献率4.35%, 30个为小贡献位点, 累积贡献率8.15%。根据遗传解析结果,RHSTI的遗传体系可以分解为3部分, 包括: 45个主效QTL, 可以解释37.58%的表型变异; 33个QTL与环境互作效应, 可以解释 12.50%的表型变异; 另外37.9% (h2)+60.11% (互作遗传力)-37.58%-12.50%=47.93%的表型变异, 可能是由未被检测到的微效QTL控制。鉴于本研究检测到的QTL与环境互作效应可能由于春季和夏季不同光温条件引起的, 差异条件较难衡量和重复, 且互作位点包含于主效位点内, 所以后续的QTL-等位变异矩阵和候选基因分析依据主效位点进行。

表3 MTZ群体根部水压胁迫耐逆指数RHSTI两环境联合方差分析Table 3 Joint analysis of variance of root hydraulic stress tolerance index (RHSTI) under two environments in the MTZ population

2.3 MTZ 群体根部水压胁迫耐逆指数QTL-allele矩阵及三亲本后代育种潜力

在 45个根部水压胁迫耐逆指数QTL上, MTZ群体总共携带 123个等位变异, 平均每个位点携带2.73个等位变异, 变幅为 2～3个。在 123个等位变异中, 等位变异效应变幅为-1.911～1.961; 其中, 58个为增效等位变异, 变幅为0.001～1.961, 65个为减效等位变异, 变幅为-1.911～ -0.001 (图2-C)。根据上述结果, 本研究建立了 MTZ群体 45×409的QTL-allele矩阵, 横坐标代表409份半同胞家系, 纵坐标表示45个QTL位点, 每个位点上的等位变异用其效应值表示, 颜色深浅表示等位变异效应大小(图2-D)。随着家系根部水压胁迫耐逆指数的升高, 家系携带增效等位变异有升高的趋势, 并且每个家系均携带有增效的等位变异和减效等位变异, 表明家系间具有较高的重组潜力。

?

根据遗传解析结果, 在本研究检测到的45个主效位点上, 3个亲本通山薄皮黄豆甲、蒙8206和正阳白毛平顶分别携带有24、21和21个增效等位变异, 另外21、24和24个为减效等位变异, 其RHSTI值分别为 1.05、0.26和 0.95。不考虑互作, 只考虑基因加性效应的情况下, 3个亲本的最优重组个体即在45个位点上均携带最优等位变异, 其RHSTI值可以高达18.92, 显示了巨大的超亲潜力(图2-E)。

2.4 大豆根部水压胁迫耐逆指数相关候选基因分析

联合大豆苗期PEG胁迫转录组数据, 在根部水压胁迫耐逆指数相关位点及其附近500 kb区间内共筛选到 38个差异表达基因, 其中 24个基因上调表达, 相对表达量变幅 2.02～4,002,450.00倍; 14个基因下调表达, 相对表达量变幅 0.04～0.50。位于RHSTI9.1位点上的候选基因Glyma.09G051700在PEG处理前后表达量差异高达 4,002,450.00倍, 其编码一个未知功能的蛋白。根据生物学功能, 38个耐旱候选基因可归为7类, 包括ABA响应因子、逆境响应因子、转录因子、发育因子、蛋白代谢因子、未知功能和其他, 分别包含5、17、3、9、8、2和8个候选基因, 其中逆境响应因子和转录因子和发育因子是最大的 3类。表明, 大豆耐旱性是一个复杂性状, 有多个具有不同功能的候选基因共同控制。具体结果见表5。

表5 转录组与GWAS相结合的根部水压胁迫耐逆指数RHSTI候选基因分析Table 5 Analysis the candidate genes of root hydraulic stress tolerance index (RHSTI) combined the results of GWAS and RNA-seq

(续表5)

3 讨论

3.1 大豆耐旱性遗传基础复杂, 不同指标、不同群体具有不同但相似的遗传体系

植物耐旱性评价主要有生理生化指标和生长形态指标两大类。生理生化指标包括脯氨酸含量、丙二醛含量、过氧化物气化酶含量和电导率等, 其是植物受到外界胁迫时的直接反应, 因表型调查困难,群体较大时较难进行; 生长指标较生理指标稳定,表型鉴定成本低, 是植物耐旱性的最终反映。因而,植物耐旱性是复杂的数量性状, 仅凭单个指标很难对其进行综合评价。脯氨酸含量作为植物抗逆生理生化研究中的一项重要指标, 已被用于水稻[29]、小麦[30]、玉米[31]等作物耐逆QTL定位研究中, 然而在大豆中文献报道较少。

为解析大豆耐旱性的遗传基础, 本研究前期构建了由2套重组自交系群体组成的NAM群体MTZ,并以生长指标评价了其耐旱性, 检测到 157个耐旱性QTL[19-20]。本研究将耐旱性指标拓展到生理指标脯氨酸含量, 评价了MTZ群体根部水压胁迫耐逆指数, 检测到 45个 QTL位点, 包含 123个等位变异,这些位点包括7个大贡献QTL和38个小贡献QTL,累计解释37.58%的表型变异。同时发现大部分位点与环境存在互作, 互作可以解释 12.50%的表型变异。与 Mueen等[19-20]利用MTZ群体检测到的相对根长、相对茎长、相对植株长度(根+茎)和相对干重等生长指标耐旱性 QTL结果相整合, 共检测到 18个QTL簇与多个指标相关, 分布在大豆3号、4号、6号、7号、10号、13号、14号、15号、17号、18号、19号和20号染色体上, 其中有3个QTL簇所包含的位点均为大贡献 QTL, 具体见表6。这些QTL簇, 尤其与多个性状相关联的QTL簇应为控制大豆耐旱性的主要染色体区段。上述结果说明, 大豆耐旱性是一个复杂性状, 不同指标具有不同的遗传体系, 共有位点仅占较小部分[16/(45+157)]。

与前人检测到的耐旱性QTL相比, 本研究检测到的45个QTL位点中, 有8个位点与前人耐旱性位点定位结果一致(表4)。RHSTI10.2和RHSTI10.3已被前人在“Benning/PI416937”[13]群体中已检测到;RHSTI14.2和RHSTI15.1位点已被前人在“Benning/PI416937”[13]群体中以叶片萎蔫程度为指标检测到。RHSTI17.1已被前人在“科丰 1号/南农 1138-2”[14]、“KS4895/Jackson”[12]和“93705 KS4895/Jackson”[13]等群体中检测到, 另外 3个位点RHSTI2.2、RHSTI11.1和RHSTI11.2已被前人在多个群体 QTL元分析[18]中检测到。表明, 大豆耐旱性是一个复杂数量性状, 遗传基础复杂, 不同的群体具有不能的遗传基础; 能被多个群体检测到的位点仅占少数,应为控制大豆耐旱性最为主要的QTL。

表6 多个耐旱指标评价MTZ群体检测到的QTL簇(±500 kb)Table 6 QTL clusters of multiple drought tolerance indicators within 500 kb in MTZ population

3.2 大豆耐旱性涉及一套具有不同功能的候选基因体系

借助高密度的分子标记, 本研究可以将根部水压胁迫耐逆指数位点定位到较小的染色体区段内,便于候选基因的预测。结合PEG胁迫下转录组数据,共预测到38个耐旱性候选基因, 可归为7类, 包括ABA响应因子、逆境响应因子、转录因子、发育因子、蛋白代谢、未知功能和其他, 其中逆境响应因子、转录因子、发育因子是最大的3类。38个根部水压胁迫耐逆指数候选基因中, 有 6个基因位于大贡献位点上, 其中RHSTI2.2位点上的候选基因Glyma.02G196500属于 BLH1类转录因子, 参与ABA刺激、耐盐和生长发育等多个生物学过程;RHSTI9.3位点上的候选基因Glyma.09G071600编码中茉莉酸信号蛋白, 参与了ABA信号、缺水响应和花的发育等 27个生物学过程;RHSTI11.2位点上的候选基因Glyma.11G207000参与应对高温和蛋白磷酸化等生物学过程;RHSTI15.1位点上的候选基因Glyma.15G222200参与了植物应对缺水的生物学过程;RHSTI18.1位点上的候选基因Glyma.18G018600参与了应对高温和低磷等生物学过程。表明, 大豆耐旱性是一个复杂性状, 有多个具有不同功能的候选基因共同控制, 位于大贡献位点上的候选基因应为控制大豆耐旱性的主要基因, 具有深入研究价值。

4 结论

PEG模拟干旱条件下, MTZ群体及其亲本在根部水压胁迫耐逆指数上存在显著差异, 且群体表型变异丰富。根部水压胁迫耐逆指数在单个环境下具有较高遗传力, 在不同环境间存在较强互作, 说明MTZ群体根部水压胁迫耐逆指数受遗传和环境的共同影响。利用限制性二阶段多位点全基因组关联分析方法, 对MTZ群体遗传解析发现, 根部水压胁迫耐逆指数遗传基础复杂, 由小贡献和大贡献的位点共同控制, 同时位点与环境间存在显著互作。对定位区间的候选基因进行注释, 共预测到38个候选基因, 涉及ABA响应因子、逆境响应因子、转录因子、发育因子、蛋白代谢、未知功能和其他等不同功能,其中逆境响应因子、转录因子、发育因子是最大的3类, 位于大贡献位点的 6个候选基因参与ABA、逆境和生长发育等多个生物学过程, 具有重要研究价值。上述研究结果表明大豆耐旱性是一个由多位点、多基因控制的复杂数量性状, 且与环境存在互作, 遗传基础复杂。