Williams 82 (W82) and Huaxia 3 (HX3) were provided by South China Agricultural University. The 315 soybean accessions used for this study were obtained from soybean research institutions in China ( Supplementary Table 1 ). Accessions were planted in December in field sites located in Guangzhou and Hainan, and these were used for disease evaluation and genetic analysis.

In December 2017, accessions were sown in the field of Hainan Experimental Station, China. Disease resistance was evaluated in February 2018, and the data were recorded as Y2018. In December 2018 and December 2019, accessions were grown in the field at South China Agricultural University, Guangzhou. PMD resistance was evaluated in March 2019 and March 2020, and the data were recorded as Y2019 and Y2020, respectively. All plants were scored and counted for PMD resistance ( Supplementary Table 1 ). All plants were challenged with M. diffusa as described previously (Jiang et al., 2019). Each plant was inoculated with spores of M. diffusa at stage V1 by brushing with PMD-infected leaves of susceptible plants maintained in the greenhouse, as described by Kang and Mian (2010). Three weeks after inoculation, PMD incidence was evaluated for the 315 soybean accessions. Accessions were scored for PMD resistance using a scale. A plant with no PMD colonies on any leaf was regarded as a resistant (R) line and scored “L0,” while a plant with one or more PMD colonies on its leaves was rated as a susceptible (S) line, depending on disease severity, from L1 grade to L5 ( Supplementary Table 1 ) (Li et al., 2016).

The sodium dodecyl sulfonate (SDS) method was used to extract DNA from GWAS populations from young unfolded trifoliolate leaves at stage V2 (second trifoliolate) (Perry, 2004). A minimum of 30 ng/µl of DNA, with an OD260/280 of 1.8–2.0, was used for resequencing, and greater than 2 µg of DNA was used to generate small fragment libraries for pair-end sequencing.

A total of 315 soybean accessions were selected for whole-genome re-sequencing using the Illumina HiSeq X Ten platform with 150-bp pair-end reads. The average sequencing depth of all samples is 15×. A detailed sequencing information is shown in Supplementary Table 2 . All the sequencing in this study was done by Annoroad Gene Technology Co. Ltd. (Zhejiang, China).

