青贮玉米品质相关性状的全基因组关联分析

关键词 青贮玉米； 品质性状； 全基因组关联分析 ； 单核苷酸多态性位点； 优异等位变异； 候选基因

近年来，随着人们生活水平不断提高，肉类、蛋类、奶类等畜牧产品的市场需求也逐渐扩大，畜牧业也因此迅速发展［1］。畜牧业的发展必然会加大对饲料的消耗，国内饲料原料（玉米、小麦、大麦、高粱等）已经不能满足需求，仍需从国外进口一部分玉米等饲料原料［2］。

青贮玉米作为一种优质饲料资源，生物产量高、营养价值高、适口性好、消化吸收率高，其资源创制和种植推广可有效保障畜牧业的发展［3-4］。2023 年中央一号文件明确指出“大力发展青贮饲料”，据统计，2023 年全国青贮玉种植面积达333 万hm2 以上（http：//www.nahs.org.cn/）。

品质性状是评价青贮玉米的指标之一，如粗蛋白含量、淀粉含量、中性洗涤纤维（neutral detergentfiber， NDF）含量、酸性洗涤纤维（acid detergent fi⁃ber，ADF）含量和可溶性糖含量等［5］。其中，粗蛋白作为生命活动必需的基础养分，其含量越高，饲料营养品质越好［6］，越有利于牲畜的生长与发育。淀粉在青贮玉米中含量较高，是饲用价值的重要基础，更是牲畜所需能量的主要来源［7］。中性洗涤纤维（NDF）和酸性洗涤纤维（ADF）是衡量青贮玉米作为粗饲料好坏的重要指标［8］。NDF 主要包括木质素、纤维素和半纤维素，可被牲畜部分吸收利用，其可吸收利用程度与牲畜的采食量有关［9］；ADF 主要包括饲料中的木质素和纤维素，牲畜一般不能吸收利用，且ADF含量越高，消化率越低［10］。一定量的可溶性糖有利于青贮发酵过程的进行，其含量的高低是能否青贮成功的关键因素［11］。同时，体外干物质消化率是衡量青贮玉米营养物质可消化性的尺度，是评价青贮玉米营养价值的综合指标之一［12］。因此，揭示粗蛋白、淀粉、NDF、ADF、可溶性糖和体外干物质消化率的遗传机制对青贮玉米的选育具有重要意义，如利用分子标记辅助选择，通过遗传改良提高青贮玉米的营养价值。

前人对玉米粗蛋白、淀粉、NDF、ADF、可溶性糖和体外干物质消化率的遗传规律进行了一些研究。Cardinal 等［13］以源自B73 × B52 的200 个重组自交系为试验材料，在叶鞘中检测到11 个QTL 与ADF显著相关，在茎秆中检测到11 个与ADF 和NDF 均显著相关的QTL。Park 等［14］以158 个02S6140 ×KSS22 的F2 单株为定位群体，结合SSR 标记，鉴定出4 个与糖含量相关的QTL。Zhang 等［15］以298 个BHO × B73 的F2：3家系为作图群体，定位到与油分、蛋白质和淀粉含量相关的QTL 分别有6 个、6 个和5个。Wang 等［16］以368 份自交系为材料，鉴定出73 个与ADF 含量显著相关的单核苷酸多态性位点（sin⁃gle nucleotide polymorphisms，SNPs）。李灿［17］利用全基因组关联分析的方法挖掘出了2 个与ADF 关联的候选基因。李雪莹等［18］通过构建F2分离群体，利用SSR 标记对玉米总淀粉、直链淀粉和支链淀粉基因QTL 进行初步定位，共检测到17 个与玉米淀粉含量有关的QTL。卢雨晴等［19］利用全基因组关联分析对玉米自交系及测交群体的ADF 性状进行显著SNPs 筛选，定位到部分与ADF 显著相关的位点。Truntzler 等［20］利用QTL 元分析手段，对59 个与消化率性状相关的QTL 进行分析与比较，发现26 个与消化率相关的元QTL 位点。

