Identification of Novel Genomic Regions for Bacterial Leaf Pustule (BLP) Resistance in Soybean (Glycine max L.) via Integrating Linkage Mapping and Association Analysis

Bacterial leaf pustule (BLP), caused by Xanthornonas axonopodis pv. glycines (Xag), is a worldwide disease of soybean, particularly in warm and humid regions. To date, little is known about the underlying molecular mechanisms of BLP resistance. The only single recessive resistance gene rxp has not been functionally identified yet, even though the genotypes carrying the gene have been widely used for BLP resistance breeding. Using a linkage mapping in a recombinant inbred line (RIL) population against the Xag strain Chinese C5, we identified that quantitative trait locus (QTL) qrxp–17–2 accounted for 74.33% of the total phenotypic variations. We also identified two minor QTLs, qrxp–05–1 and qrxp–17–1, that accounted for 7.26% and 22.26% of the total phenotypic variations, respectively, for the first time. Using a genome-wide association study (GWAS) in 476 cultivars of a soybean breeding germplasm population, we identified a total of 38 quantitative trait nucleotides (QTNs) on chromosomes (Chr) 5, 7, 8, 9,15, 17, 19, and 20 under artificial infection with C5, and 34 QTNs on Chr 4, 5, 6, 9, 13, 16, 17, 18, and 20 under natural morbidity condition. Taken together, three QTLs and 11 stable QTNs were detected in both linkage mapping and GWAS analysis, and located in three genomic regions with the major genomic region containing qrxp_17_2. Real-time RT-PCR analysis of the relative expression levels of five potential candidate genes in the resistant soybean cultivar W82 following Xag treatment showed that of Glyma.17G086300, which is located in qrxp–17–2, significantly increased in W82 at 24 and 72 h post-inoculation (hpi) when compared to that in the susceptible cultivar Jack. These results indicate that Glyma.17G086300 is a potential candidate gene for rxp and the QTLs and QTNs identified in this study will be useful for marker development for the breeding of Xag-resistant soybean cultivars.


