Journal of Northeast Agricultural University
东北农业大学学报第52卷第3期52(3):1~122021年3月
March 2021
straw大豆抗倒伏相关性状全基因组关联分析
李文滨，吴航，刘冀，张翔超，郭志文，郑立娜，赵雪，韩英鹏，滕卫丽
（东北农业大学大豆研究所，大豆生物学教育部重点实验室，农业农村部东北大豆生物学与遗传育种重点实验室，哈尔滨
150030）
摘要：研究以150份大豆种质为试验材料开展2年3点试验，在收获期测定抗倒伏相关性状，并利用RTM-
GWAS 方法作关联分析，共检测到99个显著关联SNP 位点，其中48个效应显著位点，解释0.28%~3.38%表型变异；89个与环境互作显著位点，解释0.53%~7.04%表型变异；其中，21个位点主效表型变异解释率高于1%；与6个倒伏相关性状显著关联SNP 位点中，与茎粗显著关联SNP 位点最多，
与抗倒指数显著关联SNP 位点最少；获得22个有功能注释基因，根据基因表达量和基因功能注释，推测3个基因（Glyma .18G 149700、Glyma .04G 244100、Glyma .18G 011400）可作为与大豆抗倒伏相关候选基因。
关键词：大豆；倒伏性状；全基因组关联分析；候选基因挖掘中图分类号：S565.1
文献标志码：A
文章编号：1005-9369（2021）03-0001-12
李文滨，吴航，刘冀，等.大豆抗倒伏相关性状全基因组关联分析[J].东北农业大学学报,2021,52(3):1-12.DOI ：10.19720/jki.issn.1005-9369.2021.03.001.
Li Wenbin,Wu Hang,Liu Ji,et al.Genome-wide association analysis of lodging-related traits in soybean[J].Journal of Northeast Agricultural University,2021,52(3):1-12.(in Chine with English abstract)DOI ：10.19720/jki.issn.1005-9369.2021.03.001.
Genome-wide association analysis of lodging-related traits in soybean/
LI Wenbin,WU Hang,LIU Ji,ZHANG Xiangchao,GUO Zhiwen,ZHENG Lina,ZHAO Xue,HAN Yingpen
g,TENG Weili (Institute of Soybean Rearch,Key Laboratory of Soybean Biology in Chine of Ministry Education,Key Laboratory of Soybean Biology and Breeding (Genetics)of Chine Agriculture and Rural Ministry,Northeast Agricultural University,Harbin 150030,China)
Abstract:In this study,150soybean germplasm were ud as experimental materials for 2-year
and 3-point test,lodging-related traits were measured at harvest time,and the association analysis was conducted by RTM-GWAS method.There were 99significantly associated SNP sites detected,48loci with significant effect accounted for 0.28%-3.38%of phenotypic variation,89loci with significant interaction with environment explained 0.53%-7.04%of phenotypic variation,21of them explained more than 1%of phenotypic variation.Among SNP loci significantly associated with six lodging related traits,the SNP loci associated with stem diameter were the most,and tho associated with lodging index were the least.A total of 22genes with functional annotation were obtained.According to gene expression and gene function annotation,three genes (Glyma .18G 149700,Glyma .04G 244100,and Glyma .18G 011400)could be considrerd as candidate genes related to soybean lodging resistance.
Key words:soybean;lodging trait;genome-wide association analysis;mining for candidate gene
基金项目：国家十三五重点研发计划项目（2017YFD0101306-05）；国家自然科学基金项目（31471517）；黑龙江省“百千万”工程科技重大专项（2019ZX16B01-1）；现代农业产业技术体系专项建资金项目（CARS-04-PS04）
作者简介：李文滨（1958-），男，教授，博士，博士生导师，研究方向为大豆遗传育种与生物技术。E-mail:*****************
东北农业大学学报第52卷
大豆是重要经济和粮食作物，其营养价值高，用途广泛。倒伏是限制大豆高产、稳产和优质重要因素之一，大豆品种抗倒伏性是获得高产稳产重要保证[1-2]。倒伏是受遗传、栽培、气候和生态等多种因素共同制约的复杂性状。以往研究表明作物倒伏与植株地上部性状（如株高、茎粗、茎秆强度）[3-4]和地下部性状（如根重、根长、根体积等）有关[5-6]。在关于大豆倒伏及其相关性状研究中，周蓉和钟开珍等认为大豆倒伏性与茎秆性状（株高、主茎节数、节间长度、分枝数等）和产量性状（单株荚数、单株粒数、单株产量和百粒重）显著相关[7-8]。由于大豆倒伏性与该性状密切相关，且目前对于解析大豆抗倒性遗传基础研究报道相对较少，利用分子标记技术探索大豆倒伏性遗传特性成为大豆科研工作者研究热点。通过QTL 定位方法已发掘出许多与大豆倒伏性相关位点，但利用QTL定位通常需构建作图群体，研究时间较长、定位精度低[9-10]，而基于自然群体的关联分析具有研究时间短，检测效率高、定位精度高等优势[11]。
诸多学者已对大豆主要农艺性状开展GWAS （Genome-wide association analysis）分析，例如，Jing 等通过GWAS检测到14个与株高相关SNP位点，预测6个可能与株高相关新基因[12]。Chang等对368份大豆种质资源作GWAS分析，检测到45个与株高相关位点和43个与主茎节数相关位点[13]。张军等对190份大豆种质多个农艺性状作GWAS分析，累计检测到136个位点与目标性状显著关联，其中
14个位点与株高关联，13个位点与倒伏性关联[14]。张友谊利用224份大豆微核心种质和1749个SNP标记作GWAS分析，共检测到36个SNPs位点，其中与株高关联的13个位点，主茎节数4个，茎粗9个[15]。王自力对由573份种质建立的江淮地区大豆育种群体作GWAS分析，两年重复检测到与倒伏性关联SNP标记2个，株高130个，主茎节数92个[16]。关联分析不仅可找到与大豆倒伏性状有关遗传位点，还可预测与目标性状相关新基因。因此利用大豆自然群体作关联分析研究抗倒伏性遗传特性，对实现大豆高产稳产优质等目标具有重要意义。
本研究利用150个大豆种质资源组成的自然群体开展多年多点试验，通过RTM-GWAS方法结合SNP位点对株高、主茎节数等性状作全基因组关联分析，挖掘候选基因，为分子标记辅助大豆抗倒伏育种提供依据。
1材料与方法
1.1试验材料
供试材料选用150份不同大豆种质资源构成的自然群体。其中国外品种（系）3份、北京品种（系）8份、黑龙江品种（系）93份、吉林品种（系）20份、辽宁品种（系）3份、内蒙古品种（系）15份、新疆品种（系）2份、河北品种（系）4份、山东品种（系）1份和甘肃品种（系）1份。
1.2试验设计
试验于2018~2019年在向阳（A）、呼兰（B）、阿城（C）3个试验地点开展，采用随机区组设计，3次重复，行长2m、垄宽0.65m，株距0.06m。成熟期收获并室内考种，田间管理同大田。
1.3方法
1.3.1性状调查
根据周蓉等方法于成熟期调查倒伏性[17]，田间倒伏分级标准：1级为小区植株全部直立；2级为轻度倒伏（植株倾斜<15°）；3级为重倒伏（植株倾斜15°~44°）；4级为严重倒伏（植株倾斜>45°）。从自然群体每个资源中选择具有代表性连续5株测定株高、主茎节数、节间长度、茎粗、重心高度、茎秆强度、单株荚数、单株粒数、单株粒重、百粒重、茎秆抗倒指数等性状，3次重复，以平均值为观测值。其中重心高度和茎秆抗倒指数测定均参照陈喜凤等方法[2]，其他室内考种项目及标准均参照邱丽娟等方法[18]。
fellow是什么意思
1.3.2基因型检测与分析
采用特异性位点扩增片段测序（SLAF-q）技术，利用Illumina基因组分析仪Ⅱ获得大豆染色体均匀分布高质量单核苷酸多态性（SNPs）。提取150份大豆种质资源基因型信息，使用Seqtk（ /lh3/qtk）过滤，最终得到最小等位基因频率（Minor allele frequency，MAF）>0.05的23309个标记用于进一步关联分析。标记覆盖所有染色体区段，标记之间平均距离为28kb。
1.3.3全基因组关联分析
how long
利用RTM-GWAS方法采用两阶段策略对大豆倒伏相关性状表型值作关联分析。
世界和平英文
··2
李文滨等：大豆抗倒伏相关性状全基因组关联分析第3期第一阶段模型公式如下：y i =μ+∑j =1
J
w ij αj +∑l =1
L
x il βl +εi
（1）
式中，y i 表示个体i 表型观测值；μ表示总体平均数；J 表示用于群体结构矫正的特征向量个数，w ij 表示遗传相似系数矩阵第j 个特征向量在个体i 上的系数，αj 表示第j 个特征向量效应；L 表示测验标记位点等位基因数目，x il 表示测验标记位点第l 个等位基因对于个体i 基因型指示变量，取值0或1；βl 表示第l 个等位基因效应；εi 表示假定服从正态分布残差效应。
第二阶段模型公式如下：y i =μ+∑j =1
J
w ij αj +∑k =1K
∑l =1
L
x ijk βkl +εi
（2）
式中，K 表示总QTL 数目；L k 表示第k 个位点等位基因数目，x ik l 表示第k 个位点第l 个等位基因在个体i 上基因型指示变量，取值0或1；βkl 表示第k 个位点第l 个等位基因效应。其他符号含义与模型（1）相同。1.4
数据处理与分析
采用Microsoft Office Excel 11.1.0.10314、SPSS
3.0软件分析表型数据和作图。基于Phytozome 数据库（phytozome/）和Soyba 数据库（www.Soyba.Org/）中基因功能注释预测与大豆抗倒
伏相关候选基因。在eFP 数据库（bar.Utoronto.ca ）查询候选基因表达量数据，利用R 软件CMplot 程序包以及pheatmap 程序包绘制曼哈顿图和候选
基因组织表达热图。
2
结果与分析
2.1
大豆倒伏相关性状表型数据分析
分析多环境下大豆自然群体各农艺性状相关
性（见表1）。结果表明，在不同环境下倒伏级别与各性状相关性表现不同，除2019年呼兰点外倒伏级别与株高均表现显著相关，相关系数为0.243~0.527；两年阿城点倒伏级别与主茎节数均表现显著相关，相关系数为0.171~0.221；2018年呼兰点倒伏级别与节间长度表现显著相关，相关系数为0.275；2019年向阳点倒伏级别与茎粗均表现显著
相关，相关系数为0.290；除2018年阿城点外倒伏级别与重心高度均表现显著相关，相关系数为-0.209~0.345；两年呼兰点倒伏级别与抗倒指数均表现显著相关，相关系数为-0.254~0.178，表明倒伏受植株株高、主茎节数、节间长度、茎粗、重心高度以及抗倒指数等性状影响显著，该性状均为倒伏显著相关性状。
分析多环境下大豆自然群体6个倒伏显著相关性状表型（见表2）。结果表明，不同环境下，大豆自然群体倒伏相关性状表现均不同。各倒伏相关性状在6个环境下偏度和峰度绝对值多数<1，且具有广泛变异，其中株高变异系数为0.13~0.29；主茎节数变异系数为0.11~0.23；节间长度变异系数为0.15~0.22；茎粗变异系数为0.14~0.23；重心高度变异系
数为0.14~0.23；抗倒指数变异系数为0.18~0.44。年份/地点
Year/Site 2018A 2018B 2018C 2019A
2019B 2019C
株高PH 0.280**0.400**0.243**0.445**0.0610.514**
主茎节数NOMS 0.1550.1100.171*0.157
-0.031
0.221**
节间长度IL 0.1000.275**
-0.087-0.039-0.058-0.099
茎粗(cm)
ST 0.174
-
0.1260.1270.290**0.0770.060重心高度WH 0.285**0.345**0.0750.202*
-0.209*
0.223**
茎秆强度(N)SI 0.143
-0.0600.0810.156-0.0320.069单株荚数PNP -0.027-0.1020.1340.053
-0.058-0.107
单株粒数SNP -0.0740.0780.0950.0820.006
-0.047单株粒重SWP -0.043-0.0310.0730.018
-0.115-0.154
百粒重
100-SW 0.069
-
0.050-0.054-0.015-0.1140.060抗倒指数LR
-0.003-0.254**0.073
-0.1050.178*
-0.147
表1
大豆自然群体倒伏级别与不同农艺性状相关系数
Table 1
Correlation coefficients between lodging grade and agronomic characters in soybean natural population
注：*与**分别表示0.05与0.01水平差异显著。株高-PH ；主茎节数-NOMS ；节间长度-IL ；茎粗-ST ；重心高度-WH ；茎秆强度-SI ；单株荚数-PNP ；单株粒数-SNP ；单株粒重-SWP ；百粒重-100-SW ；抗倒指数-LR 。下同。
Note ：*and **indicate significant difference at 0.05and 0.01levels ，respectively.PH-Plant height ；NOMS-Node on main stem ；IL-Inter⁃
node length ；ST-Stem diameter ；WH-Weight height ；SI-Stem intension ；PNP-Pod number per plant ；SNP-Seed number per plant ；SWP-Seed weight per plant ；100-SW-100-ed weight ；LR-Lodging resistance.The same as below.
·
·3
宁波雅思
东北农业大学学报第52卷
大豆自然群体6个倒伏显著相关性状在6个环境下频数分布情况见图1。结果表明，该群体各倒伏相关性状呈连续性变异，在群体中均呈正态分布，表明此性状为多基因控制数量性状，适合于关联分析。2.2大豆倒伏相关性状全基因组关联分析
利用23309个单核苷酸多态性位点基因分型信息，通过RTM-GWAS程序对自然群体6个倒伏相关性状作多环境关联分析（见图2）。
性状
Traits
株高(cm) PH，Plant height
主茎节数NOMS，Node on main stem
节间长度(cm) IL，Internode length
茎粗(cm) ST，Stem diameter
重心高度(cm) WH，Weight height
抗倒指数LR，Lodging resistance
年份
Year
2018
2018
2018
2019
2019
2019
2018
2018
2018
2019
2019
2019
2018
2018
2018
2019
2019
2019
2018
2018
2018
2019
2019
2019
2018
2018
2018
2019
2019understand什么意思
2019
2018
2018
2018
2019
2019
2019
地点
Site
A
B
C
A
B
C
A
B
C
A
B
C
A
B
Cuit
A
B
C
A
B
C
A
B
C
A
B
C
A
B
C
A
B
C
A
B
C
最大值
Max.
122.20
162.20
135.80
109.80
162.80
130.20
23.00
25.00
23.50
32.40
23.60
22.80
5.55
5.21
5.70
3.47
5.53
4.43
1.26
1.91
1.51
1.34
1.27
1.33
52.50
68.00
75.00
56.60
73.00
64.80
15.41
11.83
12.90
2.10
1.17
1.35
最小值
Min.
31.20
33.00
47.25
32.30
73.20
65.30
8.75
9.20
10.40
european9.80
12.60
13.20
1.67
1.67
2.17
1.74
2.11
1.78
0.42
0.48
0.50
0.48
0.58
0.55
17.00
20.00
25.00
16.90
28.80
31.20
1.51
1.38
1.18
0.69
0.34
0.26
平均值±标准差
Mean±SD
64.27±18.93
86.81±23.24
91.31±19.61
73.04±16.99
108.17±13.88
96.75±14.71
13.91±3.19
16.35±2.92
16.09±2.51
17.94±2.61
18.25±2.06
17.88±2.19
3.45±0.62
3.16±0.69
3.66±0.64
2.54±0.38
3.55±0.71
2.99±0.48
0.76±0.16
0.92±0.21
0.81±0.14
0.82±0.14
0.91±0.13
0.84±0.12
33.73±7.53
44.17±9.40
48.53±9.35
37.30±8.46
55.09±7.69
48.40±7.28
5.06±2.22
5.05±2.21
3.92±1.60
1.18±0.29
0.77±0.14
0.89±0.17
峰度
Kurtosis
-0.13
0.54
-0.42
-0.58
1.09
-0.59
-0.17
0.28
-0.01
5.69
-
0.03
-0.61
0.52
0.16
-0.13
-0.50
-0.02
0.33
0.52
2.99
4.86
1.53talk的过去式
-0.04
2.31
-0.40
-0.33
-0.25
-0.57
0.47
-0.64
3.91
0.50
6.28
0.68
0.76
1.34
偏度
Skewness
0.56
0.52
0.07
-0.29
0.47
0.17
0.53
-0.06
0.16
0.79
-0.02
0.09
0.21
0.61
0.18
0.27
toaster
0.13
0.26
0.57
1.09
1.33
0.71
0.13
0.93
0.29
-0.08
-0.03
-
0.10
-0.31
-0.10
1.52
0.84
1.55
0.88
0.03
-0.22
变异系数(%)
CV
0.29
0.27
0.21
0.23
0.13
0.15
0.23
0.18
0.16
0.15
0.11
0.12
0.18
0.22
0.17
0.15
0.20
0.16
0.21
0.23
0.17
0.17
0.14
0.14
0.22
0.21
0.19
0.23
0.14
0.15
0.44
0.44
0.41
0.25
0.18
0.19
表2大豆自然群体倒伏相关性状表型分析
Table2Phenotypic analysis of lodging traits in soybean natural population
··4
李文滨等：大豆抗倒伏相关性状全基因组关联分析第3期抗倒指数Lodging resistance
年份（地点）Year(Site)
2018A 2018B 2018C 2019A 2019B 2019C
频率（%）F r e q u e n c y
706050403020100
0~0.40
0.40~0.800.80~1.201.20~1.601.60~2.00 2.00~2.40
重心高度(cm)Weight height
年份（地点）Year(Site)
2018A 2018B 2018C 2019A 2019B 2019C
15~2525~3535~4545~5555~6565~75
频率（%）F r e q u e n c y
50454035302520151050
图1不同环境下大豆自然群体倒伏相关性状分布
Fig.1Distribution of lodging traits of soybean natural population in different environments
结果表明，在P <0.01显著性水平上，共检测到99个显著关联SNP 位点，且这些显著关联SNP 位点在
大豆染色体上呈不均匀分布（见表3）。检测到16个与株高显著相关SNP 位点，分布于1号、2号、3号、4号、6号、7号、8号、9号、13号、15号和18号染色体上，其中有9个效应显著位点总表型变异解释率为13.98%，14个与环境互作效应显著位点总表型变异解释率为19.47%，其中表型变异解释率高于1%位点5个。检测到19个与主茎节数显著相关SNP 位点，分布于3号、4号、7号、9号、12号、13号、15号、18号、19号和20号染色体上，其中有10个效应显著位点总表型变
异解释率为14.54%，16个与环境互作效应显著位点总表型变异解释率为34.45%，表型变异解释率高于1%位点7个。检测到18个与节间长度显著相关SNP 位点，分布于2号、3号、4号、5号、10号、14号、15号、17号、18号、19号和20号染色体上，其中有7个效应显著位点总表型变异解释率为5.92%，18个与环境互作效应显著位点总表型变异解释率为37.81%，表型变异解释率高于1%位点1个。检测到20个与茎粗显著相关SNP 位点，分布于2号、3号、4号、6号、8号、10号、11号、15
号、16号和19号染色体上，其中有10个效应显著位点总表型变异解释率为7.73%，18个与环境互作
节间长度(cm)Internode length
株高(cm)Plant height
主茎节数Node on main stem
茎粗(cm)Stem diameter
频率（%）F r e q u e n c y
454035302520151050
年份（地点）Year(Site)
年份（地点）Year(Site)
频率（%）F r e q u e n c y
6050403020100
2018A 2018B 2018C 2019A 2019B 2019C
30~4545~6060~7575~9090~105105~120120~135135~150
150~165
频率（%）F r e q u e n c y
9080706050403020100
年份（地点）Year(Site)
2018A 2018B 2018C 2019A 2019B 2019C
年份（地点）Year(Site)
2018A 2018B 2018C 2019A 2019B 2019C
7~10
10~1313~1616~1919~22
22~25
2018A 2018B 2018C 2019A 2019B 2019C
频率（%）F r e q u e n c y
706050403020100
0.40~0.600.60~0.800.80~1.001.00~1.201.20~1.401.40~1.60
1~22~33~44~55~6
·
·5
