大豆叶型性状全基因组关联分析与候选基因鉴定

王 琼 朱宇翔,2 周密密 张 威 张红梅 陈 新 陈华涛,*崔晓艳,*

1 江苏省农业科学院经济作物研究所 / 江苏省高效园艺作物遗传改良重点实验室, 江苏南京 210014; 2 扬州大学园艺园林学院, 江苏扬州 225009

大豆是世界重要的经济作物, 是人类和家养动物蛋白质和油脂的重要来源[1]。近年来, 我国大豆对外依存度较高, 供需形势日益严峻[2]。单产水平低是我国大豆生产面临的最大困境, 培育高产大豆品种、提高大豆单产成为扭转我国大豆被动局面的首要任务[3]。作物的产量是复杂性状的集合, 受许多形态、生理和农艺性状的影响。叶片是植物进行光合作用的主要器官, 合适的叶片形态有利于塑造理想株型, 增强光合能力, 从而提高作物的产量[4]。

理想株型育种是实现大豆超高产的重要途径[5]。叶型是大豆株型性状的重要组成部分, 影响冠层结构、光合作用效率, 并最终影响产量[6]。通过对叶型结构进行优化, 可以实现培育高产大豆新品种的目标。彭玉华等[7]认为, 按叶型垂直分布不同, 大豆的株型可分为3 种类型, 即下阔叶-上阔叶、下阔叶-上窄叶和下窄叶-上窄叶, 其中以下窄叶-上窄叶类型的单荚粒数最多, 以下阔叶-上阔叶和下阔叶-上窄叶类型的百粒重最大。尽管下阔叶-上窄叶类型的单荚粒数略低于下窄叶-上窄叶类型, 但却显著高于下阔叶-上阔叶类型, 因此, 以下阔叶-上窄叶为目标进行选择, 可以在稳定其他产量因素的前提下, 使育成品种的单荚粒数显著提高, 从而实现对大豆品种的产量改良。

然而, 目前大豆叶型相关的研究大多集中在单一类型的叶形性状上, 关于叶型垂直分布的异形叶的研究还较少。Chen 等[8]利用由叶片形状明显不同的3 个大豆材料衍生构建的2 个重组自交系(recombinant inbred lines, RIL), 通过连锁分析共检测到8 个与异形叶性状相关的QTL, 再根据2 个作图群体3 个亲本的序列差异、RT-qPCR 分析和基因功能注释分析,鉴定出4 个异形叶候选基因, 其中包含大豆叶形调控基因Ln。

近年来, 随着高通量测序和高密度SNP 基因分型等技术的快速发展, 全基因组关联分析(Genomewide association study, GWAS)已广泛应用于大豆复杂性状的遗传解析[9-12]。本研究利用283 份大豆种质资源数据, 通过分析一系列叶形相关性状, 包括叶长(leaf length, LL)、叶宽(leaf width, LW)、叶形指数(leaf shape, LS)和异形叶指数(heterophylly index, HI)等, 明确异形叶对大豆产量形成的影响; 检测与目标性状关联的区域与位点, 探究大豆叶型相关性状形成的分子机制, 包括叶片的大小、形状及其层次分布等。本研究将有助于解析大豆理想株型的组成要素, 挖掘鉴定与大豆叶型性状关联的等位变异和候选基因, 以期为大豆产量改良提供指导。

1 材料与方法

1.1 试验材料与表型鉴定

供试群体为283 份大豆种质资源组成的自然群体, 包含野生种质19 份、地方品种52 份和212 份栽培品种, 这些大豆材料代表了丰富的遗传多样性。所有材料于2021 年6 月25 日和2022 年6 月20 日在江苏省农业科学院六合基地种植2 年。种植3 个重复, 随机区组设计, 大田常规管理。

测量上部和下部叶片的叶长、叶宽, 其中上部叶测量倒三叶, 下部叶测量第一复叶, 分别测量3 个单叶(图1)。用上部叶长除以上部叶宽得出上部叶形指数, 用下部叶长除以下部叶宽得到下部叶形指数, 再计算上下部叶形指数的比值, 得到异形叶指数。

图1 试验设计示意图Fig. 1 Schematic representation of the design in the experiment

1.2 重测序与基因分型