Introduction
Bacterial leaf pustule (BLP), caused by the causal agent Xanthomonas axonopodis pv. glycines (Xag), is one of the most destructive foliar diseases in susceptible soybean cultivars. BLP causes a significant reduction in soybean yield and quality worldwide, especially in the soybean production areas in China [1], Korea [2], Thailand [3,4], the USA [5,6], and Benin [7], where seasonal high temperature and humidity favor Xag growth and development. Thus, BLP disease is more severe in South China, especially at Yangtze-Huai River Basins, in comparison with North China [1]. Recently, the incidence of BLP disease in South China has been rising along with global warming and the frequent occurrence of storms [1,8,9]. For example, the tropical storm 'Rumbia' caused an outbreak of BLP disease in the soybean growing area of Jianghuai River Basin in China in August 2018. In addition, it was reported that 88% of soybean growing in the fields in South Korea in 1998 were severely infected by Xag. Similarly, it accounted for 20.7-34.9% of soybean yield losses in Thailand in 1983 [2,10].
Xag infects soybean through natural openings and wounds on the abaxial surface of soybean leaves with the possibility to extend to stalks and petioles [5,11]. Typical BLP symptoms include small yellow-to-brown lesions with raised pustules in the lesion centers on infected soybean leaves, causing premature defoliation [5,11]. The pustule lesions caused by Xag can always be used as the infection sites by other microbial pathogens such as bacterium Pseudomonas tabaci, the causal agent of soybean wildfire disease [5]. Moreover, the BLP disease is easily misdiagnosed as Asian soybean rust caused by Phakopsora pachyrhizi [6], and the application of inappropriate pesticides could cause unnecessary economic burden to farmers.
To reduce economic losses in soybean production, disease resistance breeding is considered as the most cost-effective, efficient, and environmentally safe approach to achieve yield stability [12]. It was thought that a single major recessive resistance gene rxp confers BLP resistance while its dominant allele and associated genes regulate the degree of susceptibility [13,14]. The recessive rxp allele was originally determined in the soybean cv. CNS (P.I. 71,569), which is a highly BLP-resistant (near immune) cultivar and was developed in Nanjing, China [13]. To date, rxp has been widely used for the development of commercial BLP-resistant soybean cultivars worldwide [13,15]. In addition, rxp also confers resistance to soybean wildfire disease, which sometimes occurs with BLP disease [13].
The dominant Rxp allele, which confers susceptibility to the BLP disease, was originally found to be linked to the malate dehydrogenase (Mdh) locus with an about 16% recombination in a population of 650 F 2 soybean plants developed from the crosses of Clark 63 and P.I. 437,477B [16]. Then, the Rxp locus was mapped to the genomic region between markers Satt372 and Satt014 on chromosome (Chr) 17 by using the hybrid offspring of resistant parent Coker 237 and susceptible parent P.I. 97,100, as well as the resistant parent Young and susceptible parent P.I. 416,937 [5].
QTL fine mapping located the recessive rxp in the 33-kb-long genomic fragment between markers SNUSSR17_9 and SNUSNP17_12 on Chr 17, which contains three genes, i.e., serine/threonine/tyrosine-protein kinase HT1 coding gene Glyma17g09770, membrane protein coding gene Glyma17g09780, and zinc finger protein coding gene Glyma17g09790; the latter two were predicted to be the key candidate resistance genes [2]. A GWAS study identified one significant SNP in each of Glyma01g40560, Glyma01g40590, Glyma11g20310, and Glyma17g09801 that was associated with BLP resistance [17]. An RNA-Seq analysis identified the differentially expressed genes (DEGs) between resistant and susceptible soybean varieties under the infection of Xag strain 8ra at 0, 6, and 12 h post-inoculation (hpi) [18]. Interestingly, no significant difference was detected in the relative expression levels of Glyma17g09780 and Glyma17g09790 between the resistant and susceptible soybean varieties at 0, 6, and 12 hpi [18]. Several minor QTLs conferring resistance to various Xag isolates were also identified on Chr 18, 05, 04, 13, 19, and 10 in the recombinant inbred lines (RILs) derived from a cross between susceptible parent 'Suwonl57' and resistance parent 'Danbaekkong' using 76 SSR markers [14]. Furthermore, a novel QTL/gene responsible for BLP resistance was located at 21.5 cM away from the SSR marker Sat_108 on Chr 10 in P.I. 96,188, a soybean variety that only exhibits pustules without chlorotic haloes following Xag infection [19]. In addition, a few disease resistance genes were reported to be involved in BLP resistance in soybean. For example, the soybean LATERAL ORGAN BOUNDARY1 (GmLob1; Glyma.05g040500) was predicted to be targeted by the virulence factor (tal2b) of Xag as its homolog CsLOB1 in Citrus sinensis has been reported as a resistance gene to citrus canker disease caused by X. citri subsp. Citri (Xcc) [20,21]. A soybean lesion-mimic mutant NT302 containing a defective hydroperoxide lyase (HPL) gene Glyma.12g191400 (GmHPL) exhibited high susceptibility to Xag strain C5, indicating GmHPL's role in Xag resistance [22].
In the present study, we assessed the BLP resistance level in a soybean RIL population and an association panel of soybean breeding lines in multiple environments. Then, we applied linkage mapping and GWAS for the identification of the BLP resistance QTLs and quantitative trait nucleotides (QTNs) in the soybean genome, leading to the identification of novel QTLs/QTNs. Our linkage mapping and GWAS revealed that BLP resistance is governed by one single major QTL qrxp-17-2 and two minor QTLs qrxp-05-1 and qrxp-17-1 that were identified for the first time in the present study. The relative expression of five potential candidate genes in the resistant soybean cultivar W82 and susceptible cultivar Jack following Xag treatments showed that Glyma.17G086300, which is located in qrxp-17-2, is a candidate gene for rxp.

The BLP Resistance of the Soybean RIL and Association Panel Lines and Their Relationship with Flowering Time
BLP resistance was assessed in the soybean RIL population under the artificial inoculation condition and in the soybean association panel under both the artificial inoculation and natural morbidity conditions. The frequency distributions of the phenotypic data of soybean resistance to BLP disease exhibited a normal distribution pattern ( Figure S1, see Supplementary Materials). Under the artificial inoculation condition, the mean infection responses (IRs) of the RIL in 2014-2015 and the association panel lines in 2014-2016 varied from 4.08 ± 2.71 to 5.65 ± 2.62 and from 2.76 ± 1.54 to 3.64 ± 1.76, respectively (Table 1). However, the mean IRs of the association panel lines under the national morbidity condition in 2018 varied from 2.54 ± 1.72 to 4.33 ± 1.88 (Table 1), indicating the severity of the BLP disease under the natural morbidity condition in 2018 ( Figure 1). In addition, the BLP resistance levels in both populations showed a continuous variation from 1 to 9 (Tables S1-S3), demonstrating that the BLP resistance in these soybean lines is a quantitative trait and regulated by multiple genes. Analysis of variance (ANOVA) of the BLP resistance in each population under either condition (spray or the natural morbidity) indicated that the effects of genotype (G), environment (E), and genotype × environment (G × E) interactions were significant at the p < 0.001 level of significance (Table 1). The only exception came from the G × E interactions in the RIL lines under the spray condition, which were significant at the p < 0.05 level of significance. Moreover, the broad-sense heritability (h 2 ) for BLP resistance in the two populations under either condition was 79.81% or greater ( Table 1), indicating that the majority of phenotypic variation of BLP resistance in these soybean lines was attributed to effects of genotype (G).
Additionally, the flowering time of the soybean association panel was recorded under the natural morbidity condition in 2018 (Table S3) [23]. The correlation analysis showed a significant negative correlation at the p < 0.05 level (r = −0.29) between the BLP resistance level and flowering time, except the correlation was insignificant between the BLP resistance level and flowering time in the 2018DT-natural environment (Table S4).

