正常供磷及低磷胁迫下油菜种子植酸浓度与含量的全基因组关联分析

刘海疆，李晓娟，谢毅文，吴海亚，鲁明星，丁广大，徐芳森，石 磊*

（1 华中农业大学作物遗传改良国家重点实验室, 湖北武汉 430070；2 华中农业大学微量元素研究中心, 湖北武汉 430070；3 武穴市农业农村局, 湖北武穴 435400；4 湖北省农业农村厅油菜办公室, 湖北武汉 430070）

油菜是我国也是世界上重要的油料作物，每年种植面积约533万hm2，总产1348多万t[1–2]。菜籽油是国产食用油的第一大来源，在我国食用油市场中具有举足轻重的地位，因此，推动油菜产业发展是我国植物油供给安全的重要保障[2]。磷是植物生长发育过程中所必需的大量营养元素之一[3]。植酸是种子中磷的主要贮存形式。植酸通常以植酸钾、植酸钙和植酸镁等植酸盐的形式积累在种子中。在种子萌发时，植酸水解成磷酸盐和肌醇，并释放出金属离子，为种子萌发提供养分。油菜种子植酸磷占总磷80%左右，缺磷时油菜种子植酸磷和总磷均显著降低[4–6]。菜籽饼是油菜榨油后的副产品，磷含量较高，其中65%是植酸磷[4]。但是，种子和饼粕中植酸并不能被人和其他非反刍类动物消化吸收，造成了磷和其他矿质营养元素的浪费[7–8]。此外，植酸进入水体还可能导致水体富营养化。因此，培育低植酸油菜品种成为油菜育种的一个重要目标，解析其遗传基础对培育低植酸油菜品种具有重要意义。

植酸的合成是由葡萄糖-6-磷酸通过肌醇-3-磷酸合成酶(MIPS)氧化环化生成肌醇-3-磷酸开始的。随后，肌醇单磷酸酶(IMP)通过去磷酸化作用将肌醇-3-磷酸转化为游离肌醇。同时，在肌醇激酶(MIK)的作用下，游离肌醇会再次被磷酸化生成肌醇-3-磷酸。之后，肌醇-3-磷酸通过不依赖脂类途径合成肌醇-1,3,4,5,6-五磷酸，而肌醇通过依赖脂类途径也形成肌醇-1,3,4,5,6-五磷酸。在肌醇多磷酸-2-激酶(IPK1)的催化作用下，肌醇-1,3,4,5,6-五磷酸生成植酸。最终，植酸通过多药耐药相关蛋白(MRP5)运输到蛋白贮存液泡(PSV)中，与金属阳离子螯合，以植酸钙、植酸镁和植酸钾等植酸盐形式储存起来[9]。目前，在拟南芥、水稻、玉米以及油菜中已鉴定出几种参与植酸生物合成的酶，如AtITPK1、AtITPK2、OsMRP5、ZmMRP4、BnPGK2、BnITPK1和BnITPK2[10–13]。除植酸合成和运输相关基因外，一些影响种子植酸浓度的磷转运蛋白基因突变体也被鉴定出来。有研究表明，水稻OsSPDT基因影响磷向籽粒的运输，从而影响种子中植酸的积累[14]。OsSPDT主要在茎节中表达，负责控制磷向籽粒的分配。与野生型相比，OsSPDT突变体种子中植酸磷降低25%～32%，但种子萌发与成熟期产量并未受到影响[14]。