利用CTAB 法[13]提取全部供试群体材料的叶片基因组DNA, 利用Illumina HiSeq 高通量测序仪进行测序, 每个样品的平均测序数据量不低于10 倍基因组覆盖深度。使用Trimmomatic 软件[14]对测序产生的原始数据进行过滤, 去除接头污染和低质量数据后, 将得到的高质量数据比对大豆参考基因组,检测供试群体SNP 标记。具体SNP 检测流程参考本课题组前期相关工作[15]。

1.3 全基因组关联分析

筛选最小等位基因型频率(minor allele frequency,MAF)≥5%、缺失率＜10%的SNP 位点, 并利用Beagle软件[16]进行基因型填充后, 结合7 个叶型相关性状的表型数据, 运用EMMAX 软件[17]中的混合线性模型(linear mixed model, MLM)进行GWAS 分析。

1.4 候选基因鉴定

根据大豆群体连锁不平衡衰减距离, 选取显著SNP 位点及其上下游一定区段(120 kb)作为该位点候选基因分析区间, 结合基因在大豆各个组织的表达量数据、基因的注释信息和已报道相关QTL 的连锁定位, 挖掘鉴定与叶型性状关联的遗传标记或候选基因。

1.5 实时荧光定量PCR

使用北京天根生化科技有限公司RNAsimple 总RNA 提取试剂盒提取根、根瘤、茎、叶片、茎端分生组织、花、荚和种子的总RNA, 使用PrimeScript RT reagent kit (TaKaRa)合成cDNA。利用Primer Premier 5.0 软件设计引物, Forward: 5'-GCCCCAGG CTGAAGAAAATC-3', Reverse: 5'-CCACTCGAGGG AAGAAGTAGC-3', 使用Tubulin基因作为内参, 正向引物序列为: 5'-GGAGTTCACAGAGGCAGAG-3',反向引物序列为: 5'-CACTTACGCATCACATAGCA-3'。预混液使用SYBR Green Realtime PCR Master Mix (Toyobo), 采用ABI 7500 Real-Time PCR System分析数据。用2-ΔΔCt法计算基因的相对表达量, 每个样品3 次重复。

2 结果与分析

2.1 叶型相关性状表型分析

283 份材料于2021 年、2022 年在江苏南京连续种植2 年, 分别测量R6～R8 期群体材料上下部叶片的叶长、叶宽, 计算叶形和异形叶指数。对2 个环境的上下部叶长、叶宽、叶形以及异形叶指数等7个叶型相关性状的表型值进行统计分析发现, 大豆种质资源叶型性状的变异幅度较大, 单个环境的上部叶长极差范围为10.57～11.00, 上部叶宽极差范围为7.20～7.63, 上部叶形极差范围是1.53～1.89; 下部叶长的极差范围为8.10～10.30, 下部叶宽极差范围为5.23～6.70, 下部叶形极差范围是1.02～1.71, 异形叶指数的极差范围为1.08～1.21 (表1 和图2)。7 个性状的表型值都呈现正态分布, 表明叶型相关性状为多基因控制的数量性状(图2)。2 个环境的性状间存在显著相关性(P≤0.001), 表明本研究考察的表型数据具有较高的可靠性和重复性(表2)。

表1 大豆群体叶型相关性状的统计分析Table 1 Statistical analysis of leaf characteristics traits of soybean association panel

表2 叶型相关性状的相关系数Table 2 Correlation coefficients of leaf characteristics traits

图2 大豆关联群体叶型相关性状的分布Fig. 2 Distribution of leaf characteristics traits of soybean association panel

2.2 不同叶型性状相关性分析

对不同叶型性状进行相关性分析(图3)发现, 上下部叶片的叶长和叶宽之间都呈显著正相关, 而叶宽与叶形性状存在极显著负相关; 下部叶长与上部叶长之间、下部叶宽与上部叶宽之间均呈显著正相关, 同时上下部叶形性状间也呈显著正相关, 说明大豆种质资源叶形垂直分布以下阔叶—上阔叶和下窄叶—上窄叶类型为主; 此外, 下部叶宽与上部叶形性状之间、上部叶宽与下部叶形性状之间均呈显著负相关, 同样说明大豆叶形分布以上下都为阔叶或窄叶为主; 值得注意的是, 下部叶长与下部叶形性状之间, 存在极显著正相关, 而上部叶长与上部叶形性状无显著相关性; 异形叶指数与上部叶形性状存在显著正相关, 而与上部叶宽之间、与下部叶长和下部叶形性状之间均呈显著负相关。

