大豆7S与11S球蛋白QTL元分析及候选基因挖掘

陈庆山，张泽鑫，刘函西，齐照明，韩思宁，董晓慧

（1.东北农业大学农学院，哈尔滨 150030；2.西北农林科技大学农学院，陕西 杨凌 712100；3.中化现代农业有限公司，上海 200120；4.黑龙江省农业科学院绥化分院，黑龙江 绥化 152000）

大豆（Glycine maxL. Merr.）富含丰富植物蛋白，是人类摄入植物蛋白主要来源。豆粕可为家禽和牲畜提供丰富的可消化氨基酸，在畜牧业中应用广泛[1]。此外，因其在饮食中的营养价值和健康益处，大豆在食品生产中的需求量逐年增加[1]，大豆蛋白成分研究日益受关注。大豆蛋白根据其功能可分为贮藏蛋白、结构蛋白和防御蛋白三类，其中贮藏蛋白以球蛋白为主，约占大豆蛋白70%[2]。根据蛋白沉降系数可将大豆贮藏蛋白分为2S、7S、11S 和15S 组分，其中7S 和11S 组分占大豆种子贮藏蛋白总量比例超过70%[3]，因此7S 和11S组分研究具有重要意义。

11S球蛋白是一种六聚体蛋白，在六个亚基中每一个亚基均由酸性（A1a、A1b、A2、A3、A4 和A5）和碱性（B1a、B1b、B2、B3 和B4）多肽通过二硫键连接组成[4]。7S 球蛋白是由α′、α和β亚基组成的三聚体蛋白，7S 球蛋白主要包括β-伴大豆球蛋白、少量γ-伴大豆球蛋白以及碱性7S 球蛋白，其中β-伴大豆球蛋白占绝大部分[5]。由于11S和7S球蛋白具有不同结构和特性，其在大豆蛋白质量和特性方面发挥不同作用。2014年，Nishinari等研究报道称11S 球蛋白具有由二硫键连接的紧密结构，其乳化和起泡能力弱于7S 球蛋白[6]。2011 年，Stanojevic 等证实11S 与7S 球蛋白比率与豆腐凝胶硬度呈显著正相关[7]。刁桂珠等以4 份7S、11S 球蛋白不同亚基组成大豆种质为试材，通过测定所制豆腐品质，发现7S 与11S 球蛋白不同亚基缺失对豆腐品质产生不同影响[8]。冯芳等将姜黄素与大豆7S、11S蛋白结合，发现两者相互作用可改善姜黄素水溶性底、稳定性差的缺陷[9]。

元分析（Meta analysis）可对多项独立研究结果作综合分析，通常用于整合不同来源研究结果，此方法可克服个别研究局限性并获得精准结论[10]。目前，元分析应用于玉米、小麦、水稻、大豆以及其他作物各种性状相关QTL 研究。Guo等利用元分析对已发表的17 个大豆包囊线虫抗性相关结果进行整合与评估，缩短并定位与包囊线虫抗性相关QTL[11]；Wu 等共搜集98 个已报道的大豆生育期相关QTL，通过元分析提取“真实QTL”，为大豆生育期QTL 精细定位和基因克隆奠定基础[12]；汪霞等采用多QTL 联合分析方法，整合Soybase 数据库中信息完全的90 个大豆株高QTL，通过元分析得到大豆株高性状相关QTL[13]。国内外已开展大量有关大豆蛋白含量相关QTL 的研究，但对于7S 与11S 球蛋白相关QTL 整合分析的研究较少。

本研究利用Soybase 网站搜集42 个大豆7S 与11S球蛋白相关QTL，利用元分析方法对其进行区间整合及验证，在得到的“真实”QTL内筛选得到调控7S 与11S 球蛋白含量的候选基因18 个，研究可为大豆7S 与11S 球蛋白含量相关QTL 精细定位及分子辅助育种提供参考。

1 材料与方法

1.1 大豆7S与11S球蛋白QTL信息收集