前人利用遗传群体，分别在水稻、拟南芥、绿豆等植物中鉴定到多个与植酸相关的遗传位点(QTL)[15–16]。在水稻中，Stangoulis 等[15]利用粳稻Azucena和籼稻 IR64杂交培育的 DH遗传群体鉴定到2个控制种子植酸含量的QTL，分别位于第5和第12条染色体上，表型解释率分别为24%和15%，其中在第5条染色体上的QTL除与种子植酸含量相关外还和种子中总磷含量密切相关，对总磷含量的表型解释率为20%。在拟南芥中，Bentsink等[16]利用拟南芥Cvi/Ler RIL 群体检测到一个位于第3条染色体上的调控种子植酸含量和叶片植酸含量的多效QTL，该QTL可以分别解释拟南芥种子和叶片植酸含量的61.8%和22.5%的表型变异。相比于连锁分析，关联分析具有定位精度高，无需构建专门的遗传群体，节省时间等优点[17]。Archana等[18]以鹰嘴豆自然群体为材料，通过全基因组关联分析的方法，鉴定到1个与植酸含量相关的位于肌醇单磷酸酶基因的(CaIMP) 5'UTR 区域的简单重复序列标记NCPGR90，该标记简单重复序列长度的变化与鹰嘴豆种子植酸含量的高低相关。Perera等[19]通过GWAS分析鉴定到候选基因肌醇-3-磷酸合酶1(OsINO1)，在植酸积累过程中，与低植酸水稻品种WRC5相比，高植酸水稻品种WRC6的OsINO1的表达量上调了1.8倍，表明INO1在水稻种子植酸的积累过程中起着重要作用。

尽管植酸相关QTL在上述作物中已有报道，但在油菜中鲜有研究。本研究以来源于世界各地的403份甘蓝型油菜为关联分析群体，结合全基因组重测序获得的530余万高质量SNP (single nucleotide polymorphism)标记，对田间试验正常磷和低磷处理油菜种子中的PA_Conc和PA_Cont进行全基因组关联分析，鉴定到与植酸合成相关的SNP位点及其候选基因，开展候选基因的关联分析和单倍型分析，鉴定优异单倍型。研究结果将有助于解析甘蓝型油菜种子植酸含量的遗传机制，为低植酸油菜育种提供理论依据。

1 材料与方法

1.1 试验材料

关联分析群体包括不同来源的油菜品种403份，其中350个半冬性、44个春性、8个冬性和1个未知种性品种。这些品种中，361个来源于中国，21个来自欧洲，8个来自日本，5个来自加拿大，4个来自澳大利亚，3个来自韩国，1个未知。

供试土壤为水稻土，其基本理化性质：pH为5.61、有机质为 7.37 g/kg、全氮为 0.91 g/kg、速效磷为 11.56 mg/kg、土壤全磷为 0.57 g/kg、土壤全钾为13.2 g/kg。

1.2 田间试验与性状调查

上述油菜自然群体于2018年10月至2019年5月在湖北省武穴市梅川镇 (29.85°N，115.55°E)种植，试验设置正常磷(P2O590 kg/hm2)和低磷(不施磷肥，P2O50 kg/hm2)两个处理，每个处理设置3次重复，采用随机区组排列。每个重复每个品种种植4行，每行8个单株，共32株，行距20 cm，株距20 cm，每17个品种为一个小区，共144个小区。播种前，正常磷处理田块的氮、磷、钾和硼肥按108、90、120 和 1.6 kg/hm2作为基肥施入。施用肥料品种分别为尿素 (N 46.6%)、过磷酸钙 (P2O512%)、氯化钾(KCl 60%)和硼砂 (B 15%)。在播种 60 天后，按 72 kg/hm2追施氮素。不施磷处理除不施磷肥外，其余肥料施用与正常磷处理相同，试验田的田间管理措施与当地油菜田常规管理措施一致。待成熟时，选取3株长势一致的植株进行取样，测定种子中植酸浓度并计算植酸含量。