虽然对青贮玉米相关营养指标性状的遗传规律进行了相关研究，也挖掘到了部分显著相关的遗传位点，但主要集中在NDF 和ADF 等性状，且未有高密度的遗传图谱，有待进一步深入研究。此外，大多数研究是分别对这些指标的相关遗传位点进行单独挖掘，同时对多个营养指标进行定位的相关研究较少；对多效遗传位点挖掘的缺乏不利于一因多效位点在分子辅助育种上的应用。同时，前人所用材料主要以温带材料为主，对Suwan 等热带材料青贮品质相关遗传机制研究较少。因此，本研究基于包含5.6 万个SNPs 的Maize SNP 50 芯片对183 份含有Suwan 和Tuxpeño 等热带资源的自交系进行基因型鉴定，结合表型数据，进行全基因组关联分析，挖掘出与玉米粗蛋白、淀粉、NDF、ADF、可溶性糖和体外干物质消化率相关的QTNs，并利用生物信息学筛选出相关的候选基因，以期为后续选育营养品质优良的青贮玉米品系提供科学依据。

1 材料与方法

1.1 试验材料

本研究的供试材料由183 份骨干自交系组成，包括Suwan、Tuxpeño、Iodent 和Lancaster 等资源，遗传背景丰富，且不同自交系间的株型、果穗等性状差异较大，如青贮品种黔青446 的亲本QB446 为热带材料，其植株叶片宽大且叶间密，生物产量大；温带自交系齐319、HCL645 等植株叶片数较少且窄，生物产量较低。

1.2 田间设计

183 份供试材料于2019 年4 月种植于贵州省农业科学院试验基地（106°39′ E， 26°30′ N），行长3 m，行距60 cm，株距25 cm，单行区，2 次重复，自然散粉，常规田间管理。

1.3 基因型鉴定

当供试材料长到5 叶期，利用离心管采集叶片，采用改良CTAB 法提取基因组DNA，检测合格后送北京康普森生物技术公司，基于包含5.6 万个SNPs标记的高密度基因芯片对供试材料进行基因型鉴定。

1.4 表型数据的测定与分析

在籽粒乳线位置1/2 时，每个材料选择5 株长势、株型等较为一致的植株进行收获，并从地上部10cm 处收割；接着整株进行粉碎并混合均匀，每行材料各取2 份1 kg 样品分别装入布袋，并于105 ℃烘箱中进行1～2 h 杀青，再在65 ℃下烘干至恒质量，然后进一步粉碎为粉末，过孔径0.25 mm 筛备用。

本研究主要测定玉米粗蛋白、淀粉、NDF、ADF、可溶性糖和体外干物质消化率；其中，粗蛋白的测定参照GB/T 6432-2018《饲料中粗蛋白的测定凯氏定氮法》的操作程序执行；淀粉的测定参照GB/T20194-2018《动物饲料中淀粉含量的测定 旋光法》的操作程序执行；NDF 的测定参照GB/T 20806-2006《饲料中中性洗涤纤维（NDF）的测定》的操作程序执行；ADF 的测定参照NY/T 145-2007《饲料中酸性洗涤纤维的测定》的操作程序执行；可溶性糖的含量采用蒽酮-硫酸法测定；体外干物质消化率采用胃蛋白酶-胰蛋白酶两步体外消化法进行测定。

参照GB/T 25882-2010《青贮玉米品质分级》对供试材料进行品质分级（表1）。利用Excel 2019软件对表型数据进行整理，用SPSS 统计分析软件进行描述性统计分析。

1.5 全基因组关联分析以及候选基因的预测

根据最小等位基因频率（minor allele frequency，MAF） gt; 0.05 且样本缺失率lt; 20% 的标准，筛选出高质量SNPs。基于R 4.2.1 软件的GAPIT 语言包，利用混合线性模型（mixed linear model，MLM）对表型和基因型进行全基因组关联分析。以Plt;0.001 为阈值，鉴定出控制玉米粗蛋白、淀粉、NDF、ADF、可溶性糖和体外干物质消化率的关键QTNs（quantita⁃tive trait nucleotides）；同时，参考Shim 等［21］的方法计算出显著SNPs 的表型变异解释率（phenotypic varia⁃tion explained， PVE）。在此基础上，结合玉米B73 参考基因组（Ref Gen V4），利用生物信息分析手段，筛选出相关的候选基因。

2 结果与分析

2.1 表型分析

测定结果显示，供试材料的粗蛋白含量6.46%～13.35%，变异系数为13.08%；淀粉含量6.71%～18.02%，变异系数为20.77%；NDF 含量40.83%～65.26%，变异系数为5.38%；ADF 含量19.46%～37.58%，变异系数为12.31%；可溶性糖含量4.31%～35.84%，变异系数为44.97%；体外干物质消化率37.05%～57.73%，变异系数为7.59%（表2）。其中，粗蛋白含量≥7%的有181份自交系；淀粉含量≥15% 的有81 份；ADF 含量≤29% 的有114 份，≤26% 的有51 份，≤23% 的有11 份；NDF 含量≤55%的有48 份，≤50% 的有7 份，≤45% 的有1 份。按照青贮玉米品质分级标准（表1），有17 份自交系的品质符合三级标准，分别为14D1310、15D551、16XH-51、17HG825、17HG837、17HG850、QB44、QB193、QB196、QB506、QB662、QB2663、QB2721、QB2778、QB2894、M03 和Mei22。其中QB506 作为金玉818的亲本，已于2019 年通过贵州省青贮玉米审定，符合该试验的结果。