Genetic Linkage Map Construction and Linkage Mapping in the Soybean RIL Population
A genetic map used for linkage mapping analysis was constructed in the soybean RILs. A total of 18.77 and 14.78 M reads and 1.58 and 1.24 Gbp bases of raw data were obtained for the two parents of the RIL population (Zhengyang and M8206), respectively. In addition, the mapped reads of two parents were 11.82 and 9.22 M and the mapped bases were 0.99 and 0.77 Gbp, respectively. The genome coverage of these two parents were 8.7 and 7.8%, respectively. For the RIL populations, a total of 66,677 SNPs were identified, and the average SNP density across the soybean 20 chromosomes was 6.85 SNPs per 100 kb (Table S5).  Using the composite interval mapping (CIM) method, we identified three QTLs, i.e., qrxp_5_1, qrxp_17_1, and qrxp_17_2, which were associated with BLP resistance in the RIL population in the 2014JP-and/or 2015JP-spray environments ( Table 2). Among them, qrxp_5_1 and qrxp_17_1 were detected in the 2015JP-spray environment and located in the physical intervals of 1-1,169,356 bp on Chr 05 and 5,158,677-5,994,063 bp on Chr 17, respectively. The qrxp_5_1 explained for 7.26% of phenotypic variation of BLP resistance with a negative effect, while qrxp_17_1 accounted for 22.26% of the phenotypic variation of BLP resistance with a positive effect ( Table 2). The qrxp_17_2 was detected in both environments and located in the physical interval of 6,777,393-6,883,854 bp on Chr 17, explaining for 74.33% of phenotypic variance in the 2014JP-spray environment, and located in 6,293,843-6,883,854 bp on Chr 17, accounting for 34.68% of the phenotypic variance in the 2015JP-spray environment ( Table 2). The qrxp_17_2 showed a positive effect to BLP resistance and was considered as the most stable QTL associated with BLP resistance in soybean ( Table 2).

Population Structure and Linkage Disequilibrium (LD) Analyses in the Soybean Association Panel
The 476 lines of the soybean association panel were genotyped using a high-density soybean array that consisted of 61,166 high-quality SNPs [24]. Clustering analysis using the neighbor joining method grouped all these lines into three groups ( Figure 2A). The principal component analysis (PCA) also grouped them into three groups ( Figure 2B). The first and second PCs explained for 17.27% and 5.02% of the total variance, respectively. The results of clustering analysis and PCA showed that the soybean lines derived from the same geographical locations and the same parental materials were generally grouped into the same subpopulations. The LD analysis showed that the r 2 decreased gradually as the physical distance increased, and the LD decay distance was estimated at 1.39 Mb, where r 2 dropped to being half of its maximum value of 0.74428 ( Figure 2C). The LD parameters, viz., "D'" and "r 2 " were calculated. The LD decay rate of the population was measured as the chromosomal distance when the average r 2 decreased to be half of its maximum value.