Version 4 of Williams 82 genome in the Phytozome database (https://phytozome-next.jgi.doe.gov/info/Gmax_Wm82_a4_v1) was used as the reference for data analysis. Read quality control to remove low-quality reads and adapter sequences was done using Fastp (version 0.23.2) (Chen et al., 2018) with default parameters. Clean reads from each sample were aligned to the reference genome using BWA-MEM 2 (release 2.2.1) (Vasimuddin et al., 2019), and then, mapped reads were sorted using SAMtools (version 1.15.1) (Li et al., 2009). GATK4 packages (release 4.2.6.1) (Mckenna et al., 2010) were employed for variation detection and haplotype analysis. Duplicate-read marking was performed using MarkDuplicates. Variants from each sample were called using HaplotypeCaller to generate files in gVCF format. The gVCF files were merged using CombineGVCFs generating a VCF file composed of emit-all-sites information. The VCF file was filtered using the recommended GATK parameters “QUAL < 30.0 || QD < 2.0 || MQ < 40.0 || FS > 60.0 || SOR > 3.0 || MQRankSum < −12.5 || ReadPosRankSum < −8.0” and “QUAL < 30.0 || QD < 2.0 || FS > 200.0 || SOR > 10.0 || ReadPosRankSum < -20.0 || MQ < 40.0 || MQRankSum < −12.5” for SNPs and InDels, respectively. Apart from this, SNPs were filtered using the following criteria: removing SNPs as missing greater than 15% and MAF lower than 0.5% among 315 accessions.

TASSEL (version 5.2.8.1) (Bradbury et al., 2007) was used to calculate the association of the PMD resistance phenotype with genetic polymorphisms in 315 soybean accessions by applying a mixed linear model (MLM). In this study, admixture software was employed to estimate the population structure matrix (Q). When K = 11, the CV error was the minimum value at 0.54741. TASSEL5 software was used to calculate the kinship matrix (K) with the parameter-method Centered_IBS. As covariates, Q and K were used to control the population in the MLM analysis. As part of the TASSEL analysis, the phenotypic variation explained (PVE) was also calculated using the mixed liner model. The CMplot package (Yin et al., 2021) was used for visualizing the Manhattan p and QQ plots. The Manhattan plot shows the p-values of significant SNPs associated with PMD resistance. The QQ plot shows the relationship between theory and expectation. SnpEff (Cingolani et al., 2012) was used for gene annotation of heterotopic points according to the physical location and structure of genes throughout the genome.

LD analysis measures the nonrandom association of pairs of SNPs. TASSEL (version 5.2.8.1) (Bradbury et al., 2007) software calculates the square (R2) of the allele frequency correlation of SNP. The core SNP set derived after suitable filtering was used for calculating linkage disequilibrium (LD) patterns and the structure of haplotype blocks of SNPs by LDBlockShow (Dong et al., 2020). Data sets (r2 and D′) of pairwise LD measurement calculations showed LD levels of SNPs.

The genotype code types and physical position in the genome of the candidate locus were obtained from genotype-calling VCF format files. Unsupervised clustering was applied to three genotype codes consisting of reference homozygote, alternative homozygote, and missing by Cluster3.0 software (Hoon et al., 2004).

V1 stage W82 and HX3 plants were infected with an M. diffusa spore suspension containing 1 × 10⁵ cfu/ml and maintained in a growth chamber at 23°C, 75% relative humidity, with a 16-h light/8-h dark photoperiod to characterize the expression of 21 candidate genes. The experiments included three replicates. Leaves were sampled at 0, 6, 12, 24, 48, and 72 h after inoculation and kept at −80°C. Total RNA was extracted using the HiPure Plant RNA Mini Kit (Magen, Guangzhou, Guangdong, China), and 1 mg of total RNA was reverse transcribed to produce first-strand cDNA using the Evo M-MLV RT Kit with gDNA Clean for qRT-PCR (Accurate Biology, Guangzhou, Guangdong, China). Candidate genes in the target region were predicted in SoyBase using the Wm82.a4.v1 reference genome. Quantitative real-time polymerase chain reaction (qRT-PCR) was performed to obtain the expression profiles of candidate genes using primers designed with Primer Premier 5.0 (Premier, Vancouver, Canada). The housekeeping gene Actin was used as a control. The specific primers for 21 candidate genes are listed in Supplementary Table 6 . The qRT-PCR was performed with a CFX96 Real-Time PCR Detection System (Bio-Rad, Hercules, CA, USA) using SYBR^® Green Premix Pro Taq HS qPCR Kit II (Accurate Biology, Guangzhou, Guangdong, China). All reactions were performed in 20-µl volumes containing 1 µl of cDNA as a template. The thermal cycling conditions consisted of 94°C for 2 min, followed by 40 cycles of 95°C for 10 s, 59°C for 30 s, and 72°C for 30 s. Three independent biological repeats were used. The qRT-PCR data were evaluated using the 2^−△△CT method (Livak and Schmittgen, 2001).

The information of mutant lines from the Williams 82 EMS-induced library was extracted from this website: http://isoybean.org/ ( Supplementary Table 7 ). Seeds for mutants were obtained from Song’s laboratory (Zhang et al., 2022). These mutant lines were grown in the greenhouse at South China Agricultural University, Guangzhou. The mutant lines were subjected to two generations of selfing to obtain homozygous lines of the corresponding candidate genes ( Supplementary Table 7 ). Each mutant plant was inoculated with spores of M. diffusa at stage V1. After 2 weeks of inoculation, the M. diffusa of the mutant leaves was observed and photographed, and these leaves were stored at −20°C for DNA extraction and sequencing.

Candidate gene protein sequence from the reference genome annotation Wm82.a4.v1 version was used to blast the RefSeq Select proteins database (https://www.ncbi.nlm.nih.gov/refseq/refseq_select/). The top 100 sequences with the highest similarity were selected for complete multiple sequence alignment. The phylogenetic tree was constructed using the NJ method in MEGA X (Kumar et al., 2018) with 1,000 bootstrap replicates.

To study PMD resistance in soybeans, we collected 315 soybean accessions from around the world and evaluated their resistance to PMD during 3 years in Hainan province (P.R. China) and an experimental field at South China Agricultural University (SCAU). We divided our screening results into five levels (L1–L5) of PMD susceptible according to the previous studies (Kang and Mian, 2010) and defined PMD resistant (level 0; L0). Representative phenotypes of highly PMD susceptible (L5), medium sensitivity (L3), and resistant (L0) are shown in Figure 1A , and the PMD resistance (L0) or susceptible level (L1–L5) of each accession is detailed in Supplementary Table 1 . During the spring of 2018, we identified 196 resistant and 97 susceptible lines at the Hainan Experimental Station ( Figure 1B ). During the spring of 2019, we identified 207 resistant and susceptible 98 lines at the experimental station at South China Agricultural University (SCAU) ( Figure 1B ). During the spring of 2020, we identified 205 resistant and 110 susceptible accessions at the SCAU experimental station( Figure 1B ). Pooling these results across years and locations, we found that 206 (66.35%) of the 315 accessions are resistant to PMD, and 109 (34.60%) are susceptible. These results indicate that soybean accessions are a rich resource for PMD resistance ( Supplementary Table 1 ).

We performed whole-genome re-sequencing of 315 soybean accessions from China, the United States, Brazil, Australia, and Africa ( Supplementary Table 2 ), with an average 15× depth of sequencing. We mapped 92.65% of the reads to the soybean reference genome (Wm82.a4.v1 version) and used GATK to call 4,900,642 high-quality SNPs and 702,316 InDels (<40 bp) at bi-allelic loci. A distribution analysis of SNPs and InDels across the 20 chromosomes of the soybean genomes shows a high density near the ends of chromosomes 3, 6, 15, 16, and 18 ( Supplementary Figure 1 ). Each chromosome has an average of 245,032 SNPs and 35,116 InDels.

We used the SNPs and InDels from the diversity collection of 315 accessions to conduct a GWAS with a mixed linear model approach to identify variant sites that are significantly associated with PMD resistance ( Figures 1C–F ). A total of 56 SNPs surpassed the genome-wide significance threshold clustered at the end of the long arm of chromosome 16 ( Figure 1G ). The SNP locus with the highest −log10 (p-value) of 19.14 is at nucleotide position 37,418,263 on Chr16 and accounts for 33.38% (R2) of the phenotypic variance. This GWAS peak is flanked by two paralogous R genes (Glyma.16G214200 and Glyma.16G214300). Two other genes (Glyma.16G213700 and Glyma.16G214641) overlap several SNPs in peaks with the second highest significance ( Figure 1G ). We defined 13 genes in this region as candidate PMD resistance genes ( Supplementary Table 3 ). Furthermore, this region encompasses 8 NLR genes potentially associated with resistance to PMD, thereby yielding a total of 21 candidate genes for subsequent comprehensive analysis ( Supplementary Table 3 ).

We used LD Blockshow to detect haplotype blocks in the 400-kb region containing the most significant SNP markers ( Figure 1G ). There is strong linkage disequilibrium across the region. The most significant SNP out of 60 significant SNPs in the block (located at Chr16: 37,416,562–37,418,263) ( Supplementary Figure 2 ) is upstream of Glyma.16G214200 and downstream of Glyma.16G214300 suggesting that PMD resistance may be associated with variants in DNA regulatory elements of Glyma.16g214200 or Glyma.16g214300.

We analyzed SNPs and InDels at the 21-gene (Glyma.16G213700–Glyma.16G215400) locus in the diversity pane of 315 soybean accessions. Unsupervised clustering was used to cluster all SNP variant types based on three genotype codes: reference homozygote, SNP/InDel mutation (Alter), and missing ( Figure 2 ). There is a mutation caused by a base substitution at position 746 bp in the coding region of the Glyma.16G214000, which is associated with 71.56% (78/109, out of 109 susceptible accessions, 78 contained the mutation) of susceptible accessions ( Figure 2 ). SNP mutations in the upstream promoter region of Glyma.16G214200 from −1,534 to −1,978 bp are associated with 95% (104/109) of susceptible accessions ( Figure 3 ). Furthermore, several NLR candidate genes with relatively high peak values are located between Glyma.16G214641 and Glyma.16G215400 suggesting their potential involvement in PMD resistance. Haplotype analysis of these genes identified an SNP at position 72 in the coding region of the Glyma.16G215100 gene, and this mutation was present in 55% (60/109) of susceptible accessions ( Figure 4 ). There is a 3-bp deletion at 226 bp in the coding region of Glyma.16G215300, and 91% (99/109) of susceptible accessions are associated with this mutation ( Figure 5 ). The haplotypes of other candidate genes are depicted in Supplementary Figures 3–17 . In summary, SNP/InDel mutations in the coding and promoter regions of the above genes are closely related to their disease resistance, and these genes may be candidate genes for resistance to PMD.

The GmRmd1 (Glyma.16G214300) gene has been previously confirmed as a resistant gene against PMD in various studies (Xian et al., 2022). In this study, the majority of accessions possessing intact GmRmd1 coding regions demonstrate PMD resistance, whereas those lacking such regions exhibit susceptibility to PMD ( Figure 6 ). These findings underscore the critical importance of maintaining an intact coding region for GmRmd1 in conferring soybean resistance against accessions. However, the integrity of the GmRmd1 coding region was not consistent with the PMD phenotype in 11 of the 315 accessions in which three of the resistant accessions exhibited incomplete GmRmd1; concurrently, eight of the susceptible accessions displayed intact GmRmd1 ( Supplementary Table 4 ). This suggests the presence of additional genetic factors contributing to resistance against PMD. Subsequently, through the analysis of haplotypes in other candidate genes (Glyma.16G213700–Glyma.16G215400) within these 11 accessions, and by integrating high-quality haplotype data ( Supplementary Tables 4 , 5 ), it was observed that among the eight susceptible accessions, haplotypes (D3 of Glyma.16G213700), (H18 and H19 of Glyma.16G213800), (H15 and H16 of Glyma.16G214000), (D18 of Glyma.16G214200), (H11 of Glyma.16G214500), (H8, H19, and D7 of Glyma.16G214529), (H9 of Glyma.16G214557), (H41 of Glyma.16G214585), (H37, H40, and D1 of Glyma.16G214613), (H32, H50, H51, and D1 of Glyma.16G214641), (H12 and H18 of Glyma.16G215100), and (H8 of Glyma.16G215300) were found to be susceptible in other accessions. Among the three disease-resistant accessions, haplotypes (H4 of Glyma.16G214529), (H29 and H36 of Glyma.16G214585), (H8 of Glyma.16G214613), (H6 of Glyma.16G215100), and (H2 of Glyma.16G215300) exhibited resistance to diseases in other accessions. In conclusion, it is plausible that Glyma.16G213700, Glyma.16G213800, Glyma.16G214000, Glyma.16G214200, Glyma.16G214500, Glyma.16G214557, Glyma.16G214585, Glyma.16G214613, Glyma.16G214641, Glyma.16G215100, and Glyma.16G215300 may represent additional genes exhibiting resistance against PMD.

Through structural analysis of the coding regions of candidate genes, it was observed that among 206 disease-resistant accessions, a total of 203 exhibited intact Glyma.16G214300; conversely, out of 109 susceptible accessions, only 101 displayed incomplete Glyma.16G214300 ( Figure 7 ). Furthermore, within the cohort of disease-resistant accessions, a significant majority (193 out of 206) possessed intact Glyma.16G214800; in contrast, among the susceptible group (consisting of 109 samples), only 82 showed incomplete Glyma.16G214800 ( Figure 7 ). Similarly, within the population of disease-resistant accessions (206 in total), a substantial proportion (185) demonstrated intact Glyma.16G215000; however, among the susceptible counterparts (comprising 109 samples), merely 68 exhibited incomplete Glyma.16G215000 ( Figure 7 ). This observation suggests a strong association between the integrity of the coding region for Glyma.16G214300, Glyma.16G214800, and Glyma.16G215000 and their disease resistance phenotypes, and the resistance gene Glyma.16G214300 has been confirmed against PMD (Xian et al., 2022) suggesting that genes Glyma.16G214800 and Glyma.16G215000 are highly likely to confer resistance to PMD as well. It is plausible that these genes form a complex network involved in conferring disease resistance by synergistically exerting protective effects against PMD.

To investigate the expression of 21 candidate genes in response to M. diffusa infection, qRT-PCR was utilized to examine the expression patterns of each gene between Williams 82 (W82, resistance to PMD) and Huaxia 3 (HX3, susceptible to PMD) ( Figure 8 ). Among these genes, the expression levels of Glyma.16G213900 and Glyma.16G214000 were predominantly observed in W82, while they were scarcely detected in HX3. In comparison to HX3, the expression levels of Glyma.16G214100, Glyma.16G214300, Glyma.16G214500, Glyma.16G214585, Glyma.16G214800, Glyma.16G215100, and Glyma.16G215300 were upregulated highly in W82. Conversely, the expression levels of Glyma.16G214529, Glyma.16G214557, Glyma.16G214641, and Glyma.16G215000 showed lower expression levels in W82 compared with HX3. The expression levels of Glyma.16G214613 and Glyma.16G214669 exhibited similar patterns between W82 and HX3. However, Glyma.16G213700, Glyma.16G213800, and Glyma.16G214700 were expressed in HX3 but not detected in W82 suggesting their potential involvement in negative regulation of PMD resistance.

Additionally, six genes (Glyma.16G213900, Glyma.16G214000, Glyma.16G214500, Glyma.16G214300, Glyma.16G214585, and Glyma.16G214800) were induced by M. diffusa in W82 but showed minimal or low expression in HX3, and the expression levels of Glyma.16G213900, Glyma.16G214000, and Glyma.16G214500 gradually increased in W82 after M. diffusa infection, reaching their peak after 48 h and subsequently decreasing. Correspondingly, the expression levels of the Glyma.16G214300, Glyma.16G214585, and Glyma.16G214800 genes peaked after 24 h followed by a gradual decline. These results indicated that the above six genes could positively regulate PMD resistance. In addition, five genes (Glyma.16G214200, Glyma.16G214669, Glyma.16G215000, Glyma.16G215100, and Glyma.16G215300) were induced to express in W82 after being infected by M. diffusa. Among them, the expression level of Glyma.16G214200 reached its peak at 24 h, while that of Glyma.16G215000 peaked at 48 h. Similarly, the highest expression levels for Glyma.16G214669, Glyma.16G215100, and Glyma.16G215300 were observed after 72 h. However, these five genes exhibited distinct expression patterns in HX3 at each time point, suggesting a potential coordinated role of these genes in the regulation of PMD, necessitating their collaboration with other resistance genes for effective control of M. diffusa infection. To summarize, a total of 14 genes may be involved in the regulation of PMD resistance, with three genes specifically implicated in negative regulation.

Mutant lines ranging from 3700 (Glyma.16G213700) to 5400 (Glyma.16G215400), derived from Williams 82 induced with EMS, were obtained from Song’s laboratory (Zhang et al., 2022), and the specific mutation sites and homozygosity of each line are depicted in Figures 9B–J and Supplementary Table 7 . Based on the phenotypic characteristics of mutant lines, it can be inferred that rmd1-1, rmd1-2, rmd1-3, 4000, 4200, 4800, 5000, 5100, and 5300 lines exhibit susceptibility to PMD, while the mutant lines of other candidate genes, such as 3800 and 5400, are resistant to PMD ( Figure 9A ). Among them, the gene Glyma.16G214300 in the mutant lines rmd-1, rmd-2, and rmd-3 has been identified as a known resistance gene against PMD (Xian et al., 2022). However, Glyma.16G214300 was not mutated in lines 3800, 4000, 4200, 4800, 5000, 5100, 5300, and 5800 ( Supplementary Figures 18 - 22 ). This suggests that apart from Glyma.16G214300, Glyma.16G214000, Glyma.16G214200, Glyma.16G214800, Glyma.16G215000, Glyma.16G215100, and Glyma.16G215300 also significantly contribute to PMD resistance.