大豆巢式关联作图(NAM)群体构建及花色和种皮色遗传分析

宋 健 熊亚俊 陈伊洁 徐瑞新 刘康林 郭庆元洪慧龙 高华伟 谷勇哲 张丽娟 郭 勇 阎 哲 刘章雄关荣霞 李英慧 王晓波 郭兵福 孙如建 闫 龙 王好让姬月梅 常汝镇 王 俊,6,* 邱丽娟,*

1 长江大学, 湖北荆州 434025; 2 中国农业科学院作物科学研究所 / 农作物基因资源与遗传改良国家重大科学工程 / 农业农村部北京大豆生物学重点实验室, 北京 100081; 3 安徽农业大学农学院, 安徽合肥 230036; 4 江西省农业科学院作物研究所, 江西南昌 330200; 5 呼伦贝尔市农业科学研究所, 内蒙古呼伦贝尔 162650; 6 农业农村部长江中游作物绿色高效生产重点实验室(部省共建) /长江大学农学院, 湖北荆州 434025; 7 河北省农林科学院粮油作物研究所, 河北石家庄 050035; 8 江苏徐怀地区徐州农业科学研究所,江苏徐州 210031; 9 宁夏农林科学院作物研究所, 宁夏银川 750021

目前常见的遗传解析方法主要包括基于双亲群体的连锁分析(linkage analysis, LA)和基于自然群体的全基因组关联分析(genome wide association study,GWAS)[1]。连锁分析利用了双亲杂交和后代自交过程中创造的人工重组事件, 然而受限于双亲群体大小以及自交过程中发生的有限重组事件, 通常定位分辨率较低[2]。利用自然群体的GWAS 分析充分利用了祖先重组(或历史重组), 结合高通量标记分型技术, 显著提高了定位效率[3]。然而, 由于自然种群中遗传相关性和种群结构的存在, 可导致假阳性关联[4], 此外GWAS 分析对低频等位基因和微效位点检测效果较差[5]。为补足遗传分析及关联分析的各自缺陷, 康奈尔大学Ed. S. Buckler 课题组从302 份玉米自交系中选择具有丰富遗传多样性的25 份材料作为父本与公共母本B73 进行杂交, 获得了包含5000 株系的玉米NAM 群体[4,6]。利用该群体对开花期、大/小斑病抗性、叶型结构、茎秆强度、株高、干旱抗性、籽粒品质(淀粉、蛋白和油分)、碳氮代谢等多性状进行了分析, 取得了较好研究结果[7-13]。除玉米外, 大麦、水稻、高粱、大豆、小麦、油菜、花生、蚕豆、豌豆等作物也相继构建了NAM 群体[14-22]。

大豆是世界上最重要的粮油作物之一, 多数具有经济价值的性状属于数量性状, 如农艺性状和产量性状等。这些性状受多因素影响, 如不同遗传位点、遗传位点间互作、环境以及环境与遗传位点间互作, 遗传调控机制十分复杂。为了更好解析复杂性状的遗传机制, 美国农业部农业研究局(USDAARS)从美国大豆成熟区组I～V 的120 个品种中选择了40 个高产、抗旱且祖先血缘多样的品种或种质(其中3 份来源于中国), 与IA3023 进行杂交获得包含5600 RIL 株系的NAM 群体[23]。利用该群体对产量稳定性、产量与环境互作、产量相关农艺性状(成熟期、株高、倒伏性、百粒重)、光合作用与水分利用效率、复杂性状相互作用、开花期、籽粒成分、品质、疫霉病抗性等性状进行了遗传解析和基因组选择分析[24-28]。印度大豆研究中心利用JS335 为共同亲本与20 个不同大豆品种进行杂交, 构建包含900个RIL 株系的大豆NAM 群体, 旨在对产量、涝渍抗性、干旱抗性、黄色花叶病毒抗性、锈病抗性和早熟等性状进行遗传解析[21]。国内, 盖钧镒院士团队利用4 个亲本(临河、正阳、通山与WSB)与蒙8206杂交最终获得包含623 个RIL 株系的NAM 群体,并开发了适用于NAM 群体研究的RTM-GWAS 分析流程[29-30]。