为进一步了解不同性状间的相关性，对玉米粗蛋白、淀粉、NDF、ADF、可溶性糖和体外干物质消化率等性状开展了Pearson 相关性分析（图1）；结果表明，淀粉含量与粗蛋白含量、ADF 和可溶性糖均呈极显著负相关（Plt;0.01）；相关系数分别为-0.27、-0.21 和-0.42；可溶性糖与ADF 呈显著相关（Plt;0.05），相关系数为0.18；而NDF 与可溶性糖呈极显著负相关（Plt;0.01），相关系数为-0.30。

2.2 亲缘关系和主成分分析

根据主成分分析结果及其系谱信息，可将183 份自交系大致分为4 个类群（图2A）。其中，Ⅰ群主要包括Suwan 种质，如T32；Ⅱ群主要包括Iodent 种质，如HCL645；Ⅲ群主要包括Lancaster 种质，如Mo17；Ⅳ群主要包括Reid 种质，如Shen5003。基于Kinship分析（图2B），发现大部分自交系间的亲缘关系系数为0.0～0.5，仅有少部分为0.5～1.0，表明该183 份自交系间的亲缘关系较远。同时，主成分和亲缘关系作为协变量用于后续进行全基因组关联分析。

2.3 全基因组关联分析

根据最小等位基因频率MAFgt;0.05 且样本缺失率lt;20% 的标准，筛选出40 498 个高质量SNPs 用于全基因组关联分析。利用MLM 模型分别鉴定出31、61、11、36、20 和42 个与粗蛋白、淀粉、NDF、ADF、可溶性糖及体外干物质消化率显著相关的SNPs（Plt;0.001）（图3）。

与粗蛋白显著相关的SNPs 分布在9 条染色体上（图3A），其中4 号染色体上有11 个，2 号染色体有5个，3 号染色体上有4 个，1 号和10 号染色体上各有3个，5 号染色体上有2 个，6 号、7 号、8 号染色体上各有1 个。粗蛋白显著相关SNPs 的表型变异解释率为5.80%～11.40%，其中位于3 号染色体上的SYN1589 表型变异解释率最高。

与淀粉显著相关的SNPs 分布在9 条染色体上（图3B），其中6 号染色体上有19 个，2 号染色体上有16 个，4 号染色体上有7 个，9 号染色体上有6 个，5 号和8 号染色体上有4 个，1 号和10 号染色体上各有2个，3 号染色体上有1 个。淀粉显著相关SNPs 的表型变异解释率为5.78%～11.38%，其中位于6 号染色体上的SYN22472 表型变异解释率最高。

与NDF 显著相关的SNPs 分布在6 条染色体上（图3C），其中5 号染色体上有4 个，1 号和2 号染色体上各有2 个，3 号、7 号和10 号染色体上各有1 个。NDF 显著相关SNPs 的表型变异解释率为5.78%～7.85%，其中位于1 号染色体上的PZE-101058310 表型变异解释率最高。

与ADF 显著相关的SNPs 分布在8 条染色体上（图3D），其中2 号染色体上有9 个，4 号和6 号染色体上各有6 个，3 号和8 号染色体上各有5 个，1 号染色体上有3 个，7 号和10 号染色体上各有1 个。ADF 显著相关SNPs 的表型变异解释率为5.81%～10.37%，其中位于2 号染色体上的PUT-163a-29947492-1846表型变异解释率最高。

与可溶性糖显著相关的SNPs 分布在6 条染色体上（图3E），其中4 号染色体上有8 个，5 号染色体上有4 个，9 号染色体上有3 个，2 号和10 号染色体上各有2 个，7 号染色体上有1 个。可溶性糖显著相关SNPs的表型变异解释率为5.78%～7.35%，其中位于10 号染色体上的SYN19780 表型变异解释率最高。