本研究中大豆7S 与11S 球蛋白相关QTL 收集自大豆数据库Soybase（http://www.Soybase.org/）。通过收集共获得国内外已报道的42 个7S 与11S 球蛋白相关QTL数据，整理其基本信息（见表1），用于后续元分析。

表1 大豆7S与11S球蛋白相关QTLs信息Table 1 Information on QTLs related to soybean 7S and 11S globulins

1.2 QTL映射

本研究通过Phytozome（https://phytozome-next.jgi.doe.gov/）网站获得Williams 82 物理图谱，其具有大量与QTL 原始图谱共同的标记，常被作为参考图谱整合不同原始图谱上的QTL。通过BioMercator ver.2.1软件[16]将大豆7S与11S球蛋白含量相关原始QTL 映射到Williams 82 物理图谱上，在不影响QTL 前提下，去除映射过程中反向标记，并使用下一个侧翼标记，同时剔除原始图谱和Williams 82 物理图谱不一致标记，提高分析过程准确性和可靠性。

1.3 大豆7S与11S球蛋白QTL元分析

元分析（Meta-analysis）用于定位真实QTL位置并确定其精确有效的置信区间。Arcade 等于2004年提供BioMercator ver.2.1 软件中元分析程序详细信息[16]。该分析产生5个模型，根据赤池信息量准则（Akaike information criterion，AIC）值确定QTL最优模型，计算连锁群上QTL 数量。在该模型中，每个原始QTL在染色体上位置均决定真实QTL物理位置，并可根据以下公式计算方差：

其中，σi2为连锁群上第i个QTL位置方差。通过以下公式计算真实QTL 95%置信区间：

AIC 值取决于每个模型的模拟，最小AIC 值QTL模型最接近通用QTL，其方差用原始QTL平均贡献率表示。

1.4 大豆7S与11S球蛋白候选基因挖掘

基于通用QTL 物理位置，从数据库基因注释数据集中选择大量候选基因，该数据库主要包括表2所示数据库。

表2 基因预测数据库和QTL相关信息整合网址Table 2 Database of gene prediction and QTL related information integration URL

2 结果与分析

2.1 大豆7S与11S球蛋白性状QTL收集与整理

从Soybase数据库中搜集到42个7S与11S球蛋白相关性状QTL 信息，所涉及的种群类型包括F6和F:7RIL。通过区间映射法（IM）、复合区间映射法（CIM）、多区间映射法（MIM）和方差分析进行映射信息收集。将42 个7S 与11S 球蛋白相关性状原始QTL 映射到Williams 82 物理图谱上，覆盖D1a、N、C1、C2、A2、K、O、J、D2、G、L 和I 染色体。其中最大R2为44.16%，LOD值为2~4.9。

2.2 大豆7S与11S球蛋白性状QTL元分析

基于元分析共得到11个7S与11S球蛋白性状真实QTL（见表3），主要分布在D1a、N、C1、C2、O、D2、L及I染色体上，定位于A2、K、J和G上的原始QTL呈分散分布，因此该染色体上无可用的真实QTL。如图1所示，原始QTL被整合，通过元分析得到真实QTL区域。每条染色体上真实QTL包含2~4个原始QTL，其中仅有一个真实QTL（MQTL17）包含4 个原始QTL；3 个真实QTL 包含3 个原始QTL，另外7个真实QTL包含2个原始QTL。所有真实QTL 图谱距离为0.34~21.58 Mb，平均为4.61 Mb。在D1a上包含3个真实QTL，每个真实QTL均包含2个原始QTL，AIC值均为37.1，其中95%置信区间最小为36.04，R2值为7.11~7.74，平均R2值为7.39；在D2上包含两个真实QTL，两个真实QTL的AIC值均为32.46。第一个真实QTL包含3个原始QTL，定位在18.12 Mb，R2值为15.03，95%置信区间范围为14.61~21.63；第二个真实QTL 包含4 个原始QTL，定位在11.97 Mb，R2值为9.67，95%置信区间范围为11.49~12.44。其余6个真实QTL分别定位于另外6条染色体上。