图3 大豆关联群体叶型性状间的相关系数Fig. 3 Correlation coefficients among leaf characteristics traits in soybean

2.3 全基因组关联分析

利用3,319,306 个高质量SNP (MAF≥0.05)进行基因分型。结合叶型性状表型数据, 使用EMMAX软件的混合模型进行全基因组关联分析, 合并属于同一关联信号的SNP, 共检测到181 个显著相关SNP 位点, 其中能够在2 个环境或多个性状中重复检测到的GWAS 信号有18 个, 分布在5、10、13、16、17、19 和20 号染色体上(表3 和附图1)。位于19 号染色体的上部叶宽关联信号19_23.9 在2021 和2022 两个环境中被重复检测到; 位于5 号染色体的关联信号在2021 年下部叶长和下部叶宽2 个性状中被重复检测到; 位于10 号染色体的关联信号能够在2021 年的上部叶形和异形叶指数性状中被重复检测;位于13 号染色体上的GWAS 信号在2022 年上部叶宽和2021 年上部叶形性状中被重复检测到; 16 号染色体上的关联信号在2021 年的上下部叶长性状中被重复检测; 位于17 号染色体的关联信号能够在2022 年上下部叶长性状中被重复检测到; 另外, 位于20 号染色体上的GWAS 信号在2022 年上下部的叶形性状中被重复检测; 值得注意的是, 位于19 号染色体上的关联信号19_45.1 在2021 年上部叶形和异形叶指数、2022 年上部叶宽和上部叶形等4 个环境或性状中被重复检测到。表明大豆叶片的发育过程, 受到复杂的遗传网络调控。

表3 叶型相关性状的关联信号Table 3 GWAS signals of leaf characteristics traits

在所有叶形性状相关的SNP 位点中, 最为显著的位点为位于2 号染色体上的上部叶宽相关位点(P=6.15E-9)。叶长性状最显著相关的 SNP 位点(Chr20:37941224)位于20 号染色体(P=1.49E-7)。对于异形叶指数性状, 最为显著相关的SNP 位点位于16 号(P=2.58E-8)和19 号(P=2.84E-8)染色体, 其中位于19 号染色体上的相关SNP (Chr19:45155943)同样也是叶形指数性状(2021 上部叶形指数性状)最显著相关的位点。综上所述, 这些显著位点与特定叶型性状关联, 特别是在不同环境中多次检测到的相关位点, 将有助于阐明大豆中叶型调控和异形叶形成的机制, 为通过叶型设计培育高产大豆品种提供优异的基因资源。

2.4 关键候选基因鉴定

利用检测到的与叶型性状相关的位点, 通过基因的表达谱数据、拟南芥中同源基因的功能注释,确定与大豆叶片发育和异形叶形成相关的候选基因。基于大豆基因组注释信息, 在19 号染色体Chr19:45151266 位点下游11.7 kb, 鉴定到上部叶宽、上部叶形和异形叶指数性状的候选基因Glyma.19G 194100。基因功能注释表明,Glyma.19G194100是LITTLE ZIPPER 3 (ZPR3)基因的同源基因, 并且该基因在幼叶中具有较高的表达量。在拟南芥中,ZPR3是茎端分生组织(Shoot apical meristem, SAM)发育过程中HD-ZIPIII 转录因子的竞争性抑制基因,该基因通过形成非功能性异二聚体起作用, 参与茎端分生组织的发育和侧生器官发生, 对茎端分生组织中干细胞发挥正常功能至关重要[18]。叶片是茎端分生组织产生的第一类侧生器官, 在植物的发育中具有重要地位。因此Glyma.19G194100可能是调控大豆叶片形状和垂直分布的关键候选基因。另外, 在19号染色体的相关位点Chr19:45155943 上游109.9 kb发现基因Glyma.19G192700, 该基因是 Growthregulating factor 4 (GRF4)基因的同源基因。生长调节因子GRF 参与调节植物生长发育的各个方面, 包括叶片大小、叶片衰老、应激反应、雌蕊发育和种子大小等[19-24]。GRF4 作为转录调节因子GRF 家族的成员, 能够调节植物的株型和种子大小、氮肥利用效率等[25-26]。拟南芥叶片在近端和远端轴上具有特定的生长极性, miR396 在幼嫩发育叶片的远端优先积累并靶向GRF 从而调控这种特定的叶片生长发育模式[27]。利用荧光定量 PCR 分析Glyma.19G 192700在大豆各组织中的表达模式(图4)发现, 该基因在茎端分生组织中高表达。SNP 注释分析表明,Glyma.19G192700第3 个外显子的第168 个碱基发生C→A 替换, 导致编码的氨基酸从亮氨酸变为异亮氨酸。说明,Glyma.19G192700可能是大豆控制叶型等性状的关键候选基因。