Association Mapping in the Soybean Association Panel
BLP resistance in the association panel was evaluated under the artificial inoculation condition by using Xag C5 in Jiangpu from 2014 to 2016. The GWAS analysis detected 38 QTNs that were significantly associated with BLP resistance under the recommended threshold of −log10(p) ≥ 4 and distributed on eight chromosomes ( Figure 3; Table S6). These QTNs explained for 3.26-6.04% of the phenotypic variation (Table S6). Among these QTNs, Gm17_7603802 and Gm17_7603992 on Chr 17 explained for the highest amount of phenotypic variation, i.e., 6.04%. A total of eight QTNs, i.e., Gm09_36501019, Gm17_7603802, Gm17_7603992, Gm17_7712768, Gm17_7721556, Gm17_7736150, Gm17_7754016, and Gm17_ 7754048, were identified in at least two environments (Table S6). BLP resistance of the association panel was also evaluated under the natural morbidity condition in Jiangpu and Dangtu in 2018. The GWAS analysis revealed that 34 SNPs were associated with BLP assistance at the suggested threshold of −log10(p) ≥ 4 and distributed on nine chromosomes ( Figure 3; Table S7). Among them, Gm17_7712768 explained for the phenotypic variation of 3.31-7.30% (Table S7). We found that 11 out of the 34 QTNs were co-localized with the QTNs identified in the same association panel under the artificial inoculation condition (Table 3). These 11 QTNs were Gm05_7667820, Gm05_7668047, Gm17_5628119, Gm17_5628133, Gm17_7603802, Gm17_7604008, Gm17_7712768, Gm17_ 7721556, Gm17_7736150, Gm17_7754016, and Gm17_7754048. A total of 13 and 11 stable QTNs associated with BLP resistance were identified in the association panel in ≥2 environments under both the artificial inoculation and natural morbidity conditions, respectively (Table 3; Tables S6 and S7). On the basis of the LD distance, we extended the significant associated region to cover 1.39 Mb upstream and downstream of the most significantly associated QTN positions (Table 3).

Co-Location Regions between Linkage Mapping and GWAS Analysis
The three QTLs identified by linkage mapping and the 11 stable QTNs identified by GWAS were co-localized in three genomic regions (Tables 2 and 3; Figure 4). QTL qrxp_5_1 and two stable QTNs Gm05_7667820 and Gm05_7668047 were localized in genomic region 1, i.e., the 1-7,668,047 bp region on Chr 05. QTL qrxp_17_1 and QTNs Gm17_5628119 and Gm17_5628133 were localized in genomic region 2, i.e., the 5,158,677-5,994,063 bp region on Chr 17. QTL qrxp_17_2 was co-localized with eight stable QTNs in the 6,777,393-7,754,048 bp region on Chr 17 of genomic region 3. Thus, these genetic regions are closely linked to the causal effects of variations in the BLP resistance in soybean, making them robust regions for the identification of candidate genes. a p, the statistical p-value for the significance of the odds ratio in the GWAS. b r 2 (%), phenotypic variance (%) explained by the peak markers. Bold QTNs indicate that the markers can be located by GWAS under both artificial inoculation and natural morbidity conditions. The significant associated region was 1.39 Mb upstream and downstream of the most significantly associated QTN positions.

Effects of the Most Significant Alleles in Three QTNs Individually or in Combination on BLP Resistance in Multiple Environments
Among the 11 stable QTNs identified by GWAS, Gm05_7667820, Gm17_5628119, and Gm17_7603802 explained for the highest phenotypic variation in the three co-localized genomic regions, respectively (Table 3). Thus, we assessed the effects of the most significant alleles G/A, T/C, and T/C in the QTNs Gm05_7667820, Gm17_5628119, and Gm17_7603802, respectively, on BLP resistance in multiple environments. The disease in-dexes of soybean cultivars containing Gm05_7667820 (G) and Gm17_7603802 (T) were significantly lower than that of soybean varieties containing Gm05_7667820 (A) and Gm17_7603802 (C) ( Figure 5). Moreover, the disease indexes of the resistant soybean cultivars containing the three resistance alleles GTT were significantly lower than that containing the three susceptible alleles ACC ( Figure 5). As a result, the cumulative effect of the three resistant alleles was able to increase the resistance for BLP disease in resistant soybean.

Prediction of Candidate BLP Resistance Gene in Soybean
Among the three BLP resistance QTLs and the 11 stable BLP resistance QTNs colocalized in three genomic regions, the genomic region 3 covers the single major recessive resistance gene rxp. To further analyze the candidate genes of rxp in the genomic region 3, we downloaded all the 120 genes located in the region from the Soybase website (Table S8). Among them, 16 genes are located in qrxp_17_2. Haplotype analysis showed that the SNPs in this genomic region were clustered into five haplotype blocks ( Figure 6). The stable QTNs associated with BLP resistance were clustered into the second haplotype block covering 21 genes located in the 170-kb-long genomic region from 7.59 to 7.76 Mb (Figure 6). The remaining 83 genes are located between qrxp_17_2 and the second LD block (Table S8). Genes in the literature that meet one of the following criteria were predicted to be candidate BLP resistance genes: having been identified by gene mapping in soybean populations under Xag treatment, being responsive to the infection of Xag strains in soybean, and the paralogous genes responsive to the infection of Xanthornonas strains in other plant species. As a result, a total of four genes (i.e., Glyma.17G086300, Glyma.17G090100, Glyma.17G090200, and Glyma.17G090400) located between qrxp_17_2 and the second LD block were identified to meet at least one of the three criteria (Table 4). Thus, these four genes together with Glyma.05G040500 (the homolog of Glyma.17G086300, located near qrxp_5_1) were selected as the candidate BLP resistance genes for further analysis (Table 4).

