大豆菌核病是20世纪80年代中期开始出现的世界范围大豆病害，主要分布在美国、加拿大、中国、阿根廷、匈牙利、日本和印度等国，特别是在北美和加拿大的寒凉大豆种植区已经成为影响大豆产量的主要病虫害之一[1-2]。它主要由核盘菌(Sclerotinia sclerotiorum(Lib.)de Bary)侵染寄主植株所引起。1986年，加拿大的豆类作物菌核病发病率达25%[3]。近年来，大豆菌核病有继续扩大蔓延、加重危害的趋势，严重威胁大豆的生产[4]。生产实践证明，培育抗(耐)病大豆品种是防治大豆菌核病最经济有效的途径。目前为止，还没有发现完全免疫菌核病的大豆资源，前人研究报道，大豆种质资源对菌核病的抗性存在着遗传差异，有的大豆资源具有部分抗性，大豆菌核病是由多基因控制的数量性状遗传。数量性状对环境条件敏感，在大田条件下目标性状的选择压力难以保证一致，导致筛选抗菌核病种质资源比较困难，因此，利用分子标记辅助选择培育抗性品种显得越来越重要。
大豆抗菌核病的QTLs定位自2000年Kim等[5]利用分子标记RFLP、SSR和RAPD进行研究以来，很多研究都是基于SSR分子标记对大豆抗菌核病进行QTLs定位，但由于一代和二代分子标记在基因组中的多态性有限，使其在QTLs定位中的利用不太理想，随着二代测序技术的发展，基于SNPs的全基因组关联研究在生物领域成为热点，自从全基因组关联研究在人类疾病研究中取得进展后，植物上通过全基因组关联研究进行病害候选基因定位也越来越受到重视。由于全基因组测序面临费用高和较多重复序列的缘故，为了提高测序性价比，基于限制性核酸内切酶识别位点相关DNA序列(RAD tags)的简化基因组测序技术得到了普遍应用，其中GBS(Genotyping by sequencing)技术由于能大幅降低基因组的复杂度，尤其适合于大样本量的研究，可快速鉴定出高密度的SNPs位点，其结果在QTLs定位中应用越来越广[6-12]。目前为止，在大豆上多选择Apek Ⅰ、Msp Ⅰ和Pst Ⅰ酶进行简化基因组后测序，基于SNPs的全基因组关联研究进行大豆菌核病QTLs定位的研究报道较少，由于定位样本量大小、表型鉴定方法、分析算法及标记多态性位点数目等因素，影响了QTLs定位的一致性。且由于定位材料及标记数目不同，研究结果也不一致[13-16] 。分子标记辅助选择培育大豆抗菌核病育种还未见报道，抗性种质资源的选择很多还是用草酸浸泡叶、根、茎后，根据植株萎蔫程度和茎中可溶性色素水平高低评价抗性强弱[17]，因此，大豆菌核病QTLs的定位对大豆抗性育种及抗性机理的研究进展相对缓慢。本研究通过对126个加拿大大豆品种的30 125个SNPs，利用PLINKv 1.07分析软件，进行SNP-phenotype和Haplotype-phenotype的关联分析，找到大豆与抗菌核病强关联的位点/候选基因，为大豆菌核病QTLs分析及大豆抗菌核病育种提供参考。
1 材料和方法
1.1 DNA提取、文库准备、测序、数据质量控制
126个大豆品种的DNA提取、文库准备、测序、质量控制等步骤参考文献[13]。
1.2 分析方法
1.2.1 表型值 用棉垫接种法接种核盘菌菌丝进行表型鉴定[18]，茎秆损伤长度介于28.6～192.4 mm，供试材料的表型值分布见图1。
图1 供试材料主茎损伤长度表型值分布
Fig.1 Distribution of phenotype on lesion length
1.2.2 数据质量控制与分析 数据质量控制是全基因组关联分析的关键之一，不同参数的选择代表了不同的遗传学意义，产生不同的分析结果，不当的参数选择会产生假阳性或假阴性结果。本研究位点间的关联程度r2值0.8，LOD取值3，置信区间为95%，MAF值取0.01。由于大豆为自交作物，偏离哈迪-温伯格平衡的SNPs很可能是与病害关联的易感位点，因此，本研究不考虑哈迪-温伯格平衡，得到30 125个SNPs(图2)。126个材料共有3.7958E6个位点，其中杂合基因型位点有248 456个，占6.546%，杂合体的质量控制在10%以下，说明了数据可信度高[19-21]。
图2 每条染色体上的SNP数
Fig.2 Number of SNPs on each chromosomes
2 结果与分析
2.1 126份大豆材料的NJ法聚类分析
用TASSEL 5.0的邻接法(The neighbor-joining algorithm，NJ)对126个供试材料进行聚类(图3)，聚为3大类，分别为67，45，14个材料。
图3 基于NJ法的聚类图
Fig.3 Clustering based on NJ using TASSEL 5.0
2.2 PC分析
用SPSS对126份大豆品种进行主成分分析，结果显示前2个主成分的方差贡献率分别为6.24%和5.44%，Ward法聚类结果见图4。126份供试材料聚为2类，分别为70，56个材料。
2.3 群体结构分析
很多研究表明，群体结构分析是对目标性状进行关联分析的前提，因为复杂的群体结构会导致基因型与表型之间的假阳性。研究发现，考虑群体结构和不考虑群体结构时检测到的与表型变异相关的多态性位点数不同[22]。本研究采用Structure 2.3.4软件模拟群体遗传结构[23]，其基本原理是：假设样本存在K个(K值未知)等位变异频率特征类型数，将供试材料归到K个群体，使得该群体位点频率都遵循同一个Hardy-Weinberg 平衡。设定群体数K的取值为1～10，将MCMC(Markov chain Monte Carlo)开始时的不作数迭代(Length of burn-in period)设置为10 000 次，再将不作数迭代后的MCMC设置为10 000 次，每个K运行20次，然后依据ΔK值最大的原则选取一个合适的K 值。本研究的ΔK值最大时的K值为2，见表1。将126个供试材料划分为2个亚群，分别为68，58个材料。
图4 Ward法聚类图
Fig.4 Clustering by Ward method using SPSS
用Structure、PCA和NJ法聚类结果有差异，特别是NJ法与其他2种方法差异较大，见表2。聚类一致度用Kappa检
这里P(A)为观察到的一致度，P(E)为由于随机性产生的期望一致率；当K< 0.4时，说明一致度不够理想；K≥0.75，说明一致度相当满意。Structure和PCA 聚类的K=0.44，这2种方法聚类结果中度一致。
表1 群体结构分析
Tab.1 Structure analysis of the number of population for K
K设置范围# K重复RepsLnP(K)均值Mean LnP(K)LnP(K)标准差Stdev LnP(K)Ln′(K)Ln″(K)ΔK120-2 589 682.505 071.628 5NANANA220-2 483 401.195 0138.781 9106 281.310 00014 089.035 000101.519 245320-2 391 208.920 0321.946 392 192.275 00016 784.045 00052.133 057420-2 315 800.690 0345.342 375 408.230 00025 109.750 00072.709 754520-2 265 502.210 04 893.658 250 298.480 0004 810.175 0000.982 941620-2 220 013.905 05 965.269 045 488.305 00057 011.050 0009.557 163720-2 231 536.650 0182 438.838 4-11 522.745 00058 591.785 0000.321 159820-2 184 467.610 088 245.016 747 069.040 000524 916.440 0005.948 398920-2 662 315.010 01 367 029.724 5-477 847.400 000930 500.375 0000.680 6731020-2 209 662.035 0354 997.861 7452 652.975 000NANA
表2 用Structure、PCA 和UPGME方法的聚类结果比较
Tab.2 Comparison of the clustering using Structure，PCA and UPGME
聚类方法Clustering method第1组代号No. of group 1第2组代号No. of group 2第3组代号No. of group 3结构Structure1,2,3,5,7,9,11,12,13,14,16,18,19,22,30,32,37,38,39,41,42,44,48,50,52,54,57,58,59,61,66,67,71,72,73,74,76,77,78,79,83,87,88,92,93,97,98,99,100,101,102,103,104,106,107,108,109,110,112,114,115,117,118,120,121,122,123,124(68份)4,6,8,10,15,17,20,21,23,24,25,26,27,28,29,31,33,34,35,36,40,43,45,46,47,49,51,53,55,56,60,62,63,64,65,68,69,70,75,80,81,82,84,85,86,89,90,91,94,95,96,105,111,113,116,119,125,126(58份)主成分分析PCA1,2,3,6,7,11,12,13,14,16,17,18,19,21,22,24,26,27,28,30,31,32,37,38,39,41,42,44,45,47,54,56,57,58,62,65,67,68,69,71,72,74,75,76,77,78,79,80,83,84,87,88,89,93,97,98,99,101,102,103,104,107,112,114,115,117,120,121,123,124(70份)4,5,8,9,10,15,20,23,25,29,33,34,35,36,40,43,46,48,49,50,51,52,53,55,59,60,61,63,64,66,70,73,81,82,85,86,90,91,92,94,95,96,100,105,106,108,109,110,111,113,116,118,119,122,125,126(56份)邻接法NJ1,3,7,9,11,14,18,19,22,23,30,33,34,37,38,41,42,46,48,50,52,55,57,58,59,60,61,63,66,67,70,71,73,74,77,78,79,85,86,88,90,93,94,98,99,100,101,102,103,104,105,106,107,108,109,112,113,114,116,117,118,119,120,122,123,124,126(67份)2,4,5,8,12,13,15,16,17,21,24,26,27,28,31,32,35,36,39,40,44,45,47,54,56,62,64,65,68,69,72,75,76,80,83,84,87,89,91,95,96,97,111,115,121(45份)6,10,20,25,29,43,49,51,53,81,82,92,110,125(14份)
2.4 连锁不平衡分析
用PLINKv 1.07命令“--r2”运行LD，过滤掉r2小于0.2的数据，共有82 911对，r2平均值为r0.73，r2大于等于0.8的有46 655对(占56.27%)，2等于1的有9 391对(占11.33%)。2个SNPs位点间相关(r2大于0.2)的平均距离为114 979.2 bp，r2小于0.8的平均距离为124 596.9 bp，r2大于等于0.8的平均距离为107 505.3 bp，r2等于1的平均距离为80 397.08 bp。两位点间距离≤500 kb(遗传学上把位点间大于500 kb认为没有连锁关系)的r2平均值为0.73。由此看出，LD随着2个位点间的距离增大而衰减。从单条染色体和全基因组r2和2个SNPs位点间距离的关系图也可以看出随着物理距离的增加，2个位点间的r2越小，即连锁强度越小，甚至没有(图5，6)。
图5 一号染色体上的LD衰减图
Fig.5 LD decay plot of chromosome 1 in soybean
图6 全基因组上的LD衰减图
Fig.6 LD decay plot of genome-wide in soybean
2.5 单体型分析
SNPs位点并不是独立遗传的，而是在染色体上倾向于以一个整体遗传给后代，成组遗传的SNPs位点在一代又一代的遗传中绝少发生重组，于是，这样的一组SNPs位点类型也就是单体型，也就是说相邻位点是以单体型为单位传递的(前提是位点间的LD很强)，当然，单体型分析在统计学上的作用主要是减少多个位点组合分析时的自由度，从而增加研究的效能，例如，如果要分析3个二等位多态位点的组合和疾病/病害的关系，理论上有9种不同的组合，但是因为这3个位点间不是独立的，在群体中可能只存在3种常见的单体型，如果3个位点间独立，也就是连锁平衡，那么群体中的单体型就会有9种，这时候做单体型分析是没有意义的，因此，进行单体型分析时最好以单体域Block为基础。在PLINKv 1.07里有2种单体域推断方法，Sliding windows和Haplotype list，Haplotype list要先计算Blocks，然后利用算出的Haplotype和Phenotype之间的关联程度进行分析，这种运算比Sliding windows灵活，因为Sliding windows只能按SNP顺序依次排列分析Haplotype，不能按实际跳跃分析。在比较了2种方法后，本研究采用“--blocks”命令推断Blocks。20条染色体共有2 863个Haplotype，它们在20条染色体上的分布见图7，其中18号染色体上的Haplotype最多(230个)，11号染色体上的最少(57个)。一个Haplotype最多的由36个SNPs组成，最少2个SNPs组成。最长Haplotype跨越200 kb(因此，后续候选基因在200 kb范围内注释)，最短跨越2 bp。
2.6 基因型与表型关联分析
在基因型-表型的关联中，为了发现较多的关联信息，把P小于1E-5的分析结果都列表显示。
2.6.1 基于Haplotype-phenotype关联分析 基于Haplotype-phenotype关联分析中，17号染色体上的4个SNPs组成的Haplotype解释表型变异最大(17.56%)。其余依次为1，3，4，20号染色体上的Haplotype(表3)。
图7 每条染色体上的Haplotype数
Fig.7 Number of Haplotypes on each chromosomes
2.6.2 基于单个SNP-phenotype关联分析 在单个SNP的基因型-表型关联分析中，α≤0.05(P≤0.05/30125=1.660E-6)时，最强关联在3号染色体的34387780位置和其临近的几个SNP位点(P值为8.669E-7)，可解释表型变异的17.80%。在单个SNP-phenotype的关联分析中发现，往往多个SNPs组成了一个Haplotype，如在3号染色体上，从SNP3818到SNP3823一共6个SNP组成了一个Haplotype，这种情况同时也在其他染色体上看到。但由于126个材料中某些材料在一个Haplotype的几个SNPs中的个别SNP位点为杂合子，所以导致单个SNP解释的表型变异不一致(表4)，但峰值SNP都在相应的Haplotype中。其次关联在20，1，4，17号染色体上(图8)。
2.7 候选基因
由于峰值SNP都出现在相应的Haplotype中，因此结合单个SNP和Haplotype的基因型-表型关联分析位点，本研究将较强关联的位点在Soybeanbase数据库(http：//soybase.org/tools.php)中比对分析，可能参与大豆菌核病防御机制的蛋白质见表5，大致可分为参与生物过程、细胞组成、分子功能及未知功能蛋白。
图8 基于单个SNP的基因型-表型关联
Fig.8 Association based on SNP-phenotype
表3 用PLINYKv 1.07进行的基于Haplotype-表型关联分析
Tab.3 Association analysis based on Haplotype-phenotype using PLINKv 1.07
染色体Chr.单体型物理位置 Position of haplotype 单体型(1=A, 2=C, 3=G, 4=T)Haplotype(1=A, 2=C, 3=G, 4=T)回归R平方R2T-统计STATP值P value175575883/5647814/5648648/573489731140.175 6-5.1401.038E-615589867/5700523/5724122/572414012230.164 8-4.9462.406E-6334371310/34371316/34371324/34385919/34387780/34387823/34387841/34387945/34387962/34395745/34395833/34407540/34446090/3445315912242411111310.158 64.8353.863E-646353086/6353873240.151 24.7006.810E-62042068110/42080972/42091969/42100739/42100805/42118002/421229080.150 9-4.6946.968E-6175575883/5647814/5648648/57348970.149 14.6618.011E-646353086/63538730.148 2-4.6458.535E-615589867/5700523/5724122/57241400.146 94.6209.455E-6
表4 用PLINYKv 1.07进行的基于单个SNP-表型关联分析
Tab.4 Association analysis based on single SNP-phenotype using PLINKv 1.07
染色体Chr.物理位置PositionSNP编号No. of SNP等位基因(抗病性/感病性)Allele(Resistant/Susceptible)回归R平方R2P值P value3343877803818T/C0.178 08.669E-73343878233819C/T0.178 08.669E-73343878413820G/A0.178 08.669E-7204209196929734T/C0.172 81.295 E-615589867219A/C0.172 51.328 E-6204211800229737C/T0.168 91.751E-63343957453823G/A0.168 41.823E-63343879453821G/A0.163 22.721E-6204210073929735T/G0.158 93.776E-6204212290829738T/C0.156 94.405E-63343879623822G/A0.150 96.941E-6463538734588C/T0.150 67.145E-617573489723593T/G0.149 18.011E-6463530864587T/C0.148 88.143E-615700523220C/G0.146 99.455E-6204292862429787T/G0.146 59.717E-6
表5 候选基因
Tab.5 Candidate genes
染色体Chr.SNP位点SNP position基因名称 Feature name 开始位置Start结束位置End基因注释 Annotation description 15589867Glyma.01g04800055895705601987功能未知的蛋白质(DUF1666)334387780/34387823/34387841/34387945/34387962Glyma.03g1291003438446034389538吡咯啉-5-羧酸还原酶334387780/34387823/34387841/34387945/34387962Glyma.03g1292003438466134388525细胞色素P450、家族86、亚家族A、多肽1334395745/34395833Glyma.03g1293003439468434396542S-腺苷-L-甲硫氨酸依赖性甲基转移酶超家族蛋白334407540/Glyma.03g1295003440400034408882真核rpb5 RNA聚合酶亚基家族蛋白334446090/34453159Glyma.03g1298003444604834446285未知蛋白质334453159Glyma.03g1299003444903034455831植物蛋白1589功能未知46353086/6353873Glyma.04g07590063506666354012Peptidyl-tRNA hydrolase Ⅱ(PTH2)家族蛋白1337847000Glyma.13g2769003784095137847094功能未知的蛋白质 (DUF616)1337886116Glyma.13g2775003788386337887320功能未知的蛋白质 (DUF2921)175575883Glyma.17g07130055736715577724未知蛋白质175647814/5648648Glyma.17g07220056462585651630纤维素合酶家族蛋白175734897Glyma.17g07330057326325736607信号识别粒子受体α亚基家族蛋白2042068110Glyma.20g1827004206754842070856未知蛋白质2042080972/42091969/42100739/42100850Glyma.20g1828004207333742102170Exocyst复合体亚基82042118002/42122908Glyma.20g1831004211219942123209FAD依赖性氧化还原酶家族蛋白
3 讨论
核盘菌侵染寄主的致病机理非常复杂，目前认为核盘菌分泌的细胞壁降解酶(Cell wall degrading enzymes)和草酸在致病中起重要作用，由于草酸氧化酶(Oxalate oxidase，OXO)和草酸脱羧酶(Oxalate decarboxylase，OXDC)能把草酸降解成CO2和H2O2，从而减少草酸的危害，降低了病原菌的侵染，研究表明把从小麦和大麦中分离得到的OXO(一种类萌芽素蛋白 germin-like proteins)和OXDC转化进大豆中能提高大豆对菌核病的抗性[25-27]。研究发现，核盘菌致病仅有细胞壁降解酶或草酸还是不能完全致病，需要其他因子协同才能具有致病力。大豆抗菌核病的许多研究都是通过鉴定接种菌丝体后，不同时间段大豆植株抗逆酶谱/防御相关基因的表达谱的变化和细胞壁胼胝质沉积来鉴定有哪些抗逆酶参与了此过程，前人用蛋白组学技术鉴定出超氧化物歧化酶(SOD)、多酚氧化酶(PPO)、过氧化物酶(POD)、苯丙氨酸解氨酶(PAL)等在菌核病侵染过程中抗感材料的活性有所变化[28]。本研究中没有关联出OXO、OXDC、SOD、PPO、POD、PAL等基因。
本研究基于Haplotype-Traits的最大关联在17号染色体上，在此区域内的纤维素合酶基因(Glyma.17g072200)被鉴定出关联。植物细胞壁主要由果胶质、纤维素、木质素等组成，是病原物侵入寄主植物的主要屏障。植物没有强大的免疫系统，所以自身的防御机制就显得尤为重要，细胞壁是植物抵御外来侵染的第一道屏障，在油菜中有研究报道，菌核病菌菌株的致病力与其纤维素酶活性呈正相关[29-30]。另外，Valera[31]从生理、解剖和分子技术对大豆菌核病的侵染过程进行了研究，感染不同阶段的大豆转录组测序研究证明感性品种的木质素生物合成与JA/ET-related defense responses关闭有关，证明了细胞壁特性与木质素酶在大豆抗菌核病中的重要作用。
基于单个SNP-Trait的最大关联在3号染色体上，在此范围内的候选基因有pyrroline-5-carboxylate(P5C)reductase(Glyma.03g129100)和cytochrome P450、family 86、subfamily A、polypeptide 1(Glyma.03g129200)。P5CR是广泛存在于原核和真核生物中的一种重要的管家蛋白，其主要功能是催化脯氨酸生物合成的最后一步反应，即在NAD(P)H的作用下，将吡咯啉-5-羧酸(P5C)转化为脯氨酸的反应。许多研究报道，在应激引起的损害中，脯氨酸含量在植物体内聚集[32-33]。植物中P450在寄主-病原物互作中的表达已经有较多证据[34-35]。
用全基因组SNP关联研究进行大豆菌核病QTLs定位显示出诸多优点，但QTLs定位工具的选择是QTLs定位研究中需要考虑的主要问题之一，在本研究中利用几种模型分析后发现，并不能断言某个软件或某种方法最适合QTLs定位，因此，在关联分析中尽可能多用几种软件或多种方法分析比较，然后对共同定位的位点进行进一步研究验证。
在与复杂性状的关联分析中，样本数和标记数多少影响定位效果，从理论上讲样本数目越多且标记覆盖基因组范围越大，其分析结果越客观。
本研究引入Haplotype-phenotype的关联可能提高真正与病害相关位点的检出率，从Haplotype关联的基因可以看出，与病害相关基因簇都包含在一个关联的Haplotype中，这些基因簇在病害与寄主应答中可能协同作用。通过对关联出的基因簇的分析，可以加深对病害应答通路的认识，Haplotype很可能作为一个功能单位在起作用。另外，在单个SNP和Haplotype与表型的关联分析中，峰值SNP都包含在其相应的Haplotype中，但由于Haplotype包含许多SNPs的遗传信息，采用单体型比单个SNP具有更好的统计分析效果，能降低自由度。
以α≤0.05为阈值，3号染色体上的候选基因为Glma.03g129100、Glma.03g129200、Glma.03g129300、Glma.03g129500、Glma.03g129800、Glma.03g129900。17号染色体上为Glma.17g071300、Glma.17g072200、Glma.17g073300。目前，通过关联分析在大豆的抗菌核病中也已经定位了一些候选基因，如何在抗病育种中有效应用这些QTLs，需要在分析的基础上进一步结合蛋白组学进行功能验证，直到成为稳定的标记才能在抗病育种中应用。
[1] Hoffman D D，Hartman G L，Mueller D S，Leitz R A，Nickell C D，Pedersen W L.Yield and seed quality of soybean cultivars infected with Sclerotinia sclertiorum[J]. Plant Dis，1998，82(7)：826-829.doi： 10.1094/PDIS.1998.82.7.826.
[2] Ploper L D.Management of economically important diseases of soybean in Argentina[C]//World soybean research conference Ⅵ Chicago，USA：University of Illinois Press，1999：269-280.
[3] Tu J C. Integrated disease control of white mould(Sclerotinia sclerotiorum)in navy bean(Phaseolus vulgaris)[J]. International Symposium on Crop Protection，1986，51(2b)： 731-740.
[4] Peltier A J，Bradley C A，Chilvers M I，Malvick D K，Mueller D S，Wis K A. Biology，yield loss and control of sclerotinia stem rot of soybean[J]. Journal of Integrated Pest Management，2012，3(2)：1-7.doi： 10.1603/ipm11033.
[5] Kim H S，Diers B W. Inheritance of partial resistance to sclerotinia stem rot in soybean[J]. Crop Science，2000，40(1)：55-61.doi： 10.2135/cropsci2000.40155x.
[6] Lai G R，Zhang J，Liu H，Dong X M，Lai J S，Huang Y Q，Chen J T，Guo J J. Construction of high density genetic map via GBS technology and QTL mapping for nutritional quality traits in maize(Zea mays)[J]. Journal of Agricultural Biotechnology，2017，25(9)： 1400-1410. doi： 10.3969/j.issn.1674-7968.2017.09.003.
[7] Daverdin G，Johnson-Cicalese J，Zalapa J，Vorsa N，Polashock J. Identification and mapping of fruit rot resistance QTL in American cranberry using GBS[J]. Molecular Breeding，2017，37：38.doi： 10.1007/s11032-017-0639-3.
[8] Lin M，Cai S B，Wang S，Liu S B，Zhang G R，Bai G H.Genotyping-by-sequencing(GBS)identified SNP tightly linked to QTL for pre-harvest sprouting resistance[J]. Theoretical and Applied Genetics，2015，128(7)：1385-1395.doi： 10.1007/s00122-015-2513-1.
[9] Sonah H , O′Donoughue L , Cober E , Rajcan I, Belzile F. Identification of loci governing eight agronomic traits using a GBS-GWAS approach and validation by QTL mapping in soyabean[J]. Plant Biotechnology Journal，2015，13(2)：211-221.doi： 10.1111/pbi.12249.
[10] Balsalobre T W A，Pereira G D S，Margarido G R A，Gazaffi R，Barreto F Z，Anoni C O，Cardoso-Silva C B，
Costa E A，Mancini M C，Hoffmann H P，Souza A P，Garcia A A F，Carneiro M S.GBS-based single dosage markers for linkage and QTL mapping allow gene mining for yield-related traits in sugarcane[J]. BMC Genomics，2017，18(1)：72.doi： 10.1186/s12864-016-3383-x.
[11] Das S, Upadhyaya H D, Bajaj D, Kujur A, Badoni S, Laxmi, Kumar V, Tripathi S, Gowda C L L, Sharma S, Singh S, Tyagi A K, Parida S K. Deploying QTL-seq for rapid delineation of a potential candidate gene underlying major trait-associated QTL in chickpea[J]. DNA Research，2015，22(3)：193-203. doi： 10.1093/dnares/dsv004.
[12] Pan L，Wang N，Wu Z H，Guo R，Yu X L，Zheng Y，Xia Q J，Gui S T，Chen C Y. A high density genetic map derived from RAD sequencing and its application in QTL analysis of yield-related traits in Vigna unguiculata[J]. Frontiers in Plant Science，2017，8：1544.doi： 10.3389/fpls.2017.01544.
[13] Bastien M，Sonah H，Belzile F. Genome wide association mapping of Sclerotinia sclerotiorum resistance in soybean with a genotyping-by-sequencing approach[J]. The Plant Genome，2014，7(1)：1-13.doi： 10.3835/plantgenome2013.10.0030.
[14] Iquira E,Humira S,Francois B. Association mapping of QTLs for sclerotinia stem rot resistance in a collection of soybean plant introductions using a genotyping by sequencing(GBS)approach[J]. BMC Plant Biology，2015，15(1)： 5. doi： 10.1186/s12870-014-0408-y.
[15] Zhao X, Han Y P, Li Y H, Liu D Y, Sun M M, Zhao Y, Lü C M, Li D M, Yang Z J, Huang L, Teng W L, Qiu L J, Zheng H K, Li W B. Loci and candidate gene identification for resistance to Sclerotinia sclerotiorum in soybean(Glycine max L. Merr.)via association and linkage maps[J].The Plant Journal，2015，82(2)：245-255. doi： 10.1111/tpj.12810.
[16] Moellers T C. Genome-wide association and epistasis studies of Sclerotinia sclerotiorum resistance in soybean[D].USA:Iowa State University，2016： 15039.
[17] 厉志，王曙明，孟凡梅，刘佳，衣志刚，张红艳，董志敏.茎中可溶性色素水平法鉴评大豆资源对菌核病的部分抗性[J].大豆科学，2016，35(3)：481-488. doi： 10. 11861 /j. issn. 1000-9841. 2016. 03. 0481.
Li Z，Wang S M，Meng F M，Liu J，Yi Z G，Zhang H Y，Dong Z M. Soluble pigment level in stems determining partly-resistance to white mold in soybean[J]. Soybean Science，2016，35(3)：481-488.
[18] Bastien M，Huynh T T，Giroux G，Iquira E，Rioux S，Belzile F. A reproducible assay for measuring partial resistance to Sclerotinia sclerotiorum in soybean[J]. Canadian Journal of Plant Science，2012，92(2)：279-288.doi： 10.4141/CJPS2011-101.
[19] Tabanqin M E,Woo J G,Martin L J. The effect of minor allele frequency on the likelihood of obtaining false positive[J]. BMC Proceedings，2009，3：S41. doi： 10.1186/1753-6561-3-S7-S41.
[20] Anderson C A，Pettersson F H，Clarke G M，Cardon L R，Morris A P，Zondervan K T. Data quality control in genetic case-control association studies[J]. Nature Protocols，2010，5(9)：1564-1573.doi： 10.1038/nprot.2010.116.
[21] Pongpanich M，Sullivan P F，Tzeng J Y. A quality control algorithm for filter SNPs in genome-wide association studies[J]. Bioinformatics，2010，26(14)：1731-1737.doi： 10.1093/bioinformatics/btq272.
[22] Price A L，Zaitlen N A，Reich D，Patterson N. New approaches to population stratification in genome-wide association studies[J]. Nature Reviews Genetics，2010，11(7)：459-463.doi： 10.1038/nrg2813.
[23] Earl D A，vonHoldt B M. Structure harvester： a website and program for visualizing structure output and implementing the evanno method[J]. Conservation Genetics Resources，2012，4(2)：359-361. doi： 10.1007/s12686-011-9548-7.
[24] Eugenio B D,Glass M. The Kappa statistic： a second look[J]. Computational Linguistics，2004，30(1)：95-101.doi： 10.1162/089120104773633402.
[25] Cobe E R, Rioux S, Rajcan I, Donaldson P A,Simmonds D. Partial resistance to white mold in a transgenic soybean line[J].Crop Science，2003，43(1)：92-95. doi： 10.2135/cropsci2003.9200.
[26] Cunha W G,Tinoco M L P，Pancoti H L，Ribeiro R E，Aragão F J L. High resistance to Sclerotinia sclerotiorum in transgenic soybean plants transformed to express an oxalate decarboxylase gene[J]. Plant Pathology，2010，59(4)： 654-660.doi： 10.1111/j.1365-3059.2010.02279.x.
[27] Calla B, Blahut-Beatty L, Koziol L, Zhang Y F, Neece D J, Carbajulca D, Garcia A, Simmonds D H, Clough S J. Genomic evaluation of oxalate-degrading transgenic soybean in response to Sclerotinia sclerotiorum infection[J]. Molecular Plant Pathology，2014，15(6)：563-575.doi： 10.1111/mpp.12115.
[28] 吕春梅，赵月，赵雪，王强，孟宪新，魏淑红，韩英鹏，李文滨，张俊华. 大豆品种Maple Arrow耐菌核病生化机制[J].中国油料作物学报，2014，36(5)：630-634.doi： 10.7505/j.issn.1007-9084.2014.05.011.
Lü C M，Zhao Y，Zhao X，Wang Q，Meng X X，Wei S H，Han Y P，Li W B，Zhang J H. Biochemical resistance mechanism of soybean cultivar Maple Arrow to Sclerotinia sclerotiorum[J]. Chinese Journal of Oil Crop Sciences，2014，36(5)：630-634.
[29] 刘小燕，高智谋，汪世军，凌萍.油菜菌核病菌致病力与纤维素酶活性关系的初步研究[J].现代农业科学，2009，16(6)：1-3.
Liu X Y，Gao Z M，Wang S J，Ling P. Preliminary study on the relationship between pathogenicity and cellulase activity of Sclerotinia sclerotiorum isolates[J].Modem Agricultural Sciences，2009，16(6)： 1-3.
[30] 毛玮，侯英敏，刘志文.核盘菌和草酸诱导下的油菜几种酶活力的变化分析[J].大连理工大学学报，2011，30(1)：39-42.doi： 10.19670/j .cnki .dlgydxxb.2011.01.010.
Mao W，Hou Y M，Liu Z W. Analysis of several enzymes activity induced by Sclerotinia sclerotiorum and oxalic acid in rapeseed[J]. Journal of Dalian Polytechnic University，2011，30(1)：39-42.
[31] Valera R E. Physiological，anatomical and molecular characterization of partial resistance against Sclerotinia sclerotiorum in soybean[D]. Guelph，Ontario，Canada： The University of Guelph，2014：1-171.
[32] Claussen W. Proline as a measure of stress in tomato plants[J].Plant Science，2005，168(1)：241-248.doi： 10.1016/j.plantsci.2004.07.039.
[33] 丁泽红，付莉莉，颜彦，铁韦韦，胡伟. P5CR基因的进化及其在木薯中的表达分析[J]. 生物技术通报，2018，34(3)：105-112.doi： 10.13560/j.cnki.biotech.bull.1985.2017-0796.
Ding Z H，Fu L L，Yan Y，Tie W W，Hu W. Evolution analysis of P5CR and expression analysis of P5CR gene in cassava[J]. Biotechnology Bulletin，2018，34(3)：105-112.
[34] 杨致荣，毛雪，杨致芬，李润植.细胞色素P450基因及其在植物改良中的作用[J].遗传，2003，25(2)：237-240. doi： 10.3321/j.issn：0253-9772.2003.02.027.
Yang Z R，Mao X，Yang Z F，Li R Z. Cytochrome P450 genes and their application in plant improvement[J]. Herditas，2003，25(2)：237-240.
[35] 戴素明，周程爱，谢丙炎，冯东昕，肖启明. 细胞色素P450表达在植物防御反应中的作用[J].石河子大学学报(自然科学版)，2004，22(S1)：184-187.doi： 10.3969/j.issn.1007-7383.2004.z1.060.
Dai S M，Zhou C A，Xie B Y，Feng D X，Xiao Q M. Role of cytochromes P450 expression in plant defence responding to pathogens[J]. Journal of Shihezi University(Natural Science)，2004，22(S1)：184-187.