甘蓝型油菜角果长度性状的全基因组关联分析

孙程明 陈 松 彭 琦 张 维 易 斌 张洁夫,* 傅廷栋

1 江苏省农业科学院经济作物研究所 / 农业部长江下游棉花与油菜重点实验室 / 江苏省现代作物生产协同创新中心,江苏南京 210014；2 华中农业大学植物科学技术学院 / 作物遗传改良国家重点实验室,湖北武汉 430070

油菜是我国重要的油料作物,种植面积约690万公顷,产油量占国产植物油总量的55%以上,是我国最大的食用植物油来源[1]。作为食用油消费大国,我国食用油自给率仅35%左右,因此,提高产量成为油菜育种最迫切的目标。角果长度是油菜重要的产量相关性状,解析其遗传基础对提高油菜产量具有重要意义。

角果是油菜重要的光合产物贮存器官,适度增加角果长度,可增加角果库容量,有利于灌浆物质向籽粒中运输,每角粒数和籽粒大小也随之提高[2]。此外,角果还是油菜成熟后期重要的光合器官,在角果完全伸展后,角果皮面积占全株光合面积的60.6%[3]；对角果遮光处理导致油菜的千粒重降低54.7%、籽粒产量降低62.1%、含油量降低44.2%[3]。冷锁虎等[2]研究表明,角果长度每增加1 cm,群体中角果皮面积占总光合面积的比例增加2.64%。因此适度增加角果长度有利于扩大角果库容量,增加全株光合面积,从而提高油菜的籽粒产量。

油菜角果长度是受多基因调控的数量性状,众多研究者采用遗传差异较大的双亲构建连锁群体,进行角果长度的QTL 定位。Udall 等[4]定位的角果长度QTL 分布在A09、C02 和C08 连锁群,其中A09连锁群的QTLln9.5贡献率为25.3%[4]；Yang 等[5]定位到的QTL 分布在A01、A02、A06、A07、A09、C02、C03 和C04 连锁群,其中A09 连锁群的QTLqSLA9.3能解释65.6%的表型变异；漆丽萍[6]在A01、A07、A09 和A10 连锁群检测到11 个QTL,其中贡献率最大的2 个QTL 均位于A09,在3 个环境中分别解释38.3%～60.0%和 11.2%～16.5%的表型变异；Wang 等[7]在A05、A06、C01、C05 和C06 检测到11 个QTL,解释3.4%～14.36%的表型变异；Fu 等[8]在A01、A05、A09 和C08 连锁群上检测到20 个角果长度QTL,效应最大的2 个QTL均位于A09；Yang等[9]在A05、A09、C02、C08 和C09 检测到7 个角果长度QTL,位于C09 的主效位点解释26.51%的表型变异。

近年来,随着高密度SNP 基因分型芯片和全基因组测序等技术的发展,全基因组关联分析(genome-wide association study,GWAS)已广泛应用于玉米[10-11]、水稻[12-13]、油菜[14-15]等作物复杂性状的遗传结构解析。本研究利用60K SNP 芯片对群体进行基因型分析,并对角果长度进行全基因组关联分析,旨在挖掘显著的SNP 位点,分析其候选基因,为油菜角果长度的遗传改良奠定基础。

1 材料与方法

1.1 试验材料

用于关联分析的496 份甘蓝型油菜包括地方品种、育成品种及高世代育种材料。其中国内资源444份,主要来自湖北、重庆、江苏、湖南、四川、陕西等油菜主产省市；国外资源52 份,主要来自德国、瑞典、朝鲜、加拿大等国家(附表1)。所有材料均由华中农业大学国家油菜工程技术研究中心提供。

1.2 田间试验与表型调查

2013年、2014年于江苏南京(13NJ 和14NJ),2015年、2016年于江苏泰州(15TZ 和16TZ),采用完全随机区组设计种植自然群体,在南京设置3 次重复,在泰州设置2 个重复。单个材料每重复种2行,每行15 株,行宽1.5 m,行距40 cm,行长1.5 m。每年10月上旬大田直播,田间管理按当地常规方式进行。在油菜成熟期,取每小区8 株长势一致的植株,选主花序中部10 个角果测量角果长度。利用R 软件对各环境的角果长度表型进行统计分析和相关性分析,并计算广义遗传力和最佳线性无偏预测值(best linear unbiased prediction,BLUP)[16-17]。