Differential Expression Analysis of Candidate BLP Resistance Genes in Soybean
In order to quantify the relative gene expression levels of each of the five candidate resistance genes in resistant and susceptible soybean cultivars following Xag (strain EB08) inoculation, we performed real-time RT-PCR (qPCR) in the leaf samples of cv. W82 (Xag resistant) and cv. Jack (Xag susceptible) at 0, 12, 24, 48, 72, and 120 hpi by using Cons4 and Cons6 as the internal control (or reference) genes (Tables S9 and S10) [25]. The average relative expression levels of Glyma.17G090100, Glyma.17G090200, and Glyma.17G090400 in W82 showed no significant difference from that of Jack across the six time points after Xag EB08 treatments ( Figure 7A-C). In comparison to Jack, however, the average relative expression levels of Glyma.17G086300 was significantly increased in W82 at 24 and 72 hpi, while that of Glyma.05G040500 was significantly decreased in W82 at 24, 48, and 120 hpi ( Figure 7D,E). , and Glyma.05G040500 (E) in BLP-resistant (W82) and -susceptible (Jack) soybean cultivars following Xag infection. The Xag strain EB08 was resuspended in 10 mM MgCl 2 buffer to a final concentration of 1 × 10 8 CFU/mL. Then, the bacterial suspension was injected on the leaves of W82 and Jack with 10 mM MgCl 2 buffer being used as the mock treatment. The inoculated soybean leaves were sampled at 0, 12, 24, 48, 72, and 120 hpi, representing different stages of Xag infection. The relative expression levels were measured using qRT-PCR with the ATP-binding cassette transporter gene Cons4 and an F-box protein family gene Cons6 being used as the reference genes. The relative expression levels were calculated using the comparative threshold Pfaffl method. Values are the mean ± the standard error of three independent biological repetitions. Asterisks denote statistically significant difference (* p ≤ 0.05, ** p ≤ 0.01) between W82 and Jack by Student's t-test.