本试验采用铁沉淀法测定油菜种子中植酸浓度[20]，具体步骤如下：用万分之一天平称量5粒种子，记录种子重量，放在1.5 mL离心管中，添加1 mL 0.4 mol/L HCl和 1 粒钨球，利用 TissueLyser ball mill打样机 (Qiagen, Germany)碾磨 3 min，以 15000 r/min的转速离心5 min。离心后移取500 μL的上清液到1.5 mL 离心管。再以 15000 r/min 的转速离心 5 min。最后吸取200 μL的离心后上清液于1.5 mL离心管备用，称为原始样品提取液。同时制备植酸标准溶液，将50% (w/w)植酸溶液稀释得到植酸溶液A (浓度 2.096 mg/mL)，将植酸溶液 A、0.4 mol/L HCl和0.2 mol/L HCl按照 1∶1∶14 的比例配制为 0.131 mg/mL植酸标准溶液，利用0.2 mol/L盐酸按照梯度稀释法，稀释得到4个浓度水平的植酸标准溶液：0.131、0.0655、0.03275、0.016375 mg/mL。向 96 孔PCR 板中，每个孔中添加 15 μL 水和 30 μL 0.2 mol/L HCl，同时添加15 μL的样品提取液。在不同的孔中，添加60 μL植酸标准溶液，每式3份。将96孔PCR板以300 r/min的转速离心20 s，向每个孔中添加 120 μL 0.02% (w/v)硫酸铁铵-0.2 mol/L HCl。96孔PCR板盖帽后，在99℃加热30 min，加热后冰浴冷却10 min。随后样品以3000 g相对离心力离心30 min。离心后将96孔PCR板中各样品和植酸标准溶液转移80 μL到96孔酶标板(平底)中。每个孔中添加 120 μL 显色剂 1% (w/v) 2,2'-联吡啶-1%(v/v)巯基乙酸，置于摇床混合10 min。利用KC4酶标仪 (Bio-Tek Instruments, USA)读取在 519 nm 的吸光度，每个孔读取3次数值，以3次重复的平均值作为吸光度数值。计算后得到种子植酸浓度，以mg/g表示，含量以mg/5粒种子为单位。利用Microsoft Excel 2007 软件初步整理表型数据，计算每个表型性状的平均值、标准差、变异系数，利用R中的psych软件包(https://cran.r-project.org/web/packages/psych/)计算每个性状的峰度与偏度，进行正态分布检验。

1.3 油菜关联分析群体基因型评估

甘蓝型油菜关联分析群体全基因组重测序共获得一千余万的高质量SNP标记[21]。本研究采用MAF> 5%和缺失率<20%为标准对该群体的所有SNP进行过滤，最终得到530万SNP位点用于后续关联分析。对SNP的命名采用以下方法：chr (染色体chromosome的缩写)开头，加上该SNP所在的染色体，再接下划线“__”，再加上该SNP在染色体对应的物理距离。如，“chrA07__7809297”表示位于A07染色体上7809297bp处的一个SNP位点。为了降低群体结构和亲缘关系对关联分析结果的影响，使用 Admixture 软件 (Alexander et al 2009)进行群体结构分析[22]。使用Tassel 5.0软件对该群体进行亲缘关系分析，计算用于矫正关联分析的K矩阵[23]。

1.4 油菜关联分析群体两个磷水平下PA_Conc和PA_Cont的全基因组关联分析