颜色性状是最基本的生物学性状和农艺性状,是大豆遗传、分子和生化过程研究中最常用的表型标记之一, 被广泛用来快速描述品种特点[31-32]。同时花色和种皮色与驯化密切相关, 野生大豆多为紫花和黑色种皮, 而栽培大豆则多为白花和黄色种皮[33-34]。近年来, 关于大豆种皮色和花色的相关研究日益增多, 尤其是在大豆的色素组成、遗传、基因调控等都开展了深入研究。大豆花色由6 个独立的基因座(W1、W2、W3、W4、Wm和Wp)控制[35-36]。在这些位点中,W1编码类黄酮3’5’-羟化酶(F3’5’H)[37-38],W3和W4编码二氢黄酮醇-4-还原酶(DFR)[35,39],W1和W3/W4功能缺失均会产生白花,w4等位基因背景下的W3位点突变(即w3w4)产生了近乎纯白色的花,而W3w4基因型产生仅在近花萼部位为紫色的白花[35,38,40-41];Wp编码黄酮3-羟化酶(F3H), 突变可以将紫花转变成粉红色[42];Wm编码黄酮合成酶(FLS),隐性突变等位基因会产生蓝色花色;W2编码转录因子[43]。除此之外, 研究还发现MYB 转录因子的调节基因GmPH4和GmPH5影响液泡酸化, 进而影响花色[44]。大豆种皮色多样性高, 具有黄、黑、绿、褐、红、虎斑、鞍挂等多种不同颜色, 其中大豆的野生祖先多为黑色, 栽培大豆多为黄色[32-33]。目前已经克隆的种皮色基因至少有7 个位点(I、T、W1、R、O、K1和G)参与控制性状的形成, 主要参与编码类黄酮途径关键酶或转录因子, 其中I编码查尔酮合成酶CHS 基因[45-46],T编码黄酮3’羟化酶F3’H 基因(Glyma.06G202300)[47]、K1编码AGO5 基因(Glyma.11G190900)[46]、R位点编码MYB 转录因子[48-50]、G位点编码CAAX 基因(Glyma.01G198500)[51-54], 此外研究发现转录因子基因GmTT2A、GmTT2B、GmMYB5A表达可以引起原花青素或花青素积累[55]。

中国是大豆发源地, 拥有超过30,000 份大豆种质资源, 具有高度遗传多样性[56-57]。为了更好发挥种质资源优势, 更深入解析大豆不同性状遗传调控机制, 本研究报道了一个由36 个亲本构建的大豆NAM 群体, 并以种皮色(数量性状)和花色(质量性状)为例, 深度挖掘了两性状的遗传调控位点。该NAM群体的构建为国内大豆遗传研究和育种实践提供了良好的基础材料。

1 材料与方法

1.1 材料

本研究所用大豆种质资源材料均来源于国家大豆种质资源库(表1)。NAM 群体构建采用中豆41 (中国农业科学院油料作物研究所培育)为公共母本, 其余35 份来自不同省份的大豆资源材料为父本, 杂交后F1代自交, F2代开始采用单粒传法, 繁殖至F11代,不同世代种植地点见构建进程表(表2)。

表2 大豆NAM 群体构建进程Table 2 Construction process of soybean NAM population

1.2 花色与种皮色表型采集与分析