与体外干物质消化率显著相关的SNPs 分布在7条染色体上（图3F），其中6 号染色体上有28 个，3 号染色体上有6 个，2 号、7 号和8 号染色体上各有2 个，4 号和5 号染色体上各有1 个。体外干物质消化率显著相关SNPs 的表型变异解释率为5.79%～11.33%，其中位于6 号染色体上的PZE-106069016 表型变异解释率最高。

进一步分析发现，有4 个遗传位点存在一因多效的特性。其中，6 号染色体上的SYN6712、PHM1190.3 和SYN7541 同时与淀粉、体外干物质消化率显著关联；SYN6712 对淀粉和体外干物质消化率的表型变异解释率分别为6.66%、5.80%；PHM1190.3 对淀粉和体外干物质消化率的表型变异解释率分别为10.72%、6.45%；SYN7541 对淀粉和体外干物质消化率的表型变异解释率分别为6.12%、7.27%。4 号染色体上的PZE-104072386 同时与ADF、可溶性糖显著关联，表型变异解释率分别为6.41%、6.18%。

2.4 多效位点的等位变异鉴定和表型效应分析

SYN6712、PHM1190.3、SYN7541 和PZE-104072386 是本研究发现的多效位点。在本研究供试群体中有164 份材料在SYN6712 携带T/T 基因型，10 份材料携带C/C 基因型，4 份材料携带T/C 基因型，其余5 份材料基因型缺失；并且在该位点携带T/T 基因型材料比携带C/C 基因型材料的淀粉含量和体外干物质消化率分别平均增加3.80% 和3.62%，差异均达到极显著水平（Plt;0.001）（图4A），说明T/T 基因型是SYN6712 的优异等位变异。在PHM1190.3 上有164 材料携带T/T 基因型，15 份材料携带C/C 基因型，2 份材料携带T/C 基因型，其余2 份材料基因型缺失；基因型为T/T 材料比基因型为C/C 材料的淀粉含量和体外干物质消化率分别平均增加4.12% 和2.76%，差异分别达到极显著（Plt;0.001）和显著水平（Plt;0.01）（图4B），说明T/T 基因型是PHM1190.3 的优异等位变异。在SYN7541 上有131 份材料携带A/A 基因型，44 份材料携带G/G基因型，其余8 份材料携带A/G 基因型；基因型为A/A 材料比基因型为G/G 材料的体外干物质消化率平均增加2.74%，差异达到极显著（Plt;0.001），而该位点变异对淀粉含量没有显著影响（图4C）。在PZE-104072386 上有121 份材料携带A/A 基因型，49 份材料携带C/C 基因型，9 份材料携带A/C 基因型，其余4 份材料基因型缺失；基因型为A/A 材料比基因型为C/C 材料的ADF 和可溶性糖含量分别平均降低2.03% 和4.39%，差异均达到极显著（Plt;0.001）（图4D）。其中，品质符合三级标准的17 份自交系在多效位点的优异等位变异均有2 个及以上（表3）。

2.5 候选基因预测

针对4 个多效遗传位点进行候选基因挖掘，结合公共数据库（www.maizegdb.org），对候选位点进行功能注释，寻找位点附近的候选基因，通过评估预测基因与表型性状之间的关联性，筛选到4 个相关联的候选基因。其中，同时与淀粉、体外干物质消化率显著相关的SYN6712、PHM1190.3 和SYN7541 分别与Zm00001d037272、Zm00001d037386 和Zm00001d037532 相关联。同时与ADF、可溶性糖显著相关的PZE-104072386 与Zm00001d051166 相关联。其中，Zm00001d037272 编码α-甘露糖苷酶，该酶是真核生物蛋白质N-聚糖修饰的关键酶；在植物细胞壁中的糖蛋白含有大量N-聚糖结构，α-甘露糖苷酶对N-聚糖的修饰会直接影响细胞壁的结构；因此，α-甘露糖苷酶的表达可能破坏了青贮玉米的细胞壁，使其软化利于发酵，提高体外干物质消化率；同时，参与糖蛋白的降解，与淀粉的形成密切相关。Zm00001d051166编码苯丙氨酸解氨酶，该解氨酶是植物代谢途径的关键酶之一，参与可溶性糖等次生代谢物和ADF 中木质素的合成。

3 讨论

3.1 青贮玉米亲本的重要性