图1 7S与11S球蛋白性状QTL整合与元分析Fig.1 Integration and meta-analysis of QTL for 7S and 11S globulin traits

?

2.3 大豆7S与11S球蛋白性状候选基因挖掘

通过对真实QTL 区段进行基因提取，共获得1969 个基因。基于京都基因和基因组百科全书（KEGG）及基因本体术语（GO）对所得基因进行富集分析（见图2），分析结果涉及多项相关分子功能，其中主要包括二磷酸腺苷结合、蛋白质二聚化活性、核糖体结构成分、核酸结合等功能。此外，通过KEGG数据库进行基因富集分析得到的主要信号通路包括剪接体通路、内质网上的蛋白质合成通路、核糖体通路、植物激素信号转导通路和氨基酸生物合成通路等。同时参考同源拟南芥基因功能筛选所得基因，最终预测18 个基因与大豆7S和11S 球蛋白形成直接或间接相关，这18 个基因均在7S 与11S 球蛋白合成或分解中发挥重要作用，且该候选基因的同源基因皆被证实与7S和11S球蛋白合成代谢途径相关。这一观察结果表明，预测的候选基因均参与调控大豆7S 与11S 球蛋白含量。

图2 Meta-QTL区段内基因功能分析Fig.2 Functional analysis of genes within Meta-QTL segments

高新起研究表明，贮藏蛋白合成主要经由内质网—高尔基体—液泡内膜加工这一经典方式[17]，其中大量转录因子如LAFL[18]（LEC1、ABI3、FUS3和LEC2）、DOF[19]和AP2[20]家族转录因子直接或间接参与调控贮藏蛋白积累。贮藏蛋白运输主要包括4 种分子机制：信号肽序列（Signal peptide，SP）、液泡分选决定因子（Vacuolar sorting determinaants，VSDs）、液泡分选受体（Vacuolar sorting receptor，VSR）和SNARE蛋白复合体[21]，最终通过筛选确定的18 个候选基因皆与上述机制相关，候选基因在拟南芥（Arabidopsis thaliana）中的同源基因及基因注释信息见表4。

?

Glyma.01G159900 所编码的蛋白被注释为高尔基体定位的GRIP结构域蛋白，在内体到高尔基体的运输中发挥重要作用[22]。同时，该基因在GO数据库中被注释为参与高尔基体内蛋白质运输过程。Glyma.01G116900 和Glyma.01G117700 被 注 释 为 与Kunitz 家族胰蛋白酶和蛋白酶抑制剂蛋白相关，而Kunitz 家族胰蛋白酶与贮藏蛋白密切相关，具有保护贮藏蛋白的功能[23]。根据GO数据库对该基因的预测功能表明，其对内肽酶活性存在负调控（GO:0004866），而内肽酶具有水解蛋白质的功能。Glyma.01G160600、Glyma.03G164800 和Glyma.06G203 600在拟南芥中的同源基因均参与植物液泡相关蛋白的合成代谢，且研究表明，液泡是贮藏蛋白合成途径之一[17]。Glyma.04G166200、Glyma.17G149800和Glyma.06G196300所编码的蛋白与植物v-SNARE家族蛋白相关（GO:0005794），研究表明v-SNARE家族蛋白参与贮藏蛋白运输[21]。同时Glyma.06G196 300在GO数据库中被预测为参与高尔基体上蛋白质合成及运输（GO:0005794），KEGG数据库显示其与细胞内膜上物质转运过程相关。Glyma.17G175500和Glyma.20G186200所编码蛋白为AP2转录因子家族蛋白，AP2转录因子家族蛋白直接或间接参与调控贮藏蛋白积累[24]。Glyma.06G207100、Glyma.06G191600和Glyma.04G158100在拟南芥中的同源基因与信号肽肽酶相关，信号肽是指导蛋白多肽链合成与穿膜转移的决定因素，而贮藏蛋白合成均在内质网上完成[21]。Glyma.01G164100 被注释在GO:0042147 通路上，该基因编码液泡蛋白分选35（VPS35）。在种子成熟阶段，球蛋白和白蛋白前体在粗面内质网上合成并由VPS分选，表明VPS35可能通过将液泡分选受体从液泡前体循环到高尔基复合体[25]，从而参与种子中将蛋白质分选到液泡的过程。同时还有研究表明VPS10、VPS5、VPS17、VPS26、VPS29和VPS35具有调节贮藏蛋白运输的功能[26]。Glyma.04G138700 编码2S 白蛋白超家族蛋白，2S白蛋白类在双子叶植物中是最主要的贮藏蛋白。Glyma.04G151500和Glyma.20G194500在GO数据库中均注释为参与液泡内外物质运输过程（GO:0007034），其编码蛋白可能影响贮藏蛋白积累与运输。