1.3 基因型检测与分析

利用Illumina 60K SNP 芯片分析基因型。在Genome Studio 软件中,剔除Call Freq ＜ 0.75,Minor Freq ＜ 0.05,AA,BB frequencies ＜ 0.03 及GenTrain Scores ＜ 0.50 的SNP,基因型杂合的SNP 按缺失值处理。将过滤后的33,218 个SNP 序列比对至油菜Darmor-bzh基因组,E-value 阈值设为10-12,保留19,167 个单拷贝、位置清楚的SNP 标记用于后续分析。

1.4 全基因组关联分析

用Structure 2.3.3 计算群体结构Q 矩阵[18],用SPAGeDi 计算材料间的亲缘关系K 矩阵[19]。将基因型、表型、Q 矩阵和K 矩阵导入Tassel 4.0 后,以Q矩阵和Q+K 矩阵作协变量,分别进行基于一般线性模型(General Linear Model,GLM)和混合线性模型(Mixed Linear Model,MLM)的全基因组关联分析[20]。关联分析的显著性阈值(Bonferroni threshold)定为1/总标记数,即-lg(P)= 4.28。当1 Mb 区间内存在多个显著SNP 时,如果两两间的r2≥ 0.1,则将这些SNP 归为一个关联位点,以最小P值的SNP作为代表。用R 软件包qqman 绘制曼哈顿图和QQ(Quantile-Quantile)图[21]。用R 软件中逐步回归分析的lm 功能分析关联位点解释的总表型变异。首先,挑选出关联区间的lead SNP,将其中一种等位基因型标记成1,另外一种标记成-1,杂合和缺失位点标记成0；利用“stepAIC”函数,通过选择最小的AIC(Akaike information criterion)信息统计量,逐步增加或剔除lead SNP,最终确定最适宜的模型,模型的R2即为总表型变异。

1.5 已报道QTL 的信息整理

目前有2 篇角果长度的关联分析研究[22-23],2 个群体分别包含300 个和157 个株系,整理其显著位点的信息并与本研究的结果比较。角果长度的连锁分析报道很多,但多数研究的标记为SSR 和AFLP等传统标记,数目少且无法准确确定其物理位置。漆丽萍[6]和Yang 等[9]分别用SNP 芯片对DH 群体进行基因型分析,本研究将其角果长度 QTL 两侧的SNP 标记比对至油菜Darmor-bzh参考基因组以获得QTL 的物理区间。此外,Fu 等[8]利用混池测序将一个角果长度主效QTL 定位至油菜A09 染色体27.21～ 28.26 Mb 区段。

1.6 候选基因挖掘

以r2= 0.1 作为衰减阈值,用Tassel 计算显著关联位点的LD 衰减距离作为其置信区间。提取所有显著位点置信区间内的基因CDS 序列,与拟南芥的基因CDS 进行BLAST 比对,E-value 阈值设为10-10,以相似度最高的拟南芥基因信息来注释油菜基因。

2 结果与分析

2.1 表型统计分析

496 份油菜资源的角果长度在4 个环境中具有广泛的变异,单个环境的角果长度极差范围是5.93～ 12.52 cm。单个环境表型平均值范围是(5.40±1.01)～(6.39±1.25)cm,变异系数范围是0.13～ 0.20(表1和图1)。由表2可知,4 个环境的角果长度存在极显著相关(P≤ 0.001),表明本研究考察的表型数据具有较高的可靠性和重复性。根据R 软件包lem4 的计算,角果长度在4 个环境的广义遗传力为86.0%。

表1 关联群体角果长度性状的统计分析 Table1 Statistical analysis of silique length of the association panel

图1 关联群体在4 个环境的角果长度分布 Fig.1 Distribution of silique length of the association panel in four environments

表2 4 个环境的角果长度相关系数 Table2 Correlation coefficients of silique length among four environments

2.2 角果长度全基因组关联分析