Discussion
Soybean BLP is a worldwide disease and causes serious soybean yield reduction under high temperature and humidity. To date, only the soybean genotypes containing the unidentified recessive resistance gene rxp have been used in BLP resistance breeding, even though soybean BLP resistance is a quantitative trait that is regulated by multiple genes and affected by the interaction between genotypes and environment in the field. Although several BLP resistance QTLs have been identified previously, little is known about the underlying molecular mechanisms of BLP resistance in soybean. In the present study, we combined linkage mapping and GWAS to identify the BLP resistance QTLs and QTNs in soybean, leading to the identification of a potential candidate gene for rxp.
The linkage mapping and GWAS methods are frequently utilized to identify disease resistance QTLs. Linkage mapping is based on the co-segregation of genetic regions over the genome in bi-parental families [26,27]. A drawback of this method is that the construction of the RIL populations is very time-consuming and laborious. Another drawback of this method is that it is only able to detect allelic diversity between the bi-parents, and its resolution highly depends on the number of recombination events. Considering QTL intervals usually extend over several cM, this limited resolution makes it extremely challenging to identify the target genes. This challenge could be overcome by using GWAS, which studies nature populations composed of numerous varieties [28]. GWAS can be used to detect the loci with multiple alleles over a genome. It can also be used to narrow down the detected QTL regions [28]. This is why GWAS has been used to identify the target genes associated with traits in various plant species [23,29,30]. However, the effectiveness of GWAS was limited by factors such as the false positives, linkage disequilibrium, complex population structures, and the alleles with the substantial effects on trait phenotypes but low frequencies in the populations [31,32]. Thus, linkage mapping and association analysis can be used together to facilitate genetic architecture analysis and the identification of target genes underlying large QTLs [33,34].
Earlier studies have reported four QTLs/marker-trait association (MTA) genomic regions on Chr 01, 10, 11, and 17 in soybean, correspondingly, with the single recessive resistance gene rxp being located on Chr 17 [2,17,19]. In this study, we identified three genomic regions that contain three BLP resistance QTLs and 11 stable BLP resistance QTNs and were significantly associated with BLP resistance by combining linkage mapping and GWAS analysis (Figure 3; Tables 2 and 3). Moreover, two new minor QTLs, qrxp_5_1 (in genomic region 1) and qrxp_17_1 (in genomic region 2), were identified on Chr 05 and 17, respectively, for the first time (Table 2). These two new QTLs were further confirmed by the co-localization of QTNs Gm05_7667820 and Gm05_7668047 with qrxp_5_1 and QTNs Gm17_5628119 and Gm17_5628133 with qrxp_17_1 (Tables 2 and 3). Moreover, genomic region 3 containing qrxp_17_2 and eight QTNs was located in the same genomic region as rxp (7.30 Mb) (Figure 4). We also evaluated the effects of the most significant SNPs Gm05_7667820, Gm17_5628119, and Gm17_7603802, which were distributed in the three genomic regions of soybean, on the improvement of BLP resistance. The average resistance of soybean lines containing all the three resistance alleles significantly increased when compared with that containing the three susceptible counterparts ( Figure 5). Thus, BLP-resistant soybean cultivars can be developed through polymerization breeding. Furthermore, the QTLs and QTNs identified in this study can be used for marker-assisted selection of BLP-resistant soybean lines.
Instead, we identified four candidate resistance genes for rxp in the present study. We further tested the relative expression levels of these four candidate genes together with Glyma.05G040500 (the homolog of Glyma.17G086300, located near qrxp_5_1) in W82 (resistant) and Jack (susceptible) at 0, 12, 24, 48, 72, and 120 hpi under the infection of virulent Xag strain EB08. In order to overcome the dilution effect, we selected eight locations along the veins of each leaf for pathogen injection. In addition, we designed six sampling time points after treatment to prevent false negative results caused by insufficient sampling time points. According to the qPCR results, no significant difference was found between W82 and Jack at 12 hpi, which is consistent with the result of the reported RNA-Seq analysis [18]. However, we found that the expression of Glyma.17G090400 increased 10 times in W82 and Jack under infection than mock treatment at 72 hpi. Interestingly, the relative expression levels of two Lob genes Glyma.17G086300 and Glyma.05G040500 showed a significant difference between W82 and Jack treated with EB08. When compared with that in Jack, the former significantly increased in relative expression levels in W82 at 24 and 72 hpi, while the latter significantly decreased in relative expression levels in W82 at 24, 48, and 120 hpi (Figure 7). Glyma.17G086300 and its paralog Glyma.05G040500 were predicted to be targeted by virulence factor tal2b, a transcription activator-like effector (TALE) of Xag [20]. Xanthomonas pathogens can secrete TALEs to activate multiple susceptibility genes in many plant spaces. For example, it can secrete TALE PthA to activate the expression of the canker susceptibility gene CsLOB1 in citrus [21]. Taken together, Glyma.17G086300 is a candidate gene for rxp.