3 讨论与结论

大豆7S 与11S 球蛋白含量不仅对大豆营养品质和加工特性起关键作用，且影响大豆种子生长与发育。目前，有关大豆7S 与11S 球蛋白性状相关QTL 研究较多，然而这些QTL 来自不同群体，生长在不同环境条件下，并用不同分析方法进行测试，存在QTL 定位区间过大、定位区间重叠等问题。因此需整合原始7S 与11S 球蛋白性状QTL以确定准确位点。

随基因组测序完成，研究学者们构建比遗传图谱精度更高的高密度物理分子图谱，物理图谱可准确描述已知DNA 在基因组中位置。同时物理图谱具有降低数据丢失风险、简化数据处理步骤、加快试验过程等优点。Schmutz 等将物理图谱和遗传图谱与大豆全基因组序列相结合，创建染色体规模的大豆序列组装[27]。本研究运用大豆高密度物理图谱作为参考图，解决遗传图谱或物理图谱与作图种群之间无足够数量共同标记时的映射问题。并通过基因注释获得大量相关候选基因，其中包括玉米植株形态、产量、抗病性及耐低磷等相关性状；水稻产量；大豆油分含量、白霉菌病；小麦面包中谷物膳食纤维含量等。本研究中，通过元分析从42个原始QTL中获得11个真实QTL，并利用基因注释在1 969 个基因中筛选得到18 个参与调控大豆7S 与11S 球蛋白含量的候选基因，研究表明，注释的18 个候选基因均直接或间接与大豆贮藏蛋白合成及运输相关。

Qi 等整合不同环境下多个试验大豆种子油、脂肪酸及蛋白质相关QTL 数据，通过元分析得到11个Meta-qtl，结合RNA-seq技术对不同种子发育阶段差异表达基因进行比较和鉴定，最终确定7个与大豆油和蛋白质积累相关的关键基因，其中Glyma.01G164100 编码液泡蛋白分选相关蛋白，其参与种子发育后期贮藏蛋白运输，对大豆贮藏蛋白含量有直接影响[28]。Huang 等通过对192 个大豆蛋白质QTL 进行共线性分析得到6 个热点区域，从QTL 热点区间检测到13 个与种子蛋白相关候选基因，其中Glyma.01G164100 在高蛋白品种中高表达，结合富集分析结果预测该基因与大豆蛋白含量相关[2]。王亚英以中黄13 和黑农51 杂交后构建的F2 群体为研究对象，通过精细定位确定Glyma.06G203600与大豆生育期性状相关[29]，表明Glyma.06G203600 可能通过大豆生育期影响大豆蛋白含量。Peng 等采用RNA-seq 技术对不同种子发育阶段的差异表达基因进行比较和鉴定，发现Glyma.20G186200 在高油低蛋白和高蛋白低油材料中表达量存在显著差异，表明该基因可能与大豆蛋白质和油含量相关[30]。

综上，因7S 与11S 球蛋白相关QTL 分析的研究尚少，本研究收集的部分QTL 涉及群体数量较少，单个QTL 提供的信息较为有限，试验结果具有一定局限性，有待后续研究。