黄秋葵嫩果多糖含量性状全基因组关联分析

Genome-Wide Association Analysis of Polysaccharide Content in Young Fruits of Okra

摘要: 植物多糖是具有保健功能的生物活性成分，黄秋葵鲜嫩果具有很高的多糖含量。对黄秋葵多糖含量进行全基因组关联分析(Genome-wide Association Study, GWAS)，挖掘相关基因位点，可以为利用分子育种手段培育高多糖黄秋葵新品种奠定理论和材料基础。对180份黄秋葵种质进行GBS建库测序(Genotyping by sequencing, GBS)，结合多环境下黄秋葵多糖含量的表型数据进行全基因组范围的关联分析，以明确与多糖含量显著关联的SNP位点，挖掘位点附近可能与多糖含量相关的基因。通过4个环境共检测到46个显著关联的SNPs位点，其中Cluster1671-67622和Cluster86-210279位点至少在两个环境下稳定出现。在这2个显著关联的位点附近的100 Mb区域内发掘出黄秋葵多糖含量性状关系密切的候选基因，经过分析筛选确定UGT-糖基转移酶和MYB转录因子为候选基因。UGT-糖基转移酶和MYB转录因子与黄秋葵嫩果中多糖的合成、积累和分解紧密相关。本研究结果为阐明多糖的遗传调控机制提供了新的视角，亦为黄秋葵嫩果多糖含量的遗传改良奠定了一定的理论基础。

Abstract: Plant polysaccharide is a biological active component with health care function. The fresh and tender fruit of okra has a high content of polysaccharide. Genome-wide Association Study (GWAS) on the polysaccharide content of okra and the mining of relevant gene loci can lay a theoretical and material foundation for the cultivation of new varieties of okra with high polysaccharide by molecular breeding methods. Based on the accurate determination of polysaccharide content of 180 okra Germplasms in multi environment, the whole genome sequencing (Genotyping by sequencing, GBS) of the tested materials was carried out, and the genome-wide association analysis was carried out combined with the phenotypic data of polysaccharide content in multi environment, to identify the SNP sites significantly associated with polysaccharide content. Meanwhile, excellent genes which may be related to the content of polysaccharide near the sites were mined. A total of 46 SNPs were detected in four environments, and two of them were selected as candidate sites. In the 100 Mb region near these two significant loci, candidate genes closely related to the polysaccharide content traits of okra were found. UGT-glycosyltransferase and MYB transcription factor were finally identified as candidate genes. UGT-glycosyltransferase and MYB transcription factor are closely related to the synthesis, accumulation, and decomposition of polysaccharides in okra tender fruit. The results of this study provided a new perspective to clarify the genetic regulation mechanism of polysaccharides, and laid a theoretical basis for the genetic improvement of polysaccharide content in okra tender fruit.

1. 引言

黄秋葵(Abelmoschus esculentus L.)为锦葵科(Malvaceae)秋葵属(Abelmoschus Medicus)一年生草本植物。其用途十分广泛，嫩果、花均可供食用，种子可榨油，丰富的茎秆纤维也具有较高利用价值 [1]。黄秋葵的果荚中含有丰富的糖类物质，主要是半乳聚糖、阿拉伯聚糖、果胶类多糖、及少量的糖蛋白混合物。多糖作为天然的生物大分子，具有多种生物活性，如提高免疫力、增强身体耐力、抗疲劳、降低血糖，以及抗肿瘤等 [2]。因此，黄秋葵因其较高的营养价值和显著的保健和食疗效果越来越受到人们的关注和喜爱 [3] [4]。培育营养价值高和保健功效好的新品种是目前黄秋葵的育种目标和方向，当前可供黄秋葵育种利用的多糖功能性主效基因匮乏，严重影响了高多糖黄秋葵品种的培育的发展。目前国内外对黄秋葵多糖的提取及加工工艺报道较多，但对黄秋葵高多糖种质资源的鉴定评价与关键基因挖掘并无深入研究。