在关联分析中，群体结构(Q)及品种间的亲缘关系(K)会导致关联分析结果存在一定的假阳性[24]。本研究采用TASSEL 5.0软件中的只矫正群体结构(Q)的一般线性模型和同时矫正群体结构(Q)和亲缘关系(K)的混合线性模型来进行关联分析。用于筛选显著关联SNP位点的阈值设为1/n，n为SNP标记数。根据关联分析所计算的所有–lg(P)观察值和期望值，利用R包GGplot2 (https://cran.r-project.org/web/packages/ggplot2/)和 CMplot (https://cran.rproject.org/web/packages/CMplot/)分别绘制Manhattan 图和 Quantile-Quantile散点图 (QQ plot)。

1.5 候选基因预测

根据 Wang等[25]的研究方法，将目标SNP上下游300 kb距离范围内的基因作为初步候选基因，结合拟南芥同源基因的功能注释，筛选确定候选基因。

1.6 候选基因关联分析及单倍型分析

根据 Li等[26]的研究方法，利用软件bcftools对目标候选基因BnaA07.MRP13的CDS区域及其前面2 kb范围的SNP进行提取(http://samtools.github.io/bcftools/bcftools.html)，将所获得的SNP按照最小基因型频率 (minor allele frequency，MAF) > 5% 和缺失率 (missing rate) < 20% 为标准进行过滤，获得高质量SNP标记用于后续的关联分析。使用GLM模型对候选基因内SNP和低磷胁迫下的PA_Cont进行候选基因关联分析，挖掘BnaA07.MRP13内与目标性状相关联的显著SNP，随后利用软件Haploview.4.2对候选基因内的显著SNP进行单倍型分析。

2 结果与分析

2.1 甘蓝型油菜关联群体正常磷和低磷处理种子中植酸浓度和植酸含量的遗传变异

低磷和正常磷处理甘蓝型油菜403个品种种子中PA_Conc和PA_Cont均呈连续性的正态分布，表现出广泛的遗传变异(表1、图1)。其中，正常磷处理PA_Conc变异范围为17.52～81.32 mg/g，均值为36.87 mg/g，变异系数为21.37%；低磷处理PA_Conc变异范围为 17.47 ～49.05 mg/g，均值为 31.91 mg/g，变异系数为24.57%，与正常供磷相比，PA_Conc的均值下降13.5%。正常磷处理每5粒种子PA_Cont变异范围为 0.30～0.99 mg，均值为 0.65 mg，变异系数为18.85% (表1)。低磷处理每5粒种子PA_Cont的变异范围为 0.33～0.85 mg，均值为 0.58 mg，变异系数为18.91% (表1)与正常供磷相比，下降10.8%，PA_Cont的均值下降10.8%。与正常磷处理比较，低磷胁迫PA_Conc和PA_Cont分别下降4.96 mg/g和0.07 mg/5粒种子，但PA_Conc变异系数明显增加(表1、图1)。

表1 正常磷处理及低磷处理甘蓝型油菜关联分析群体种子植酸浓度和含量Table 1 Seed phytate concentration and content in an association panel of Brassica napus under sufficient and deficient phosphorus supplies

图1 正常磷处理和低磷处理甘蓝型油菜关联分析群体种子植酸浓度和含量的频率分布直方图Fig. 1 Frequency distribution of phytate concentration and content of an association panel of Brassica napus under sufficient and deficient phosphorus supplies

2.2 正常磷和低磷处理甘蓝型油菜种子植酸浓度和含量全基因组关联分析

采用GLM和MLM模型对正常磷和低磷处理油菜种子的PA_Conc和PA_Cont进行全基因组关联分析(图2)。其中，利用GLM模型鉴定到正常磷处理与PA_Conc显著关联的SNP位点42个，分别分布在 A01、A02、A03、A04、A05、A06、A09、A10、C01、C02、C04、C05、C06、C07、C08、C09等16条染色体上，解释表型变异范围为8.90%～13.89% (表2)；与 PA_Cont显著关联的SNP位点1个，解释表型变异为7.33% (表2)。低磷处理鉴定到与PA_Conc显著关联的SNP位点6个，其中3个SNP位点分别分布在A03和A10两条染色体上，解释表型变异范围为9.36%～14.28% (表2)；与PA_Cont显著关联SNP位点7个，其中6个SNP位点分别分布在A02和A07两条染色体上，解释表型变异范围为8.37%～12.10% (表2)。

图2 正常磷和低磷处理油菜种子植酸浓度和含量全基因组关联分析的曼哈顿图和QQ图Fig. 2 Manhattan and QQ plots of GWAS on seed phytate concentration and content of Brassica napus under sufficient and deficient P supply

表2 正常磷和低磷处理甘蓝型油菜种子植酸浓度和植酸含量显著关联的SNP位点Table 2 Significant SNPs associated with phytate concentration and content in Brassica napus seed under sufficient and deficient phosphorus supplies

续表2 Table 2 continued

相比于GLM模型，MLM模型对群体结构和亲缘关系均进行了矫正，鉴定到的显著关联的位点较少，但MLM模型鉴定到的23个显著SNP位点均与GLM共定位(表2)。其中，正常磷处理与PA_Conc显著关联的SNP位点20个，分别分布在A03、A05、A06、A09、C01、C02、C04、C05、C06、C07等10条染色体上，解释表型变异范围为9.60%～13.89% (表2)。MLM模型没有检测到与正常磷PA_Cont显著关联的SNP位点，低磷处理鉴定到1个位于A03染色体上与PA_Conc显著关联的SNP位点，该位点可解释的表型变异为14.28% (表2)；2个与PA_Cont显著关联的SNP位点，解释的表型变异分别为11.18%和11.22% (表2)。

2.3 候选基因预测

基于甘蓝型油菜Darmor-bzh参考基因组序列，对显著关联SNP位点前后300 kb范围内的基因进行分析，同时结合拟南芥同源基因的基因功能注释，最终确定14个参与植酸合成和转运途径的基因作为候选基因(表3)。由GLM和MLM模型共同检测到的植酸浓度和含量显著关联的23个SNP位点的候选基因有结合盒转运蛋白MRP13(BnaA07g07310D，BnaA07g07320D)、MRP12(BnaA07g07330D)，肌醇1,3,4磷酸-5/6-激酶蛋白(BnaA10g17710D)，磷脂酰肌醇磷酸酶、磷脂酶C (BnaC09g33980D，BnaC 09g34000D，BnaC09g34010D)，磷酸甘油酸变位酶家族蛋白(BnaUnng04100D)以及磷转运相关的蛋白(BnaC06g37190D)等的编码基因。此外，由GLM模型单独鉴定到的油菜种子植酸浓度和含量显著关联SNP位点的候选基因有液泡转运蛋白(BnaC07g 30170D)、磷脂酰肌醇-3,4-激酶(BnaC07g30200D)、磷脂酰肌醇单磷酸(BnaC04g24700D)以及类硫酸盐转运蛋白(BnaA03g04410D)的编码基因(表3)。

2.4 BnaA07.MRP13候选基因关联分析及单倍型分析

在A07染色体上由GLM和MLM模型同时鉴定到低磷处理植酸含量显著关联SNP位点chrA07__7809297，确定多药耐药相关蛋白的基因MRP13(BnaA07.MRP13)为该位点的候选基因(表3)。BnaA07.MRP13的CDS区域及其前面2 kb范围内共有19个SNP位点，将这19个SNP位点与低磷处理植酸含量进行候选基因关联分析。结果表明，共有4个SNP位点与低磷胁迫下植酸含量显著相关，且这4个SNP表现出很强的连锁不平衡(图3a)。对这4个SNP位点进行单倍型分析，共分成2种单倍型：高植酸单倍型“BnaA07.MRP13ContHap1”(TTTA)和低植酸单倍型“BnaA07.MRP13ContHap2”(CCCT) (图3b)。高植酸单倍型“BnaA07.MRP13 ContHap1”具有218个油菜品种，其5粒种子的平均植酸含量为0.62 mg，低植酸单倍型“BnaA07.MRP13ContHap2”具有72个油菜品种，其5粒种子的平均植酸含量为0.56 mg (其他油菜品种SNP缺失，不计算在内)。含有“BnaA07.MRP13ContHap1”的品种在低磷胁迫下的植酸含量要显著低于含有“BnaA07.MRP13ContHap2”的品种 (P= 0.0064)(图3b)。

图3 BnaA07.MRP13的候选基因关联分析(a)和单倍型分析(b)Fig. 3 Candidate genes association analysis (a) and haplotype analysis (b) of BnaA07.MRP13

3 讨论

3.1 甘蓝型油菜种子植酸浓度和植酸含量显著关联的SNPs及其候选基因

磷肥的施用量显著影响作物种子植酸浓度和含量[27–28]。研究表明，种子中植酸浓度与磷肥的施用量、磷肥在土壤中的有效性及溶解度显著相关[29]。本研究正常磷与低磷处理分别检测到42和6个与种子植酸浓度显著关联的SNP位点，分别检测到1和7个与种子植酸含量显著关联的SNP位点(表2)。但是，在正常磷和低磷胁迫下没有鉴定到相同的SNP位点。在 A07染色体上鉴定到1个低磷处理PA_Cont显著关联的SNP位点和由 3个邻近的 SNP位点组成的显著 SNP 簇(表2)。这3个SNP位点中，chrA07__7809297和chrA07__7809300同时被GLM和MLM模型检测到，预测BnaA07.MRP12(BnaA07g07330D) 和BnaA07.MRP13(BnaA07g 07310D，BnaA07g0732D)为该位点的候选基因(表3)。拟南芥AtMRP5编码ABC 转运子家族的多药耐药相关蛋白，负责将合成的植酸运输至蛋白贮存液泡[30]。与野生型比较，水稻OsMRP5突变体植株籽粒植酸浓度降低35.8%～71.9%，无机磷浓度显著提高，且无机磷的增加量要高于植酸浓度下降量[31]。BnaA07 g07330D、BnaA07g07310D和BnaA07g0732D均属于MRP家族基因，很有可能参与油菜种子中植酸的运输(表3)。低磷处理PA_Cont显著关联SNP鉴定到硫酸盐转运蛋白家族SULTR4;1(BnaA03g04410D)(表3)。水稻OsSPDT基因编码硫酸盐转运蛋白，与野生型比较其突变体种子植酸浓度降低20%，而产量无显著降低，这有利于减少水稻种子中植酸磷通过人和动物流向环境，减轻环境污染[14]。此外，低磷处理PA_Cont显著关联SNP位点还鉴定到液泡膜铁转运家族蛋白(VIT) (BnaC07g30170D)、肌醇或磷脂酰肌醇磷酸酶(BnaA02g34000D)以及参与将肌醇-6-磷酸(InsP6)合成肌醇-7-磷酸(InsP7)的激酶(BnaUnng04100D) (表3)。这些候选基因及其功能需要通过CRISPR-Cas9等分子生物学技术进一步验证。

植酸含量是种子植酸浓度和种子质量的乘积，因此，种子植酸含量与植酸浓度和种子质量均相关。如果种子植酸含量的差异主要由植酸浓度决定，两者的显著SNP位点可能一样；如果与两者均有关，或者主要由种子质量决定，两者的显著相关的SNP位点可能就不一样(表2)。

3.2 甘蓝型油菜低植酸品种的选育

作物种子中的磷是全球磷循环的主要驱动力，通过降低种子中植酸的浓度可以减少有限磷资源的浪费[32]。Matthus等[33]利用高效液相色谱法研究油菜种子的磷组成，发现油菜种子中的植酸浓度要显著高于其他作物。研究人员利用CRISPR-Cas9技术构建了甘蓝型油菜BnITPK1和BnITPK4基因的双突变体，结果表明，与野生型比较，双突变体种子植酸浓度显著降低27.2%～35.3%；此外，双突变体种子萌发、千粒重、株高、含油量以及产量与野生型无显著差异[12–13]。这表明，通过编辑油菜中BnITPK1和BnITPK4可以实现低植酸育种。在本研究中，通过低磷胁迫油菜种子植酸含量鉴定到植酸合成最后一步的MRP家族基因BnaA07.MRP13。BnaA07.MRP13基因的候选基因关联分析及单倍型分析结果表明，与低磷胁迫下植酸含量显著相关的4个SNP构成2种单倍型：“BnaA07.MRP13ContHap1”(TTTA)和“BnaA07.MRP13ContHap 2”(CCCT)。其中，“BnaA07.MRP13ContHap2”单倍型属于低植酸单倍型，对于后续油菜低植酸育种研究具有重要意义。本研究油菜自然群体正常磷处理PA_Conc的变异幅度为17.52～81.32 mg/g，低磷处理PA_Conc变异范围为17.47～49.05 mg/g。这些低植酸品种也可以作为低植酸育种的种质资源。

4 结论

甘蓝型油菜关联分析群体正常磷和低磷处理种子植酸浓度和植酸含量均具有广泛的遗传变异。全基因组关联分析共检测到56个与种子植酸浓度和植酸含量显著关联的 SNP 位点，筛选和鉴定了14个与油菜种子植酸合成相关的候选基因，通过候选基因关联分析以及单倍型分析鉴定了油菜种子高和低植酸单倍型“BnaA07.MRP13ContHap1”(TTTA)和“BnaA07.MRP13ContHap2”(CCCT)，这为解析油菜种子植酸浓度和含量的遗传机制奠定了基础，为低植酸品种的培育提供了优异种质和基因资源。但是，这些候选基因的功能及其决定高和低植酸单倍型的分子机制还不清楚，在今后的研究中，可利用CRISPR-Cas9 等技术创制候选基因的突变体，深入开展基因功能研究，明确这些基因影响种子植酸含量的生物学机制。