The mung bean RIL population of 200 lines was developed through the single-seed descent method from a cross of D2945 (smaller seed size and lower HSW) and D4702 (larger seed size and higher HSW), two landraces that originated from the association population of 558 Chinese mung bean landraces reported previously (Han et al., 2022) and were obtained from the Hubei Provincial Medium-term Gene Bank for Crops (Wuhan, China).

Two hundred mung bean plants of the F₆ generation, along with two parental lines, were evaluated across four locations from 2021-2022: Ledong (18.47° N, 108.54° E) in Hainan Province during 2021 (2021LD), Ezhou (30.40° N, 114.89° E) in Hubei Province during 2022 (2022EZ), Gucheng (32.29° N, 111.52° E) in Hubei Province during 2022 (2022GC) and Wuhan (30.58° N, 114.03° E) in Hubei Province during 2022 (2022WH). In 2021LD, mung bean lines were sown in late October and harvested in early January of the following year, benefitting from the tropical conditions of Ledong, Hainan Island. For 2022EZ, sowing occurred in mid-April during the spring season, with harvest occurring in early July. In both 2022GC and 2022WH, planting occurred in late June, with harvests conducted in October. Each location included two replicates. Every plot comprised a single 2-meter row with 11 plants spaced 20 cm apart within the row and 30 cm between rows. The weights of 100 randomly selected healthy, mature, and dry seeds were recorded via an electronic balance, and the average of three technical replicates was considered. The HSW used in subsequent analyses was derived from the average of these replicates per location.

For both the parental lines and the 200RILs, genomic DNA was extracted from young leaves of three two-week-old plants via the cetyl trimethyl ammonium bromide (CTAB) method, as described by Chen and Ronald (1999). A minimum of 2 µg of genomic DNA from each accession was used to construct sequencing libraries following the manufacturer’s protocol (Annoroad, Beijing, China). The libraries, with an average insert size of approximately 350 bp, were sequenced on a DNBSEQ-T7 platform, which produced 150 bp paired-end reads. The raw sequence data were cleaned via SOAPnuke v2.0.5 (Chen et al., 2018) to eliminate adaptor contamination and low-quality reads.

After performing the above filtering steps, good-quality reads were aligned to the high-quality chromosome-level mung bean M5311 reference genome assembly (Li et al., 2022a). For generating genomic alignments, BWA v0.7.17 (Li and Durbin, 2009) was employed, and the resulting SAM files were converted, sorted and indexed to the BAM format SAMtools v1.9 (Li et al., 2009). For SNP calling, GATK v4.1.8 (McKenna et al., 2010; parameters: –filter-expression QD< 2.0 || MQ< 40.0 || FS > 60.0 || SOR > 3.0 || MQRankSum< -12.5 || ReadPosRankSum< -8.0) was employed, and stringent filtering was applied to the SNPs on the basis of three criteria: (1) genetic diversity between the two parents; (2) support from more than four reads; and (3) conformity to a Mendelian segregation ratio of 1:1, as indicated by a chi-square test (P< 0.001). To facilitate mapping procedures, a genetic map was constructed with Joinmap v4.0 software by applying the Kosambi mapping function (Stam, 1993). A logarithm of odds (LOD) threshold of 5.0 was set.

QTL analysis was performed via the R/qtl package (Broman et al., 2003). The composite interval mapping (CIM) approach was employed to identify QTLs associated with HSW. The threshold LOD scores for each location’s HSW were determined through 1000 permutations (P< 0.05). QTL confidence intervals were determined via the 1.5 LOD-drop method. Additionally, a linear QTL model was applied to estimate the additive effects of the QTLs and the proportion of phenotypic variation explained.