Plant Materials
A soybean RIL population and an association panel of soybean breeding lines were used in this study. The RIL population consisting of 126 F 2:9 lines was developed by using the soybean cvs. Meng8206 and Zhengyang, which are resistant and susceptible to Xag strain C5, respectively [39], as the parents, and used for linkage mapping. These RIL lines along with the two parents were planted at Jiangpu Experimental Station, Nanjing, Jiangsu Province (latitude 33 Seeds were sowed in a randomized complete blocks design (RCBD) with three replications with 1 × 0.25 m hill plots and 50 cm row spacing. Five soybean lines with five biological replicates per line were planted in each single row plot with 20 cm distance between accessions. Three seedlings of each line were kept and grown to maturity after seed germination. Field management was conducted under normal conditions.

Pathogen Inoculation
BLP resistance was assessed in these soybean plants under both artificial inoculation and natural morbidity conditions. The artificial inoculation was carried out by using a highly pathogenic Xag strain C5 [1] to evaluate BLP resistance in the soybean RIL and association panel lines at Jiangpu in 2014, 2015, and 2016, which were named as '2014JPspray', '2015JP-spray', and '2016JP-spray', respectively. The C5 strain was restreaked on nutrient agar (NA) medium (BBL; Becton Diskinson and Co., Cockeysville, MD, USA) from a 30% glycerol stock at −80 • C and incubated at 28 • C for 24-48 h [1]. Then, a single colony was grown overnight in the nutrient broth (NB) liquid medium (BBL; Becton Diskinson and Co. Ltd., Cockeysville, MD, USA) on a rotary shaker (220 rpm) at 28 • C. The bacterial culture was diluted to a final concentration of 3 × 10 8 CFU/mL using sterilized distilled water and sprayed onto both sides of the soybean leaves at the early flowering stage with a high-pressure atomizer [40]. Ten days later, the inoculated plants were sprayed again with the same strain.
The natural morbidity condition in 2018 was used to evaluate BLP resistance in the association panel at Jiangpu and Dangtu Experimental Stations due to an outbreak of BLP disease in the soybean growing area of Jianghuai River Basin in China in August 2018, which resulted from the stormy weather brought by the tropical storm 'Rumbia'. These inoculation environments were named as '2018JP-natural' and '2018DT-natural', respectively.

Disease Assessment
To evaluate the resistance of the soybean RIL and association panel lines under the artificial inoculation condition, we counted the disease spot numbers on the inoculated leaves of soybean plants at 14 days post-inoculation (DPI). The infection responses (IR) in the soybean genotypes under the artificial inoculation condition were grouped at five levels: Level 1, no apparent disease spots; Level 3, the existence of 5-20 localized disease spots; Level 5, the existence of 20-50 irregular spots; Level 7, the existence of >50 irregular disease spots that formed a small lesion covering 10-25% of leaf area; and Level 9, a large lesion covering more than 25% of leaf area (Figure 1).
When we considered the severity of the BLP disease under the natural morbidity condition in 2018, the five levels used to evaluate the susceptibility of soybean lines under the natural morbidity condition were Level 1, 0-20 disease spots; Level 3, 20-50 irregular disease spots; Level 5, >50 irregular disease spots that form a small lesion covering 10-25% of leaf area; Level 7, a large lesion covering 25-50% of leaf area; and Level 9, a large lesion covering more than 50% of leaf area (Figure 1).
Three biological replicates of each line were planted in a randomized complete block design with three blocks per environment. The average resistance levels of the three replicates per block in each environment were used for linkage mapping and association analysis.

Phenotypic Data Analysis
Descriptive statistics, such as mean, standard deviation (SD), maximum and minimum trait values, coefficients of variation (CV%), and skewness and kurtosis, in both soybean populations for BLP disease response were calculated using the origin pro 2018 software (Origin Lab, Northampton, MA, USA). The analysis of variance (ANOVA) for BLP disease resistance level of both soybean populations was performed to evaluate the effects of genotype (G), environment (E), and genotype × environment interaction (G × E) in the joint environments using the PROC GLM of SAS 9 (SAS Institute, 2010, Inc., Cary, NC, USA). Broad-sense heritability (h 2 ) of BLP resistance of both soybean populations for combined environments was estimated as the following equation: h 2 = σ 2 g/ (σ 2 g + σ 2 ge/ n + σ 2 e/ nr), where σ 2 g is the genetic variance, σ 2 ge is the genotype × environment interaction variance, σ 2 e is the error variance, n is the number of environments, and r valid points is the number of replications within an environment [41].

Genetic Linkage Map Construction and Linkage Mapping in the Soybean RILs
RAD-seq (restriction-site-association DNA sequencing) [42] was used to genotype the 126 individuals of the RIL population and the two parents. Briefly, the genomic DNA was extracted from approximately 1 g of young leaves of individual plants of the RIL population and the two parents using the cetyltrimethylammonium bromide (CTAB) method [43]. The DNA fragments of 400~700 bp in length obtained by TaqI digestion were sequenced on Illumina HiSeq 2000 instrument following the standard protocol for MSG (multiplexed shotgun genotyping), and 90 mer paired-end reads were generated. The reads from RIL lines were aligned against the Williams 82 reference genome (Glyma.Wm82.a1.v1.1) [44] using the SOAP2 software (version 2.21) [45]. The SNP calling was performed using Real SFS software [46] on the basis of the Bayesian estimation. Finally, the high confidence SNPs were obtained by following the filtering criteria as follows: 50 < depth < 2500, a probability of site mutation = 95%, and every SNP loci being separated by at least 5 bp. Bin maps were then constructed for the 126 RIL lines. Consecutive SNPs were examined with a slightly modified sliding-window approach (window size: 15 SNPs, step size: one SNP) [47]. Recombination breakpoints were determined as the window sledded along the chromosome. Windows with 11 or more SNPs from either parent and windows with fewer SNPs from a single parent were considered to be homozygous and heterozygous, respectively. The consecutive 30-kb-long intervals with the same genotype in each line were joined into a bin using a PERL script according to the breakpoint information [48]. Finally, the genetic map of bin markers was constructed for the RIL population using JoinMap 4.0 [49].
Composite interval mapping (CIM) model of WINQTLCartographer2.5 [50] was used to detect QTLs for BLP resistance in RILs in each environment. A bin map with 2600 bin markers was constructed, covering 2630.2 cM. The phenotypic variance explained by a single QTL was calculated by using WINQTLCartographer2.5 software with the PVE = (VG/VP) × 100% (PVE: phenotypic variation explanation; VG: genetic variance of QTL; VP: phenotypic variance). The walk speed was set at 1 cM, and the window size was set at 10 cM. A log likelihood of 2.5 was used as the threshold for the presence of QTLs [50]. The QTLs were named according to the normal nomenclature [49].

Genotyping, Population Structure, and LD Analysis in the Soybean Association Panel
Genotyping was conducted as described in Li et al. [26] and Karikari et al. [24]. Briefly, the genomic DNA of each plant in the association panel was extracted using the CTAB method [43] and genotyped by the restriction site-associated DNA sequencing (RAD-Seq) technology. After Taq I digestion, the DNA fragments of 400-600 bp in length were obtained and then sequenced using Illumina HiSeq 2000 instrument (Illumina, San Diego, CA, USA). After being filtered, the high-quality SNPs with a >5% minor allele frequency (MAF) were used for principal component analysis (PCA), kinship analysis, LD analysis, and GWAS analysis. The genotyping data are available in the NCBI database (PRJNA648781, https://www.ncbi.nlm.nih.gov/bioproject/?term=prjna648781, accessed on 20 May 2020) and the National Center for Soybean Improvement website (http://ncsi.njau.edu.cn/info/ 1149/2070.htm, accessed on 20 May 2020).
Population structure analysis including the PCA of whole-genome SNPs and the construction of the neighbor-joining tree for the 476 lines in the association panel was performed using the TASSEL software version 5.2 [51] as described in Karikari et al. [24].
The linkage disequilibrium (LD) between pairwise SNPs was analyzed by calculating the LD parameter value, i.e., the squared Pearson correlation coefficient (r 2 ) using the RTM-GWAS V1.1 software [52]. A vcf format genotype file including the 60,315 high-quality filtered SNPs created by the TASSEL software was used as the input file in the RTM-GWAS V1.1 software. For LD analysis, the maximum inter-locus distance was set at 5 Mb, and the minimum r 2 threshold was set at 0.05. The calculated r 2 against the physical distance between pairs of SNP markers was visualized using the origin pro 2018 software (Origin Lab, Northampton, MA, USA) with scatter plot and polynomial fitting. The LD decay rate was the physical distance where the r 2 dropped from its maximum value to the half [51].

Association Analysis and Haplotype Block Analysis
Association analysis was performed by running mixed linear model (MLM) including the principal component analysis matrix (PCA) and the kinship matrix (K), i.e., MLM (PCA+K) using the TASSEL software version 5.2 [51]. Manhattan and quantile-quantile plots (QQ-plot) were created by using the qqman R package to visualize the GWAS result [53]. The distribution of p-values and the significant loci associated with BLP resistance over the whole genome in the GWAS panel were shown in the Manhattan graph. The significant threshold −log10(p) = 4.0 was used to control the genome-wide type I error rate [53]. The SNPs detected in at least two environments were considered as being relatively stable. A haplotype block analysis was performed using the Haploview software (version 4.2) [54].

Gene Annotation and Candidate Gene Prediction
The names and physical locations of the identified genes were downloaded from the Soybase database (https://soybase.org, accessed on 20 May 2020). Gene annotations were conducted according to the NCBI database. All the genes that respond to Xag infection were considered as potential candidate genes.

qPCR
Total RNA extracted from leaf samples of soybean cvs. W82 (resistant) and Jack (susceptible) following the treatment of Xag strain EB08 at 0, 12, 24, 48, 72, and 120 hpi in Zhao et al. [25] was used for cDNA synthesis and qPCR analysis of the relative expression levels of the five candidate genes, i.e., Glyma.17G090100, Glyma.17G090200, Glyma.17G090400, Glyma.17G086300, and Glyma.05G040500 in both cultivars. Three biological replicates were designed for each treatment, and three leaves of each biological replicate were infected as three technical replicates. Primer design and stepwise optimization of qPCR conditions including annealing temperatures, primer concentrations, and cDNA concentration ranges were conducted as described in Zhao et al. [25] so that R 2 ≥ 0.99 and E = 100 ± 5% were achieved for the best primer pair for each gene. The relative expression levels were calculated using Cons4 and Cons6 as the two reference genes [25] and the Pfaffl method [55].