图4 候选基因Glyma.19G192700 在大豆各组织中的表达模式Fig. 4 Relative expression pattern of candidate genes Glyma.19G192700

3 讨论

叶型包括叶片的形状和垂直分布不仅能够影响植物对光能的吸收和利用、光合效率, 也可以改变群体的冠层结构和冠层空间光分布, 从而影响通风透气性和光能利用效率[28-29], 并最终影响产量。本研究通过对283 份大豆种质资源的7 个叶型相关性状进行连续2 年的考察, 并通过全基因组关联分析检测到181 个相关位点。其中能够在2 个环境或多个性状中重复检测到的位点有18 个, 能被重复检测到的位点较少, 说明大豆叶型性状可能受环境影响较大。

Tischner 等[30]通过对3 个大豆品种相互杂交得到的3 个RIL 群体进行QTL 分析表明, 大豆的种子败育、叶形和繁殖等相关性状由多个显著数量性状位点QTL 控制, 占遗传变异的50%左右。这些QTL分别与影响雌雄不育的Ms1、Ms6和St5基因以及影响小叶数目和叶片形状的Lf1和ln基因相关。其中,以窄叶命名的ln基因具有多效性, 长期以来ln一直被认为还与大豆的单荚粒数性状紧密相关[31]。窄叶由隐性基因ln控制, 其显性等位基因Ln控制阔叶的形成。Ln是拟南芥JAGGED(JAG)基因的同源基因,大豆叶片由阔叶(Ln)到窄叶(ln)的转变与JAG基因中EAR 基序的氨基酸替换有关[32-33]。本研究在20 号染色体Chr20:36152820 位点上游322.7 kb 发现Ln基因(Glyma.20G116200), 相关性状为2022 年下部叶形指数(P=7.84E-7)。值得注意的是, 经典的开花抑制因子E1也能够通过直接抑制叶片发育相关基因TCP14和TCP29参与大豆叶片发育的过程[34]。在本研究中,E1基因(Glyma.06G207800)位于6 号染色体上的相关位点Chr06:20059560 (2021 上部叶宽,P=8.37E-6)的下游, 距离为147.5 kb。Chen 等[8]利用具有异形叶性状的冀豆17 分别与卵圆形叶的冀豆12和披针形叶的绥农14 构建的2 个RIL 群体, 鉴定到异形叶性状相关的候选基因4 个, 分别位于1、18、19 和20 号染色体上。其中19 号染色体上的候选基因Glyma.19G194300位于本研究中鉴定的相关位点Chr19:45151266 附近。生长调节因子GRF 具有保守的QLQ 和WRC 结构域, 并且它们中的大多数可以与GRF 互作因子(GRF-Interacting-Factor, GIF)相互作用。此外, 大多数GRF 受miR396 的调控, miR396、GRF 及其相互作用因子GIF 共同组成复杂的调控网络。在水稻中, OsmiR396c-OsGRF4-OsGIF1 分子调控模块调节水稻的籽粒大小和产量[25]。本研究中预测的候选基因Glyma.19G192700, 作为拟南芥GRF4基因的同源基因, 可能具有一因多效的作用, 除与异形叶等叶型性状相关外, 可能也与粒重等性状相关。

4 结论

本研究初步分析了大豆异形叶等叶型性状的遗传特征, 利用全基因组关联分析检测到181 个叶型性状相关SNP 位点, 其中能够在2 个环境或多个性状中重复检测到的位点 18 个。最显著相关位点Chr19:45155943 位于19 号染色体, 在其附近鉴定到候选基因Glyma.19G192700、Glyma.19G194100, 分别被注释为 Growth-regulating factor 4 (GRF4)和LITTLE ZIPPER 3 (ZPR3)基因的同源基因。此外, 在其他多个位点附近发现大豆已知叶形和叶片发育相关基因Ln和E1。