根据《大豆种质资源描述规范和数据标准》, 在大豆开花期观察开花颜色, 记录并进行统计分析。种皮色表型采集则是对NAM 群体中父本种皮色为黑色、绿色、褐色、虎斑等与母本黄色种皮不同的9 个群体, 包括N001、N002、N003、N007、N014、N016、N017、N032 和N035 收获的种子进行拍照, 在同一属性设置(P 档, ISO300, 曝光补偿-2)环境下,将相机(Canon EOS 550D)固定在70 cm 高度照相,利用Tomato Analyzer-color Test4.0 软件[58]读取大豆种子照片, 软件的颜色测试模块记录所选对象的每个像素的RGB 值并计算, 将种皮颜色转化为反应颜色的色彩属性值, 包括红(R)、绿(G)、蓝(B)、光度(Luminosity)、明暗度(L)、红绿色(a)、黄蓝色(b)、色调(Hue)和色彩饱和度(Chroma) 9 个色彩属性[59-60]。3个重复, 取平均值。利用R 语言包pastecs 和Hmisc对9 个色彩属性数据进行描述性统计和相关性分析,其中采用Pearson 相关系数法进行相关性分析。

1.3 NAM 群体基因型鉴定、聚类及PCA 分析

NAM 群体于2022年在安徽省宿州市 埇桥区安徽农业大学试验站(33.69°N, 117.10°E)种植, 于V5时期取幼嫩叶片提取基因组DNA, 利用200K SNP芯片完成基因型鉴定[61], 由北京康普森生物技术有限公司完成DNA 提取及芯片基因型鉴定。