近些年，随着黄秋葵所含功能营养保健成分不断被发掘，黄秋葵作为一种保健型蔬菜已经在欧美、东南亚及非洲等地大面积种植，相关基础研究也引起了学者们的广泛关注。对黄秋葵的基础研究主要集中在遗传育种 [5] [6]、栽培生理 [7] [8]、营养成分 [9] 和药用功效 [10] [11] 等方面，而有关分子生物学的研究特别是多糖含量相关功能基因研究方面的研究甚少。王旭等(2014)利用克隆和RT-PCR技术，首次获得了黄秋葵查尔酮合成酶基因(CHS)cDNA的全长序列 [12]。张少平等(2017)通过高通量测序掌握了紫色黄秋葵叶片基因表达的情况，对紫色黄秋葵叶片的代谢途径进行了初探 [13]。李和平等(2018)利用RNA-seq技术分析黄秋葵果实基因表达情况，以及了解各代谢途径关键基因的表达模式，为黄秋葵遗传分析提供了一定的依据 [14]。对于植物多糖含量相关基因的挖掘已有一些报道，韩文韬(2019)以海带全基因组和转录组数据为基础，利用生物信息学方法，在整个生物界范围内对岩藻多糖合成途径中的重要基因家族进行比较基因组学研究，推测FucXylTs可能参与了褐藻中岩藻多糖的侧链的生物合成 [15]。范华(2017)发现在盐胁迫条件下，发菜多糖分泌量增加。对盐胁迫及正常培养条件下的发菜样品进行转录组测序分析，筛选出部分可能参与盐胁迫下发菜多糖代谢相关的差异表达基因 [16]。耿腊(2021)等以前期收集的全球119份大麦基因型为材料，对不同大麦材料籽粒β-葡聚糖含量进行GWAS分析，共鉴定到10个与糖合成、转运及分解相关的酶类基因 [17]。

但是，对黄秋葵果实多糖含量相关基因的全基因组范围挖掘至今还未见报道。本研究拟利用高通量测序技术对180份黄秋葵核心种质进行全基因组测序，对其群体结构和遗传进化进行深入分析，并结合多环境下多糖表型数据进行候选基因挖掘，并开发功能标记，获得高多糖黄秋葵新材料。为揭示黄秋葵多糖合成的分子机制提供重要线索，为分子育种培育高多糖黄秋葵新品种奠定理论和材料基础，对加快开发利用黄秋葵产品、丰富区域性菜篮子工程、提高黄秋葵产品的市场竞争力具有重要意义。

2. 材料与方法

2.1. 试验材料

通过国内考察收集和国外引种，收集保存黄秋葵180份核心种质，其中来自国内8个省(自治区)的108份，来自国外10个国家和地区的72份。分别于2019年5月和2020年5月播种于湖南省望城白箬铺基地(28˚15′N, 112˚43′E) (该环境以下简称为2019白箬铺)，江苏省南通市南通大学试验基地(31˚41′N, 120˚12′E) (该环境以下简称为2020南通)，湖南省长沙市中国农业科学院麻类研究所沅江实验基地(28˚50′N, 112˚20′E) (该环境以下简称分别为2019沅江和2020沅江)试验基地。每份材料按双行区种植，3次重复，行长10 m，行宽1 m，行距30 cm，穴播，穴距50 cm，田间管理同大田生产。取花后6 d同一结位的黄秋葵嫩果烘干测定多糖含量。

2.2. 试验方法

2.2.1. 黄秋葵嫩果多糖含量红外光谱测定

分别称取冻干的黄秋葵果实样品1 mg，加入100 mg烘干衡重的KBr，置于研钵中研磨直至极细，用压片机压片，用红外光谱仪在4000~400 cm⁻¹范围内进行红外扫描。

2.2.2. 基因组DNA提取

采集180份黄秋葵核心种质材料的幼嫩叶片于液氮中速冻，用于基因组总DNA的提取。DNA浓度稀释至50~100 ng/µL，NanoDrop 2000紫外可见光度计检测DNA浓度和质量，确保符合测序分型的要求。

2.2.3. 基因组修饰