我国青贮玉米产业起点晚、资源少、规模小，有待进一步发展。现阶段，我国青贮玉米育种存在的主要问题之一是缺乏优异的青贮亲本，选育符合育种目标的亲本为当前的难题之一。本研究发现17 份自交系的品质符合青贮玉米品质三级标准，且17 份自交系在多效位点的优异等位变异均有2 个及以上。其中，QB506 的粗蛋白、淀粉、NDF 和ADF 含量分别为10.98%、15.27%、54.81% 和28.43%。笔者所在团队利用QB506 作为亲本之一，组配的金玉818 已于2019 年通过贵州省青贮玉米审定。同时，金玉818的另一个亲本T32，其粗蛋白含量为10.35%、淀粉含量为16.28%。此外，我们选育的另一个青贮玉米黔青446 的双亲（QB446、QB1545）同样有较高的粗蛋白和淀粉含量，QB446 的粗蛋白和淀粉含量分别为10.75% 和14.72%，QB1545 的的粗蛋白和淀粉含量分别为10.96% 和12.97%。

热带、亚热带玉米种质具有生物产量高、品质优、抗逆性强、持绿期长等特点，符合当前青贮玉米的育种目标，对我国的青贮玉米品种选育起着关键性作用。我国育成的大部分青贮玉米品种含有一定的热带血缘［22-24］，例如京科青贮301、北青贮410、中农大青贮67、雅玉青贮8 号、雅玉青贮26 和雅玉青贮04889 等［25］。优良自交系是品种的保障，为了进一步利用热带种质，国内科研单位相继选育出了一批粮饲兼用自交系，如Suwan 种质S37、YA3729、YA3737，巴西热带种质YA8201、YA8702，Tuxpeño种质YH-1 等，这些自交系的成功选育为我国青贮玉米的发展拓宽了新的途径。常海滨［26］认为热带、亚热带玉米与温带材料组配是青贮玉米的杂种模式之一。同样，采用“热×温”模式组配的金玉818 和黔青446 为品质优、产量高、抗性好的强优势组合，适宜在贵州等西南地区种植［27-28］。

3.2 多效遗传位点分析

遗传位点的多效性对青贮玉米相关品质性状研究起着重要的作用，多效遗传位点的鉴定有助于更好地理解不同品质性状间的遗传相关性。Park 等［14］利用02S6140 和KSS22 构建作图群体，发现4 号染色体的umc1088-bnlg1265 区域同时与玉米的淀粉、葡萄糖和蔗糖有关。Wassom 等［29］发现在8 号染色体上的bnlg162–umc1130 区域同时与玉米籽粒的油分、蛋白和淀粉相关。Wang 等［16］利用368 份玉米自交系的茎秆品质性状进行全基因组关联分析，发现在1 号染色体的chr1.S_18202551、2 号染色体的chr2.S_198724413 和5 号染色体的chr5.S_79867527 均同时与ADF 和体外干物质消化率相关，4 号染色体的PZE-104075114、7 号染色体的chr7.S_2739958 和8号染色体的chr8. S_150908620 均同时与NDF 和ADF 相关，6 号染色体的PZE-106111799 同时与NDF 和体外干物质消化率相关，6 号染色体的chr6.S_155653406 同时与NDF、ADF 和体外干物质消化率相关；其中PZE-104075114 与本研究发现的多效位点PZE-104072386 相隔4.15 Mb，chr6. S_155653406 与本研究发现的多效位点SYN6712、PHM1190.3 和SYN7541 分别相距40.52、35.42 和30.84 Mb，可见，6 号染色体可能存在着大量与青贮玉米品质性状相关的多效遗传位点，有待进一步深挖，进而构建高密度青贮玉米品质性状相关的遗传图谱，利用分子辅助选择产量高、品质优的青贮玉米品系。

挖掘多效遗传位点的优异等位变异是利用分子标记辅助选择育种的首要条件之一［30-31］；因此，α-甘露糖苷酶的表达可能破坏了玉米植株的细胞壁，使其软化利于畜牧消化，提高体外干物质消化率；同时，参与糖蛋白的降解，可能与淀粉的形成密切相关［32］。Zm00001d037386 编码跨膜蛋白，跨膜蛋白是细胞与外界联系的桥梁，参与细胞膜内外物质交换、能量和信号的传递［33］，可能与淀粉等物质的吸收有关。Zm00001d037532 编码6-磷酸果糖-2-激酶/果糖-26-二磷酸酶，参与糖的代谢［34］，可促进营养物质的降解，利于牲畜消化，提高体外干物质消化率。Zm00001d051166 编码苯丙氨酸解氨酶，该酶作为植物代谢途径的关键酶之一，参与次生代谢物的合成，如可溶性糖和ADF 中木质素等［35］。这些候选基因的发现，可为品质相关性状的功能基因克隆提供一定的科学依据。