聚类分析采用最大似然法, NAM 群体亲本和NAM 群体采用的基因型过滤标准为完整度＞99%和次要等位基因频率＞10%, 分别得到72,219 个SNP 和49,476 个SNP。SNP 数据先通过Tassel 转换为phylip格式, 然后分别采用 FastTree[62]和 IQtree2[63]构建NAM 群体及亲本最大似然树。NAM 群体最大似然树可靠性检验采用 Bootstrap 法(1000)。最后采用itol (https://itol.embl.de/)和Adobe illustrator进行可视化。

采用GCTA (1.93.2beta)软件[64]默认参数分析完成PCA 分析, 基因型采用全部芯片鉴定的NAM 群体个体, 过滤标准为完整度＞80%和次要等位基因频率＞5%。前3 个主成分(PC1～PC3)的二维和三维可视化利用ggplot2 和scatterplot3d 包绘制。

1.4 GWAS 分析

在PCA 分析所用基因型数据基础上, 分别筛选出父母本在花色和种皮色存在差异的分离群体的基因型数据, 利用Beagle[65]进行基因型填补, 然后利用4 种方法对种皮色和花色进行QTL 定位, 包括限制性两阶段多位点全基因组关联分析(Restricted twostage multilocus genome-wide association analysis,RTM-GWAS)[30]、三方差组分多位点随机SNP 效应混合线性模型(3 variance-component Multi-locus random-SNP-effect mixed linear model, 3VmrMLM)[66]、TASSEL (Trait analysis by association, evolution and linkage)[67]和EMMAX (Efficient mixed-model association eXpedited)[68]。4 种方法所用表型数据一致, 与基因型数据对应。

RTM-GWAS 构建了 SNP 连锁不平衡块(SNP linkage disequilibrium blocks, SNPLDB), 并以显著水平0.01 为阈值, 采用多位点模型检测QTL。类似地, 3VmrMLM[66]利用多位点模型进行QTL 定位,似然函数比值对数值(logarithm of odds, LOD)阈值设为3, 其通过R 语言包3VmrMLM 实现, method参数选用“Single_env”。基于混合线性模型(Mixed linear model, MLM)模型, 利用TASSEL 5.0[67]软件和EMMAX[68]软件中emmax-kin-inter64 命令分别计算亲缘关系矩阵, TASSEL 利用PCA 结果中的前5 个主成分和亲缘关系矩阵作为协变量, EMMAX 仅以亲缘关系矩阵为协变量, 二者均并以-log10(0.1/m)(m为SNP 总数量)检测显著位点并采用R 语言包CMplot 绘制曼哈顿图和QQ 图。

1.5 候选基因分析

对关联分析的区间范围, 参照GlycinemaxWm82.a2.v1 基因组(http://phytozome.jgi.doe.gov/pz/portal.html#!info?alias=Org_Gmax)提取基因, 基因功能根据SoyBase (http://soybase.org/)和Phytozome(http://phytozome.org/)注释结果, 同时结合文献检索基因的功能注释。根据已公开可用的大豆转录组数据集[69], 分析候选基因在W82 不同组织的特异性表达。表达量数据使用log10(x+1)进行归一化, 其中x为基因表达水平, 通过计算每千个碱基的转录每百万映射读取的fragments (Fragments Per Kilobase of transcript per Million reads mapped, FPKM)获得。聚类热图则利用R 语言包pheatmap 绘制。

2 结果与分析

2.1 大豆NAM 群体亲本及群体组成

本研究基于大豆种质产量表现、地理远缘与遗传多样性等标准, 从来源于全国19 个省市的大豆种质中选择36 份材料作为NAM 群体亲本, 其中中豆41 作为公共母本, 其余35 份种质作为父本。NAM群体亲本涵盖了北方春、黄淮春、黄淮夏、南方春和南方夏等5 种不同的生态类型, 其中14 份来源于黄淮地区夏大豆, 占38.9%, 其次是北方春大豆占比27.8 %, 来源于黄淮的春大豆最少, 仅占5.6%(表1 和图1-A)。从种质类型来看, NAM 亲本以选育品种为主, 包含27 份, 约占亲本总数的77%, 其余为地方品种。地方品种亲本主要为南方夏大豆(4),其次为黄淮夏(3), 北方春最少(2)。

图1 大豆NAM 群体亲本地理来源、分类与群体大小Fig. 1 Geographic origination, classification of NAM parents and population size

从性状分布来看, NAM 亲本蛋白含量范围为30.9%～48.2%, 油分含量范围为16.1%～25.2%, 蛋脂总和在56.0%～67.0%之间, 百粒重介于9.3～30.2 g 之间(图1-E)。亲本中高蛋白(＞45%)品种8 个, 其中包含1 个地方品种(邳县四粒糙), 其余为育成品种; 高油(＞21.5%)品种3 个, 均为育成品种(表1 和图2-A)。除此之外, 亲本中还包含特异种质1 份, 为β 亚基(7S 贮藏蛋白)低含量材料。父本中种皮色包括黄、绿、黑、褐、虎斑等5 种颜色, 分别包含26、2、3、3 和1 份品种。对花色而言, NAM 群体父本中白花品种12 个, 紫花品种23 个。

图2 NAM 群体遗传结构分析Fig. 2 Genetic structure of NAM population

大豆NAM 群体由35 个RIL 群体组成, 不同RIL 原始群体大小并不一致, 分布在46～576 之间(图1-F), 平均群体大小约200 RIL, 合计7010 株系。为了不同群体间的相对平衡, 本研究从不同NAM 亚群体中分别随机选取44～72 RIL 株系(小于72 株系个体全部选取)组成包含2420 RIL 株系的NAM 群体。

2.2 NAM 群体遗传结构分析

为明确NAM 群体亲本及NAM 群体遗传结构,本研究对所有NAM 群体亲本及选择的RIL 株系进行了基因型鉴定, 共鉴定标记159,035 个SNP, 经过完整度和次要等位基因型频率过滤后保留112,912个SNP。利用NAM 群体基因型分别构建NAM 群体亲本与NAM 群体最大似然树发现, 黄淮地区NAM群体亲本与东北育成品种以及来自南方的NAM 群体亲本可分别聚类, 而来自东北地区的亲本则与南方大豆亲缘关系较远(图 2-A), 可能与生产上黄淮地区大豆育种利用亲本材料来源广泛而南方地区大豆育种极少直接利用东北地区大豆作为育种亲本有关。鉴定的全部NAM 群体个体进行聚类分析发现, 来源于不同RIL 群体的个体基本聚类在一起,不同RIL 群体间结构清晰(图2-B)。主成分分析发现, 第1 至第3 主成分分别解释基因型的4.53%、4.24%和2.63%, 散点图(图2-C, D)显示, 来自不同NAM 群体RIL 亚群体的个体分布相对集中, 具有明显的群体结构特征, 表明大豆NAM 群体遗传结构清晰可溯。

2.3 NAM 群体花色、种皮色表型分析

对NAM 群体表型鉴定分析发现, 在花色分离群体包含1531 RIL 株系, 其中白花816 份, 紫花663份, 白紫花52 份。种皮色分离群体包含561 株系,Red、Green、Blue、Luminosity、L、a、b、Hue、Chroma 9 个颜色属性的分布范围分别为 56.74～224.83、66.14～225.37、1.41～206.22、60.34～203.82、26.85～88.92、-17.16～16.28、-9.08～39.43、26.88～242.16、2.47～42.64 (表3)。对颜色属性进行相关性分析发现, Red 与Green、Blue、Luminosity、L、b、Chroma 之间的相关系数在0.727～0.960 之间, 均呈极显著正相关(P＜0.01), Green 与 Red、Blue、Luminosity、L、b、Chroma 之间相关系数在0.683～0.994 之间, 均呈极显著正相关, Blue 则只与Red、Green、Luminosity、L 呈极显著正相关, a 和Hue 与其他属性相关度低于0.6 (表4)。

表3 种皮色9 个属性表型数据统计分析Table 3 Statistical analysis of 9 different attributes of seed coat color

表4 种皮色9 种不同属性值相关性分析Table 4 Correlation analysis of 9 different attributes of seed coat colors

2.4 NAM 群体全基因组关联分析

利用质控后得到的124,029 个高质量的SNP 基因型数据, 使用EMMAX 方法进行种皮色和花色全基因组关联分析, 以5.09 (-log10(P))为阈值筛选显著关联的SNP 位点。结果发现, 花色检测到1 个主要位点, 位于13 号染色体, 显著相关的SNP 主要分布在 Chr13: 17,215,500..17,425,536 之间,220.04 kb, 最显著的SNP 为Chr13_17,315,536。种皮色检测到8 个有显著SNP 连续分布的QTL 位点,分别位于1、5、6、8、9 和12 号6 条染色体上。为了进一步验证关联位点的可靠性, 本研究采用Tassel、RTM-GWAS 和3vmrMLM 等3 种不同方法分别进行全基因组关联分析, 发现上述8 个关联位点均在3 种以上方法中关联到, 同时在多种不同颜色属性表型中重复鉴定(表5、图3 和图4)。其中qSCC8-1、qSCC9-1在9 种不同颜色属性中均关联到,关联定位重叠区间分别位于 Chr08: 7,088,758..9,481,774 和Chr09: 45,305,565..45,822,016 区间。qSCC5-1和qSCC6-1, 除Blue 以外其余8 个颜色属性均可关联到, 分别位于 Chr05: 36,175,387..36,290,717 和 Chr06: 18,409,621..19,314,852 。qSCC12-1位于Chr12: 4,401,461..4,591,803, 被7个颜色属性关联到。而在 Chr06: 3,940,810,608..42,012,070 区间也关联到一个位点qSCC6-2, 与5个颜色属性性状相关,qSCC1-1仅与a、Hue 相关。同时, Tassel、RTM、3vmrMLM 3 种方法又共定位到一个位于 10 号染色体的位点, 重叠区间位于Chr10: 43,550,478..44,132,066, 关联到6 个颜色属性。RTM、3vmrMLM 共同定位到3 个位点, 分别位于11 号和13 号染色体(表5)。

图3 基于EMMAX 的花色和种皮色全基因组关联分析Manhattan 图和QQ 图Fig. 3 Manhattan-plots and QQ-plots for genome-wide association analysis of flower color and seed coat color by EMMAX

(图4)

2.5 关联区间候选基因分析

对花色位点候选基因进行分析,qFC13-1位点与W1位点重合, 推断花色候选基因为编码F3’5’H 的Glyma.13G072100, 该基因位于Chr13: 17,312,471..17,317,198, 基因上芯片检测到 2 个 SNP, Chr13:17,315,536: A-G, Chr13: 17,316,431: C-A, 其中Chr13: 17,315,536 显著性最强, A 等位基因型在亲本品种中为12 个白花4 个紫花, G 等位基因型19 个紫花1 个白花, 与花色显著相关(P＜0.001)。

对通过不同方法共定位到的12 个种皮色位点进行分析, 其中9 个位点由3 种不同方法共定位, 4个已克隆经典位点和5 个新位点, 已经克隆的位点候选基因具体可见表6。为进一步挖掘与种皮颜色相关的候选基因, 分析发现5 个新位点共包含138个候选基因, 对其中功能可能与种皮颜色相关的基因[33,36,46,55,70-72], 比如转移酶基因包括UDP-糖基转移酶、谷胱甘肽转移酶、ABC 转移酶, 类黄酮可能相关基因功能细胞P450 基因、查尔酮异构酶基因、转录因子MYB 结构域、bHLH 结构域、WD40 结构域、锌指蛋白等基因进行表达谱数据分析, 发现部分在种皮特异性高表达或者特异性不表达的基因10个, 可能为候选基因, 具体表达模式见图5。其中qSCC1-1位点候选基因Glyma.01G075200在种皮高表达, 编码 60S 核糖体蛋白, MFPT repeat-like;qSCC5-1位点有 2 个基因在种皮中相对高表达,Glyma.05G171700、Glyma.05G172400(RCC1, 染色体凝聚调节因子)、Glyma.05G173300编码羰基水解酶在种皮和种脐中不表达, 同时该位点附近存在Glyma.05G153200(Chr05: 34,687,008.. 34,693,243)编码CHS 查尔酮合成酶;qSCC6-2位点有2 个基因在种皮中高表达,Glyma.06G242200编码MRKY 转录因子、Glyma.06G247200编码 W40 转录因子;qSCC12-1有3 个基因在种皮中特异高表达,Glyma.12G060700编码锌指蛋白转录因子在种皮和种脐中高表达,Glyma.12G061600在种皮特异性高表达, 编码二磷酸核酮糖羧化酶;qSCC10-1有1 个基因在种皮中特异高表达,Glyma.10G204800编码亮色花青素还原酶。此外, 位于11 号和13 号染色体的3 个位点为2 个方法重复鉴定, 共包含747 个基因, 对候选区间基因进行功能注释分析及表达分析发现5 个基因可能相关,qSCC11-1区间内Glyma.11G077200(Chr11: 5,785,880..5,787,797)编码查尔酮异构酶基因, 在种皮、种脐和种子中特异性不表达,Glyma.11G027700(ANS, 白花青素双加氧酶)基因位于位点附近, 在种皮和种脐中高表达, 除此之外,Glyma.11G054600编码 bHLH 转录因子和Glyma.11G060200编码细胞色素P450, 这2 个基因在种皮中高表达;qSCC11-2区间内Glyma.11G231300(bHLH 转录因子)在种子、种脐和种皮中特异性不表达,Glyma.11G236300(bZIP 转录因子)在种皮和种脐中特异性不表达, 同时区间附近有类黄酮途径GST 基因(Glyma.11G212900)。qSCC13-1区间内Glyma.13G256900在种皮中高表达, 编码MADS 转录因子。

图5 种皮色、花色及色素合成相关候选基因表达模式[69]Fig. 5 Relative expression pattern of candidate genes related to seed coat color, flower color, and pigment synthesis[69]

表6 关联定位区间候选基因Table 6 Candidate genes within association region

3 讨论

3.1 花色关联位点和候选基因分析

本研究利用124,029 个SNP 对NAM 花色分离群体的1531 个株系进行全基因组关联分析, 使用EMMAX 方法主要定位到1 个主效位点, 位于13 号染色体(Chr13: 17,215,500..17,435,536), 包含29 个显著SNP, 其中最显著SNP 为Chr13: 17,315,536, 与已克隆的花色主效位点W1一致,W1为编码F3’5’H基因[36],Glyma.13G072100, 物理位置 17,312,471..17,317,198, 可以看到最显著位点Chr13: 17,315,536即位于该基因上。同时, 使用RTM-GWAS 分析同样位于13 号染色体Block_13_17,310,400_17,410,106与花色最显著相关, 且遗传贡献率PVE 高达65.66%,说明该位点在花色调控中起主要作用。该位点显著相关位点Chr13: 17,315,536, 包含A-G 2 个等位基因,A 等位基因型在亲本品种中为12 个白花4 个紫花, G等位基因型19 个紫花1 个白花, 与花色显著相关。同时可以发现亲本中花色可能有其他位点的调控,但是由于国内对花色的鉴定主要分为紫色和白色,本研究表型性状分类也是2 类, 因此使用EMMAX关联定位时仅定位到最主效的位点, 也是其他3 种方法定位到的最主效位点。而国外则将花色细分为9 种, 也因此鉴定出Wm、Wp、W3、W4、W2这些位点[31,35,39,44]。同时, Tassel、3vrmMLM、RTM-GWAS由于其计算模型不同分析出了更多的相关SNP, 多为1 个区间范围仅关联到1～2 个相关SNP, 可靠性需要进一步分析验证。具体分析发现3VmrMLM 方法关联的其他位点显著 SNP, Chr08: 7,175,086、Chr09: 46,432,006 、 Chr09: 48,088,955 、 Chr11:34,008,772、Chr17: 38,128,201、Chr17: 39,714,045,这几个 SNP 分别位于类黄酮途径 CHS (Chr08:8,380,000..8,520,000) 、 Myb 转录因子(Glyma.09G235100, Chr09: 45,758,761..45,760,773)、ANS(Glyma.11G027700, Chr11: 1,992,155..1,993,544)、DFR (Glyma.17G252200/W4, Chr17: 40,652,089..40,657,097)附近区间。对花色克隆相关基因的表达模式分析发现,W1基因Glyma.13G072100在花中高表达, 其他位点W3、W4、Wp、Wm候选基因在花中也均高表达, 同时发现W3[35]和Wp[42]候选基因在种皮和种脐中也高表达, 可能与DFR 和F3H 均参与类黄酮代谢途径, 影响色素沉积。

3.2 种皮色关联位点和候选基因分析

本研究利用124,029 个SNP 对9 个种皮色分离群体的561 个株系进行全基因组关联分析, 3 种以上方法重复鉴定到的位点有9 个, 其中qSCC1-2、qSCC6-1、qSCC8-1、qSCC9-1这4 个位点与经典位点G、T、I、R区间重叠。其中qSCC8-1(I)位点的相关性最高, 4 种不同方法且在9 种不同颜色属性性状中均定位到, 该区间显著性最高的连续多个SNP 位点主要分布在Chr08: 8,300,000～8,690,000 之间, 恰好与8 个串联CHS基因范围重叠, 其中最显著位点Chr08: 8,400,595 即位于Glyma.08G109500(CHS)基因下游,CHS通过转录后基因沉默调控色素沉积[73-74]。qSCC9-1(R)和qSCC6-1(T)同样被4 种不同方法关联到, 且分别与9 种和8 种颜色属性关联定位区间重叠, 说明种皮色的主要调控位点与前人研究发现的I、R、T位点一致, 主要影响种皮颜色的形成[46,71]。qSCC1-2为G位点, 与绿色种皮的形成密切相关, 在全基因组关联分析或者栽培野生杂交大豆群体中检测到[33,51,54]。同时, 基于3 种以上方法共定位到了5 个新位点, 分别位于1、5、6、10 和12 号染色体上。位于11 号和13 号染色体上的3 个新位点被2 种不同方法共定位, 其中qSCC11-1区间与Yuan 等[34]利用栽培野生大豆群体定位的种皮色相关位点qSC11一致。另外几个种皮色位点均关联到多个颜色属性, 而qSCC1-1只与a、Hue 颜色属性相关, 与表型分析发现9 个颜色属性的相关性相符合(表4)。

大量研究表明, 种皮色的相关基因主要与类黄酮生物合成途径密切相关, 在各物种中成功克隆的基因主要有3 类, 第1 类是编码类黄酮合成途径相关酶的基因, 这些基因编码查尔酮合成酶(CHS)、查尔酮异构酶(CHI)、黄烷酮3-羟化酶(F3H)、类黄酮3’5’-羟化酶(F3’5’H)、花青素合成酶(ANS)等; 第2类是编码类黄酮转运途径酶的基因, 包括编码类黄酮糖基转移酶、谷胱甘肽-S-转移酶、ABC 转运相关蛋白等; 第3 类是一些参与类黄酮合成途径的调控基因及对结构基因转录进行调控的转录因子, 包括MYB 转录因子、WKRY 结构域蛋白、WD40 重复蛋白、MYC(bHLH)类蛋白、锌指蛋白类。本研究关联的4 个经典位点, 其中3 个即编码类黄酮合成途径相关酶CHS 和F3’H 以及MYB 转录因子。新关联位点中, 分析发现qSCC10-1候选区间包含基因Glyma.10G204800编码LAR2 即亮色花青素还原酶/亮色花青素还原酶, 参与类黄酮代谢途径[55]。qSCC11-1区间内Glyma.11G077200基因属于类黄酮途径CHI 家族, 同时参与类黄酮途径的基因ANS(Glyma.11G027700, 1,992,155..1,993,544)也位于该区间附近。qSCC11-2位于参与类黄酮途径转运相关的谷胱甘肽转移酶GST (Glyma.11G212900, 30,575,282..30,577,501)附近。通过表达分析筛选候选基因, 发现11 个在种皮中高表达的基因以及4 个特异不表达的基因可能为种皮色的候选基因, 这些基因主要为类黄酮途径编码酶LAR、细胞色素P450、糖基转移酶、bHLH/MRKY/WD40/锌指蛋白转录因子(表6 和图5)。其中Glyma.10G204800、Glyma.12G060700、Glyma.06G202300、Glyma.08G062000、Glyma.11G027700不仅在种皮中表达较高, 同时在种脐中也高表达,这和Glyma.06G202300(T)、Glyma.08G062000(O)可能同时调控种脐或其他部位色素积累相关是相符合的[71]。同时类MYB 转录因子在参考基因组W82 中基本都不表达或表达量较低, 因此候选区间的转录因子是否为候选基因不能仅通过表达模式去除, 需要进一步分析验证。

3.3 大豆NAM 构建为复杂性状遗传解析奠定基础

本研究报道了一个国内大豆NAM 群体, 该群体包含35 个RIL 亚群体, 亲本遗传来源广(19 个省市), 生态型多样(5 个不同生态类型), 尤其是引入地方品种作为群体亲本, 提高了群体遗传多样性, 有利于提高定位精度和优异等位变异的挖掘。本研究利用该群体分别进行了花色和种皮色的遗传分析,均得到了较好的定位效果, 表明NAM 群体的成功构建能有效进行质量性状和数量性状的遗传分析,为进一步研究大豆复杂性状奠定了材料基础。进一步将利用NAM 群体开展产量、农艺性状以及品质性状相关表型鉴定及遗传分析, 以期加深对各性状遗传机制的理解。

4 结论

本研究利用35 份不同地区来源的代表性种质与中豆41 (公共母本)杂交, 构建了一套大豆NAM群体, 该群体遗传结构清晰。利用4 种方法对亲本间花色和种皮色具有显著差异的NAM 群体进行全基因组关联分析, 共定位到1 个主要位点qFC13-1与花色显著关联, 该位点与W1位点重合。定位到12 个位点与种皮色显著相关, 其中3 种以上方法共定位9 个位点, 2 种方法共定位3 个位点, 4 个位点与前人研究一致, 8 个可能为调控种皮色的新位点, 可以进一步研究。构建的NAM 群体适于大豆质量性状和数量性状的遗传分析, 为大豆复杂性状的遗传解析和育种实践提供了良好的基础材料。

致谢:中国农业科学院油料作物研究所周新安研究员为本项目提供中豆41 种子。