利用限制性内切酶Hae III对黄秋葵基因组DNA进行酶切。酶切体系：黄秋葵基因组DNA 500 ng，NEB Buffer 41 µL，限制性内切酶0.12 µL，补水至50 µL，37℃酶切反应15 h。酶切结束后，用纯化试剂盒纯化，并用50 µL EB回溶。对酶切产物进行平端化修复，并且对5’末端进行磷酸化的修饰。在5’磷酸化的平端片段3’末端添加一个“A”，以便后续与5’端突出的“T”Solexa接头进行互补的连接，从而提高了接头的连接效率，同时阻止了Solexa接头的自连。反应结束后再用纯化试剂盒纯化，并于10 µL EB回溶。最后连接Solexa的接头，从而使不同的DNA片段具有了相同的末端序列，从而可以在后续够顺利的杂交到流路池上，进行扩增、测序。

2.2.4. PCR扩增与序列测序

根据酶切预测软件PAProC (http://www.paproc.de/)的分析结果，切胶回收500~580 bp的DNA片段，进行PCR扩增。扩增体系如下：纯化获得的DNA产物8 µL，Forward-引物1.5 µL，Reverse-引物1.5 µL， DNA聚合酶20 µL，补充dd水至40 µL。扩增程序如下：95℃预变性30 s；95℃变性10 s，62℃复性30 s，72℃延伸40 s，20个循环；最后再经72℃延伸6 min。反应结束后使用纯化试剂盒纯化，再于30 µL EB回溶。最后将纯化获得的PCR产物进行精细定量，然后在芯片表面进行桥式PCR，使DNA片段扩增成为单分子的DNA簇。单分子的DNA簇扩增完成后，将所得产物移入IlluminaHiSeq2000中，进行测序。

2.2.5. 高通量测序

根据DNA条形码序列，将原始的reads按照其个体分类，并剔除低质量的reads，双端均比对到基因组上的为可靠的reads，其可以用以定义SNP标签。根据序列比对和错误校正的结果，选取样本平均深度达到4以上的组别，用来定义SNP标记。统计基因组上的每100 Kb范围内SNP标记个数，得到SNP标记在scaffolds上的分布，进行基于180个样本的群体内部SNP检测，统计每100 Kb范围内SNP个数。

2.2.6. 群体结构分析

为了估计样本的群体结构，根据LD分析结果，选取相互间不连锁的SNP标记，采用STRUCTURE2.2对黄秋葵材料进行了基于数学模型的亚群划分，并获得每个材料的Q值。

2.2.7. 全基因组关联分析

利用Plink对获得的SNP进行质量控制，标准是剔除最小检出率 < 85 %和最小等位基因频率 < 5 %的SNP。使用SPSS19.0对多糖含量数据进行统计分析，利用Admixture [18] 分析群体结构，使用TASSEL3.0 [19] 的一般线性模型(GLM)和混合线性模型(MLM)进行全基因组关联分析，通过SPAGeDi 1.4软件计算样品间亲缘关系K。使用TASSEL 3.0软件生成各性状全基因组关联分析的Quantile-Quantile (QQ)图和曼哈顿图。

3. 结果与分析

3.1. 不同环境下多糖含量的描述性统计

在表型鉴定数据的基础上，结合地理来源等信息，采取逐级取样策略，构建了180份不同地理来源的黄秋葵自然群体。为解决目前育成品种遗传基础狭窄和满足多用途应用提供了材料保障。在4个不同环境条件下，黄秋葵嫩果中的多糖含量的最小值为0.9834 mg/g，最大值为3.3315 mg/g，变异系数最大为16.62 %。以上结果表明供试黄秋葵种质资源的变异非常丰富，所得结果见表1。

3.2. 群体结构和进化分析

分析180份黄秋葵资源的群体聚类和群体结构(图1和图2)，分别假设研究群体的分群数(K值)为0~10，进行聚类，并对聚类结果进行交叉验证，根据交叉验证错误率的谷值确定最优分群数，发现当K = 3时ΔK值最大，所以该群体可以分为3个亚群，表明所研究群体最可能来自于3个原始祖先。基于群体SNP，通过MEGA 5.0软件的neighbor-joining算法，对180份核心种质进行了进化分析(图3)，结果显示，180份材料可以分为3大类，第1类主要为日本引进品种及含有其血缘的种质(红色表示)；第2类主要为印度及东南亚引进品种及含有其血缘的种质(蓝色表示)；第3类主要为国内野生种及相关种质(黄色表示)。

Table 1. Variation of polysaccharide content in okra core germplasm under different environments

表1. 不同环境下黄秋葵核心种质多糖含量的变异情况

Figure 1. Population structure of 180 core germplasm of okra

图1. 180份黄秋葵核心种质的群体结构

Figure 2. Selection of K value in population structure of 180 okra core germplasm

图2. 180份黄秋葵核心种质群体结构的K值选择

Figure 3. Evolutionary analysis of germplasm resources of okra

图3. 黄秋葵种质资源的进化分析

3.3. 黄秋葵嫩果多糖含量的全基因组关联分析

3.3.1. 阈值确定

利用GEC软件计算579,654个SNPs的有效标记数(En)，结果表明，En为261,342，软件建议的P = 1.16 × 10⁻⁷ (1/En)，该值作为关联分析的显著阈值。

3.3.2. 关联分析结果

对4个环境下的嫩果多糖含量的GWAS结果进行分析，针对不同模型对假阳性的控制来看，Q模型最差，K模型和Q + K模型优于Q模型。从对假阴性的控制效果看，Q + K模型比K模型稍过严格。因此，选用表现相对好的K模型的GWAS结果用于后续的位点挖掘和基因预测。基于K模型，在P ≤ 1.16 × 10⁻⁷下，4个环境一共检测到46个SNPs与嫩果多糖含量显著关联，涉及到28个位点(图4，表2)。在2019年望城环境下，共鉴定到2个显著SNP，其曼哈顿图和QQ图见图4(A)。距离这2个SNPs位点最近的基因分别是油菜黄质脱氢酶和GTP结合蛋白(表2)。在2019年沅江环境下，共鉴定到27个显著SNP其曼哈顿图和QQ图见图4(B)，距离这27个位点最近的基因功能分别是GRAS家族基因(1个)、胞质烯醇化酶(2个)、UGT-糖基转移酶(4个)、蛋白二硫异构体酶(1个)、植物细胞器识别域(3个)、P450蛋白(1个)、转运蛋白(3个)、信号转导蛋白(1个)、核糖体蛋白(1个)、DNA结合蛋白(5个)及未知功能基因(5个) (表2)。在2020年沅江环境下，共鉴定到14个显著SNP其曼哈顿图和QQ图见图4(C)，距离这14个位点最近的基因功能分别是UGT-糖基转移酶(1个)、翻译后修饰蛋白(1个)、光系统I反应蛋白(1个)、DSBA氧化还原酶(1个)、信号转导蛋白(1个)、GDSL脂肪酶(1个)、MYB转录因子(5个)及未知功能基因(3个) (表2)。在2020年南通环境下，共鉴定到3个显著SNP其曼哈顿图和QQ图见图4(D)，距离这3个位点最近的基因功能分别是GABA转录因子(1个)，MYB转录因子(1个)及细胞分裂蛋白(1个)。

A：2019望城；B：2019沅江；C：2020沅江；D：2020南通。

Figure 4. Manhattan and Quantile-Quantile diagrams of different environments

图4. 不同环境条件下的曼哈顿图和Quantile-Quantile图

3.3.3. 候选基因确定

由于该群体的连锁不平衡衰减距离是50 Mb (r² = 0.1)，因此，在峰值SNP的上下游各50 Mb，即100 Mb范围内搜索候选基因。在每个关联显著的SNPs位点周围100 Mb的区域内寻找可能的候选基因，并且结合候选基因在多个环境下的出现情况，最终确定了2个候选基因。Cluster1671-67622位点在2019年沅江和2020年沅江均被检测到，且其最近的基因功能被注释为UGT-糖基转移酶。另外Cluster86-210279在2020年沅江和2020年南通均被检测到，该位点最近的基因被注释为MYB转录因子。因此初步确定UGT-糖基转移酶和MYB转录因子为候选基因(表2)。

Table 2. Significantly associated locus information and genes in different environments

表2. 不同环境下显著关联位点及基因功能注释信息

4. 讨论

GWAS以遍布整个基因组的单核苷酸多态性分子标记为基础，对复杂性状进行直接的关联分析，被认为是一种确定影响重要性状的分子标记的有效方法 [20] [21]。随着测序技术的飞速发展，利用这种普通分子标记进行关联分析研究呈现出诸多缺点。基于SLAF-seq的GWAS已在玉米、水稻、大豆等多种作物中全面开展 [22] [23] [24] [25]。基于SLAF-seq测序的GWAS挖掘植物优异等位基因技术已经非常成熟，效果十分显著。与水稻、大豆、玉米、棉花等作物相比，黄秋葵GWAS及优异等位基因挖掘的研究还远远不够，并且基于全基因组测序的GWAS在黄秋葵中还未见报道。本研究首次利用二代测序技术/SLAF-seq在全基因组范围内筛选与黄秋葵多糖含量性状显著关联的SNPs位点，并确定相应的候选基因。

本研究绘制了黄秋葵多糖含量性状不同环境条件下的QQ图以检验群体校正的准确性，QQ图结果显示，观测值与期望值基本相符，仅在末端存在少许偏离位点，说明群体结构校正的结果准确可靠，关联分析结果并没有因为群体分层而产生假阴性。本研究利用K模型及4个环境下的表型数据对黄秋葵的多糖含量性状进行GWAS，共筛选得到46个效应位点与黄秋葵多糖含量性状显著关联(P < 1.16E−07)。利用两个环境条件下同时关联到的显著关联的SNPs位点进行筛选，最终得到与多糖含量显著关联的SNPs位点有2个。由于这2个SNPs位点可以稳定出现，说明它们分别对黄秋葵的多糖含量性状具有重要的影响，可以作为影响黄秋葵多糖含量性状的重要候选遗传标记。对于植物多糖及单糖的全基因组关联分析报道寥寥无几，其中耿腊等(2021)利用MLM和GLM模型针对大麦籽粒β-葡聚糖含量进行全基因组关联分析，合并2个模型重叠位点后共得到44个显著位点，其中HORVU5Hr1G022710基因在2个模型、2个地点均被鉴定到，故被认为是与β-葡聚糖含量显著相关的候选基因 [17]。该研究结果并未与本研究筛选到相同的候选基因，究其原因主要是由于两者研究的物种和性状都存在差别。

在关联显著的候选SNPs标记周围100Mb区域内寻找可能的候选基因，结果找到的UGT-糖基转移酶和MYB转录因子等2个基因，初步筛选为黄秋葵多糖含量性状的候选基因。在整个多糖产生的代谢途径中，UGT-糖基转移酶参与多糖重复单元的形成，其主要作用于活性糖基，使其从糖基供体转移至糖基受体，最终形成糖苷键，因而对于多糖的生物合成具有重要的调控作用 [26] [27]，随着分子生物学的快速发展，目前借用基因工程手段，包括基因的超表达及基因敲除等，多个物种都对糖基转移酶的功能进行了深入的研究 [28]。本研究恰好筛选到UGT-糖基转移酶基因，因而本研究将糖基转移酶作为重要的候选基因之一。MYB转录因子是植物中数目庞大、功能注释最全面的转录因子之一，在很多生命及代谢过程中扮演着非常重要的角色，MYB转录因子已成为目前植物表达网络调控及基因功能研究的热点。植物MYB转录因子的主要生物学功能有参与植物的初生及次生代谢过程 [29] [30]，参与植物细胞形态和模式建成 [31] [32]，参与植物多个生长发育过程 [33] [34] [35]，参与植物逆境胁迫应答 [36] [37] [38]。而本研究发现MYB转录因子很有可能参与植物多糖的合成，具体功能有待进一步的研究验证。

5. 结论

利用SLAF测序的方法对180份黄秋葵种质资源的多糖含量性状进行GWAS，通过4个环境共检测到46个显著关联的SNPs位点，其中2个位点被筛选为候选位点。在这2个显著位点附近的100 Mb区域内发掘出黄秋葵多糖含量性状关系密切的候选基因，经过分析最终确定UGT-糖基转移酶和MYB转录因子为候选基因。以上显著关联的SNPs位点和候选基因信息可以作为黄秋葵多糖含量性状改良的理论依据。

基金项目

感谢南通市科技项目(JC2020156)；湖南省自然科学基金面上项目(2021JJ30768)；中国农业科学院麻类研究所横向课题(IBFC-KYSY-080)；国家自然科学基金项目(31801409)的资助。

参考文献

NOTES

^*通讯作者。