利用MLM 模型,当阈值为-lg(P)= 4.28 时,检测到15 个显著SNP。合并存在连锁不平衡的SNP(r2≥ 0.1),得到7 个显著位点,分布在A07、A09、C02、C04 和C09 染色体(表3和图2)。效应最大的位点Bn-A09-p29991443位于A09染色体[-lg(P)= 13.59],对表型变异的贡献率为13.89%；携带其等位基因型A 的材料(48 份)角果长度均值为6.62 cm,比携带基因型T 的材料(389 份)角果长度均值增加0.89 cm。商品种如华油6 号、华油12 号、中双5 号、中双7 号、中双11 号、宁油6 号、湘油13 号等,地方品种包括南川长角、湖北白花油菜等均携带该位点的优异等位基因型(附表2)。其余位点的-lg(P)值范围是4.29～ 5.80,可解释3.21%～5.10%的表型变异。利用一般线性模型计算,7 个位点联合解释25.01%的表型变异。与单环境的关联结果比较,5 个BLUP 位点在单环境中被重复检测到,3 个位点在2 个环境中得到验证,Bn-A09-p29991443 位点在3 个环境中被重复检测到。

通过GLM 模型共检测到136 个显著SNP,合并存在连锁不平衡的SNP 后得到25 个显著位点,分布在A02、A03、A06、A07、A09、A10、C02、C04、C07、C08 和C09 染色体(表4和图2)。效应最大的位点和 MLM 相同,在 GLM 模型中 Bn-A09- p29991443 的-lg(P)= 16.27,对表型变异的贡献率为12.86%。其余的位点-lg(P)值范围是4.31～9.38,可解释2.91%～7.76%的表型变异。利用一般线性模型计算,25 个位点联合解释41.77%的表型变异。与单环境的关联结果比较,23 个BLUP 位点在单环境中被重复检测到,其中8、6 和5 个位点分别在2、3和4 个环境中得到验证。比较2 个关联模型的结果,5 个MLM 位点与GLM 重叠,合并重叠的位点共得到27 个位点。利用一般线性模型计算,27 个位点联合解释44.04%的表型变异。

图2 油菜角果长度全基因组关联分析(BLUP)Fig.2 Genome-wide association study of rapeseed silique length(BLUP)

表3 MLM 角果长度显著关联位点(BLUP)Table3 Significant GWAS loci of silique length in MLM(BLUP)

表4 GLM 角果长度显著关联位点(BLUP)Table4 Significant GWAS loci of silique length in GLM(BLUP)

2.3 与前人报道QTL 的比较

综合MLM 和GLM 的结果,与已报道的角果长度QTL 结果比较,本研究发现的7 个位点与前人研究结果一致(表3和表4)。周庆红等[22]、漆丽萍[6]、Fu 等[8]、Dong 等[23]通过不同的定位策略检测到效应最大的位点均落在A09 染色体27.26～28.61 Mb区段。本研究中效应最大的位点Bn-A09-p29991443也同样落在该区间。位于 A07 的位点 Bn-A07- p18319586 物理位置是20.22 Mb,与Dong 等[23]检测到的角果长度显著位点(范围为19.76～20.57 Mb)重叠,紧邻漆丽萍[6]检测到的QTLqSL.A7.2(范围为21.04～22.52 Mb)。位于C08 的位点Bn-scaff_16361_ 1-p241830 物理位置是27.75 Mb,与Yang 等[9]检测到的QTLcqSL.C08a(范围为25.92～27.79 Mb)重叠。此外,还有4 个位点与前人检测的QTL 重叠或紧邻。这些结果证实了本研究结果的可靠性。

此外,综合MLM 和GLM 的结果,本研究共检测到20 个新位点。多个位点效应较高,如位于A07的Bn-A07-p8572785 表型贡献率为9.38%,位于C04的Bn-scaff_20817_1-p60579 表型贡献率为7.98%。多个位点在多环境中被检测到,可信度高,例如位于C04 染色体的位点Bn-scaff_27914_1-p9919 和Bn- scaff_20817_1-p60579 及位于C07 染色体的位点Bn- scaff_16069_1-p3731985,这3 个位点在4 个环境中均被检测到。这些位点值得进一步验证和研究。

2.4 候选基因挖掘