Markers were developed in the qHSW1 candidate interval via a penta-primer amplification refractory mutation system (PARMS) derived from whole-genome resequencing of the parental lines to further delineate and fine map the target genomic region. Each PARMS marker constituted a set of two forward primers and one reverse primer to efficiently determine genetic effects around the target locus. In addition, a FAM fluorophore tail sequence (5′-GAAGGTGACCAAGTTCATGCT-3′) or a HEX fluorophore tail sequence (5′-GAAGGTCGGAGTCAACGGATT-3′) was attached to the 5′ ends of both forward primers. These primers were synthesized by Gentides Biotech Co., Ltd. (Wuhan, China). The PCR plates were read via a TECAN Infinite M1000 plate reader, and SNP calling, as well as plot generation, was carried out via the online software SNpdecoder (http://www.snpway.com/snpdecoder/ accessed on 2 November 2022) with manual adjustments. PCR was performed following the protocol outlined by Li et al. (2022b). Detailed information about the PARMS markers is provided in Supplementary Table S2 .

Using the high-density genetic map, a RIL designated D120, which harbored a heterozygous fragment in the qHSW1 region within an otherwise homozygous background, was identified. To refine the mapping of qHSW1, plants carrying the heterozygous qHSW1 region, selected from the selfed progeny of D120, were further self-pollinated to establish a large F₂ population. Recombinant individuals from this F₂ population were identified via PARMS markers flanking qHSW1. Fine mapping was performed via a progeny-testing approach, as described by Han et al. (2020), in which the selected recombinant plants were selfed to produce homozygous recombinant plants (HRs) and nonrecombinant plants (HNRs) through marker-assisted selection. Significant phenotypic differences (P< 0.01) between HRs and HNRs indicated that the recombinant plant carried qHSW1 within its heterozygous region. Through progeny performance evaluation of all the recombinants, the qHSW1 locus was further delineated to a further narrowed interval.

Descriptive statistics and analysis of variance (ANOVA) for the HSW were estimated across the four locations via IBM SPSS Statistics v29.0.0. The broad-sense heritability (h² ) was determined according to the method described earlier by Hallauer and Miranda (1998), as described below:

where σ ɡ 2 represents the genetic variance, σ ɡ e 2 represents the interaction of genotype with the environment, σ ϵ 2 denotes the error variance, n represents the number of environments with replications, and r represents the number of replications per environment. The values for   σ ɡ 2 , σ ge 2 , and σ ϵ 2 were obtained from variance components via a generalized linear model (GLM). Correlation analysis was performed via Origin 2021 software. Differences between groups were evaluated with a two-tailed Student’s t test.

The seed morphology and HSW in all four environments (2021LD, 2022EZ, 2022GC, and 2022WH) differed significantly between the two parental mung bean accessions D2945 and D4702 ( Figures 1A, B ). The HSW in the D2945×D4702 RIL population displayed a continuous and approximately normal distribution pattern ( Figures 1C-F , Table 1 ), indicating that HSW was a quantitatively inherited trait in this study. The overall broad-sense heritability for the HSW trait was 0.83 across the four environments. Significant correlation values ranging from 0.65 to 0.83 ( Table 2 ) were also observed. Additionally, genotype, environment, and their interaction effects were found to be significant in exhibiting the HSW trait within the RIL population ( Table 3 ). These findings imply that genetic variation predominantly governs the phenotypic diversity of HSW.

A total of 1,673,185,938 clean reads were obtained from whole-genome resequencing in this study. Among them, the numbers of reads for D2945, D4702, and RILs were 83,967,545, 59,874,725, and 1,529,343,668, respectively. The sequence depths for D2945 and D4702 were 26.16× and 18.65×, respectively, whereas the average sequence depth for the RILs was 2.38×. Detailed sequencing data for the two parental accessions and RILs used in this study can be found in Supplementary Table S1 .

SNP identification yielded 889,575 SNPs between the two parental accessions. After filtering, 3521 high-quality SNPs were utilized to construct a linkage map. Finally, a linkage map comprising 1194 SNPs randomly distributed across 11 linkage groups was developed via the Kosambi mapping function integrated in Joinmap v4.1 software ( Figures 2A, B ), with an overall length of 2,493.82 cM and an overall marker-to-marker distance of 2.09 cM ( Table 4 ).

Combining the genotypes and phenotypes of the RILs, QTL mapping was performed through the CIM function in R/qtl, with the LOD threshold determined as 3.25 through a permutation test consisting of 1,000 iterations (P< 0.05). A total of 5 QTLs for HSW were identified through single-environment QTL analysis; these QTLs were located on chromosomes 1, 2, 3, 7 and 11 ( Figure 3 , Table 5 ), with LOD values ranging from 3.41 to 17.66, and explained 2.46% to 26.15% of the phenotypic variation (R²). The positive alleles of these five QTLs were all from the high-HSW accession D4702. Among them, qHSW1, qHSW3 and qHSW11 were detected in all four environments, and qHSW2 and qHSW7 were environment specific. Moreover, qHSW1 and qHSW11 explained more than 10% of the observed phenotypic variation ( Table 5 ). These findings suggest that qHSW1 and qHSW11 are stable major QTLs. The qHSW1 has shown significant effects in this and previous studies, indicating its potential importance in trait variation. Thus, we choose the qHSW1 for fine mapping.

For fine mapping of HSW1, a RIL, named D120, which was heterozygous in the qHSW1 region with a homozygous background, was screened out. Twenty plants with a heterozygous region of qHSW1 that were selected from the selfing progeny of D120 were self-pollinated to develop a large F₂ population. All the seeds of the large F₂ population were planted in the field, which resulted in 3,715 plants. A total of 34 recombinants were identified from this population via the PARMS markers M1-1 and M1-2 flanking the qHSW1 region. These 34 recombinants were further classified into eight distinct crossover events via four newly developed PARMS markers (M1-3 to M1-6) within the M1-1 to M1-2 region. Eight recombinants were selected from the 34 recombinant plants to represent these eight distinct crossover events. Genotype identification and phenotypic difference analysis were subsequently conducted on homozygous recombinant plants and homozygous nonrecombinant plants selected from the selfing progeny of each of the eight recombinants. Progeny performance evaluation revealed significant differences in HSW between homozygous recombinant plants and nonrecombinant plants corresponding to the four recombinants (R1, R2, R7, R8) with the heterozygous qHSW1 genotype in the M1-3 to M1-5 region. However, no significant difference in HSW was detected between the homozygous recombinant plants and the homozygous nonrecombinant plants corresponding to the four recombinants (R3, R4, R5, R6) with the homozygous qHSW1 genotype in the M1-3 to M1-5 region. These findings indicated that qHSW1 is located within the approximately 506 kb (Chr1:53868960~54374894 bp in M5311 RefGen) chromosomal region between M1-3 and M1-5 ( Figure 4 ).

To estimate the genetic effect of qHSW1, a set of NILs, designated qHSW1 ^C2945 and qHSW1 ^C4702, were generated by selfing D120 with genotype selection. The investigation and analysis of important agronomic traits in qHSW1 ^C2945 and qHSW1 ^C4702 revealed significant differences in HSW between NILs. The HSW of qHSW1 ^C2945 was approximately 0.54 g lower than that of qHSW1 ^C4702, and in addition, the number of seeds per pod of qHSW1 ^C2945 was significantly (P< 0.05) greater than that of qHSW1 ^C4702 ( Table 6 ). No significant differences were observed in traits such as pod width, pod length, protein contentof seed, starch content of seed, flowering time, plant height, or branch number ( Table 6 ).

To evaluate the potential utility of qHSW1 in mung bean breeding, marker M1-3, which were successfully classified into two distinct groups, was used to exploit allelic variations present within the D2945×D4702-derived RIL population. For validation, a subset from the RIL population (approximately 94 lines) was randomly selected for genotyping with the M1-3 marker. Among these, 50 lines presented the same genotype as D2945, whereas 40 lines presented the same genotype as D4702 ( Figure 5A ). Association analyses revealed that the M1-3 genotype was significantly associated with HSW (P< 0.05) in the D2945×D4702 RILs, where the D4702 allele had predominantly higher HSW (mean 5.97 g) contributions than the D2945 allele (5.43 g) ( Figure 5B ). Thus, marker M1-3 can be effectively employed in marker-assisted selection programs for HSW in mung bean.

Many of these loci had minor effects and were mapped to wider genomic/physical intervals (Somta et al., 2022). Enhancing the mapping resolution necessitates enormous population sizes and a higher density pool of good-quality markers. Furthermore, maintaining adequate replications while minimizing phenotyping errors, along with a high-resolution genetic map, would increase the precision of QTL mapping procedures. Here, we employed 200 mung bean RILs and developed a genetic map based on high-quality SNP markers that were genotyped via a whole-genome resequencing strategy. Two major QTLs, qHSW1 and qHSW11, were consistently identified for HSW across four environments. Further literature analysis revealed that qHSW1 corresponds to the QTL Sd100w8.1.1 described by Isemura et al. (2012) and aligns with the QTL Sd100wt8.1 identified by Kajonphol et al. (2012). Additionally, it partially overlaps with the QTL qSDWT8.1 reported by Muktadir et al. (2014). Collectively, these findings indicate that qHSW1 is likely a significant contributor to the natural variation observed in HSW within mung bean germplasm.

Marker-assisted selection employs molecular markers that are intimately or tightly linked to the target traits, thus offering opportunities for the selection of plants possessing favorable alleles that could significantly impact the targeted traits (He et al., 2014). Marker-assisted selection has been extensively applied in multiple economic crops, such as rice (Gouda et al., 2020), maize (Xu et al., 2020) and wheat (Sun et al., 2020), but rarely in mung bean. Owing to the limited advancements in functional genomics research on mung bean, only a small number of closely linked, reliable and MAB-compatible markers related to phenotypic traits have been identified (Somta et al., 2022). Consequently, there are few reports of MAS in mung bean.

In conclusion, qHSW1 exhibited stable genetic effects across four different environments over multiple generations. The major effect of the HSW QTL qHSW1 identified in this study was able to explain ~21% of the phenotypic variation and varied in HSW by 0.54 g ( Tables 5 , 6 ); thus, introgression of this QTL is an efficient approach for mung bean HSW improvement. This QTL has no additional influence on pod width, the protein content of seeds, flowering time, plant height, or branch number; and represents easy-to-use and specific features for MAS. Through a progeny test strategy, we successfully narrowed qHSW1 down to a 506 kb region. The M1-3 marker was tightly linked to qHSW1, thus allowing deciphering of exact qHSW1 coordinates and exclusively introgressing the candidate locus to develop an ideotype mung bean accession. As a next step, further fine mapping of qHSW1 will be performed, and functional markers for the candidate gene for qHSW1 will be developed.