基于中双11 号基因组注释信息,在Bn-A09- p29991443 位点下游153 kb 和184 kb 处分别找到油菜角果长度已知基因ARF18和BnaA9.CYP78A9(表5)。ARF18编码一个生长素响应因子,通过生长素信号途径调节角果皮的细胞生长[24]。BnaA9.CYP78A9编码一个细胞色素P450 单加氧酶,通过调节细胞的伸长影响角果长度[25]。基于Darmor-bzh基因组注释信息,在Bn-A09-p3051349 位点上游252 kb 处找到一个候选基因BnaA09g05500,该基因的拟南芥同源基因FUL编码一个MADS-box 转录因子,该基因突变导致植株角果无法伸长,果实不能正常开裂[26]。在Bn-scaff_ 20817_1-p60579 位点上游293 kb 处找到一个候选基因BnaC04g50960,该基因的拟南芥同源基因EOD3编码一个细胞色素P450 单加氧酶,其功能缺失突变体角果长度缩短、种子变小[27]。在Bn-scaff_16069_1-p1668600位点下游490 kb 处找到一个候选基因BnaC07g36530,该基因在拟南芥中的同源基因DOF4.4属于Dof 转录因子家族,RNAi 株系角果增长、产量增多[28]。此外,本研究还找到别的候选基因,如拟南芥已知基因GID1b和ga20ox1在油菜中的同源拷贝[29-30]。

表5 角果长度关联位点候选基因信息 Table5 Information of candidate genes of silique length GWAS loci

3 讨论

角果长度是油菜重要的农艺性状,适度增加角果长度有利于扩大角果库容量,增加植株光合面积,从而提高油菜的籽粒产量。因此,众多研究者围绕角果长度开展了一系列的定位工作,在油菜19 条染色体上均检测到 QTL,对表型变异的贡献率为3.4%～65.5%[4-9,22,23]。本研究对496 份油菜资源的角果长度进行4 个环境的考察,通过全基因组关联分析挖掘到27 个显著位点,分布在A02、A03、A06、A07、A09、A10、C02、C04、C07、C08 和C09 共11 条染色体,对表型变异的贡献率为 2.91%～ 13.89%。多个研究者在A09 染色体27.26～28.61 Mb区段检测到主效QTL,最高可解释超过60%的表型 变异[5,6,8,22-23]。本研究中效应最大的位点 Bn-A09- p29991443 同样落在该区间。此外,本研究还检测到另外6 个与前人结果共定位的位点及20 个新位点,所有位点联合解释44.04%的表型变异。可将携带不同有利等位基因的材料杂交,利用与上述QTL 紧密连锁的分子标记进行辅助选择,聚合更多有利等位基因,培育角果长度理想的油菜新品种。

目前油菜中已克隆的角果长度基因较少,Liu等[24]和Shi 等[25]分别报道在A09 染色体克隆了角果长度基因ARF18和BnaA9.CYP78A9。2 个基因均落在本研究最显著的位点Bn-A09-p29991443 所在区段,相距28.1 kb,处于紧密连锁状态。可能由于2个基因效应叠加,致使该区间的显著性很高。此外,在对应的 C08 同源区段检测到微弱的关联信号(GLM,-lg(P)= 3.72),表明ARF18和BnaA9.CYP78A9位于C08 的同源拷贝,对角果长度的表型贡献很小,与Liu 等[24]和Shi 等[25]的研究结果一致。Bn-A07-p8572785 是用MLM 和GLM 检测到显著性较高的位点,但在附近未找到角果长度已知基因。分析其附近基因的注释信息,找到一个候选基因BnaA07g10150,该基因位于Bn-A07-p8572785 上游339 kb 处,其拟南芥同源基因PIN7编码一个生长素运输基因。生长素调节细胞的伸长生长,ARF18通过生长素途径调节角果长度,因此BnaA07g10150的功能值得进一步研究。Bn-scaff_20817_1-p60579 是GLM 模型中显著性较高的位点,根据基因注释信息在该位点附近发现重要候选基因BnaC04g50960,其拟南芥同源基因EOD3同时调控角果长度和籽粒大小[27],后续研究可以通过转基因或CRISPR 敲除验证该基因的功能。

4 结论

用MLM 检测到7 个位点,联合解释25.01%的表型变异；用GLM 检测到25 个位点,联合解释41.77%的表型变异。合并共同的位点后得到27 个位点,其中7 个位点与前人QTL 共定位,其余20 个是新位点。效应最大的位点Bn-A09-p29991443 位于A09 染色体,在其附近检测到油菜已克隆的角果长度基因ARF18和BnaA9.CYP78A9。此外,在其他多个位点附近发现拟南芥已知角果长度基因GID1b、FUL、EOD3、DOF4.4和GA20ox1的同源拷贝。

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