Introduction

Mungbean (Vigna radiata (L.) R. Wilczek) or green gram holds immense promise for global food security and sustainable agriculture1. Mungbean (2n = 22) is a self-pollinated crop having a short life cycle (60–70 days to maturity), and a small genome size (579 Mb)2. Mungbean grains contain nearly 20–30% protein and 60–70% carbohydrates along with numerous vitamins, and minerals3. When consumed with cereals, they provide an ideal balance of amino acids required by the human body, especially in Southeast Asian countries3. Mungbean is highly preferred for its protein because it is easily digestible and has no antinutritional components. It has detoxifying, hypolipidemic, hypoglycaemic, anti-cancerous, antihypertensive, and immunomodulation properties3,4. Despite its significant economic importance, research on mungbean is mostly underfunded Internationally, especially in the field of genomics5.

Globally, mungbean is cultivated in ~ 7.3 m ha and is predominantly grown in South, East, and Southeast Asia with an average productivity of 721 kg/ha. India and Myanmar occupy the top positions with a combined contribution of 60% followed by China (16%), and Indonesia (5%) to world production6. India alone accounts for more than 70% of world acreage (~ 5.6 m ha) and produces a meager ~ 3.68 m tones7. The low mungbean productivity in India is mainly due to the use of poor-quality seeds, non-availability of suitable ideotypes for different cropping systems, low harvest index, poor crop management practices, and the susceptibility of crops to various biotic and abiotic stresses8,9.

Grain yield is the primary focus of any breeding program, but its complex nature with high environmental interaction can only be dissected with effective phenotyping and genotyping. Yield of mungbean crop is determined by several parameters like number of pods per plant, number of seeds per pod, plant height, number of leaves, chlorophyll content, maturity duration, grain weight, resistance and tolerance to various biotic and abiotic stresses10. Comprehending the genetics and genomics of yield-related agronomic traits is essential for incorporating these traits into elite crop varieties.

Mapping traits using bi-parental population is a popular approach due to its power and ability to identify the linked QTLs (Quantitative Trait Loci) when allelic frequencies are low compared to Genome-Wide Association Studies (GWAS). The success of QTL mapping depends on factors such as population size, marker type, and marker density11,12. In mungbean, limited reports are available that associate markers with complex yield-contributing traits13,14,15. In addition, these studies have mainly used low-throughput and poorly dispersed markers like SSRs, which do not provide accurate information about the number and ___location of QTLs. The use of RILs is advantageous due to high recombination and were used to map various traits in mungbean10,16.

High-density markers like single nucleotide polymorphisms (SNPs) are the most abundant form of DNA variation which can be identified by genotyping-by-sequencing (GBS)17. The use of high-density linkage maps can enhance the precision of effect estimates for detected QTLs, provide accurate QTL positions, and facilitate the identification of candidate genes underlying complex traits. Recently GWAS has been gaining importance, but in comparison to QTL mapping, it may lack precision due to multiple testing, prone to false positives, and can only identify common variants that influence a trait, despite the advantage of mapping with high resolution18.

Traditional breeding methods have played a central role in advancing mungbean breeding; however, their efficiency can be enhanced by integrating molecular approaches. The identification of QTLs associated with various traits, along with their candidate genes can unravel the genomic basis of specific traits. The development of user-friendly markers within the QTLs regions holds great promise for expediting mungbean breeding efforts. This study aims to identify the QTLs governing various yield-contributing traits and to develop the InDel markers closely linked to the QTLs for their use in marker-assisted selection (MAS).

Results

Variation, heritability, and correlation analysis

The phenotypic data for the parental genotypes and RILs displayed a remarkable range of diversity across all traits (Table S1; Table S2). The number of leaves per plant ranged from 8.3 to 37.3 (mean = 20.3); plant height varied from 49 to 100.6 cm (mean = 67.59 cm); 100 seed weight spanned from 2036 to 4632 mg (mean = 3178.05 mg); number of pods per plant ranged between 2.8 and 80.3 (mean = 30.29) and grain yield varied from 3.8 to 43.6 g (mean = 15.77 g). The skewness (range: -0.9 to 1.9) and kurtosis (range: -0.6 to 1.2) indicated a nearly symmetrical distribution around the mean (Fig. 1).

Fig. 1
figure 1

Distribution of yield-related traits in the RIL population. (a) Number of leaves; (b) Plant height in cm; (c) SPAD value measured after 25 days of sowing; (d) 100 seed weight in mg; (e) Number of pods; (f) Total yield of three plants in g. Where blue line indicates mean.

Full size image

Both genotypes and seasons showed a highly significant effect on the traits as revealed by ANOVA (Table S3). It means, the genetic composition of RILs and environmental conditions significantly influenced the trait expression. Additionally, significant interaction was noted between genotype and season for NL, PH, and SPAD values. Although, replications were significant for NL and GY, yet their contribution to the total variation was minimal. The broad-sense heritability estimates ranged from 64% (chlorophyll content) to 98% (seed weight). A significant correlation was recorded between total yield and NL (0.45**), NP (0.86**), and SPAD value (25 DAS 0.26** & 45 DAS 0.21*). Significant correlation was also recorded between SPAD values and NL (0.17*), NP (0.24**), and SW (0.16*) (Figure S1).

Linkage map

The linkage map was generated using 1374 filtered SNPs with Kosambi function. SNP markers were found evenly distributed across 11 chromosomes, with an average density of one SNP per 250 kb. Chromosome 7 exhibited the highest number of markers (278), while chromosome 4 had the lowest (69). Chromosome 1 exhibited lowest marker density (2.9 SNPs/Mb) and chromosome 3 displayed the highest density (8.7 SNPs/Mb) (Table 1, Figure S2-S4).

Table 1 Chromosome-based distribution of 1374 polymorphic SNPs which are used for QTL mapping.
Full size table

QTL mapping for grain yield-related traits

QTL mapping was conducted for six yield-related traits namely, plant height (PH in cm), number of leaves per plant (NL), number of pods per plant (NP), chlorophyll content (as SPAD value at 25 and 45 days after sowing), seed weight per 100 seeds (SW in mg), and total grain yield (GY in g) using Multiple Interval Mapping (MIM) and Composite Interval Mapping (CIM) approaches. BLUPs were used for the partitioning of phenotypic variance into its genetic and environmental components. Based on the LOD permutations, a threshold of 3.0 was set for identifying a QTL. The MIM detected 5 QTLs across all the traits, while CIM could identify 18 QTLs, including 5 QTLs which were detected by MIM (Table 2; Figs. 2 and 3; Figure S5-S16). These QTLs were located across 9 chromosomes (except on chr 1 and 8; Fig. 4) having LOD scores ranging from 3.05 to 9, and PVE ranged from 9 to 24%. Chromosome 7 had the highest number of QTLs (6 No), followed by chromosome 11 (3 No), chromosomes 2 and 4 (2 each), and chromosomes 3, 5, 6, 9, and 10 (1 each) (Table 2). The additive effect (‘a’, Table 2) represents the contribution of a QTL allele to the trait variation. A higher absolute value indicates a stronger genetic influence, while a positive effect indicates that the allele from Pusa Baisakhi enhances the trait value, and a negative direction suggests the allele from PMR-1, contributes to the increase. Yield-related QTLs were inherited from both parents but the QTLs associated with higher seed weight, pod number, and yield predominantly carried favourable alleles from Pusa Baisakhi, while those linked to taller plant height, increased chlorophyll content, and greater number of leaves were inherited from PMR-1, aligning with their respective agronomic characteristics.

Fig. 2
figure 2

Composite interval mapping of Yield-related traits in Mungbean. (a) Number of leaves; (b) Plant height; (c) SPAD value; (d) Seed weight; (e) Number of pods; (f) Total grain yield.

Full size image
Fig. 3
figure 3

Effect Plot of important QTLs identified for Yield-related traits in mungbean. Where AA represents PMR-1 and BB represents Pusa Basakhi genotypes.

Full size image
Fig. 4
figure 4

Distribution of QTLs for yield-related traits on chromosomes of mungbean identified using composite interval mapping (CIM).

Full size image
Table 2 QTLs identified for yield-related traits in an RIL population using MIM and CIM.
Full size table

Gene ontology and candidate gene analysis

Genes located within the QTL intervals were subjected to gene ontology and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis. To identify the candidate genes, polymorphic variants detected by WGRS of the parental genotypes were annotated. WGRS produced nearly 220 million raw reads (151 bp), with an average depth of 34× coverage. Very high variants were recorded between the parents (1.2 million; density: 1 variant per 300 bp), and of these, 3% were Insertions/Deletions (InDels). Nearly 66,000 variants fell within the functional class, of which 53% were missense, 1% were nonsense, and 46% were silent mutations. Genes exhibiting non-synonymous mutations in the exons were identified which include premature start codon gain variants, conservative in-frame deletions and insertions, disruptive in-frame deletions and insertions, missense mutations, frameshifts, exon loss, premature stop codons, start loss, and start gain.

In the QTL region, 298 genes showed non-synonymous variations in the exons, and of these, 6 candidate genes each for PH and GY, 4 each for NL, SW, and NP, and 1 for SPAD value were identified (Table 3). Some of the important genes were identified as candidates like pentatricopeptide repeat-containing protein (NL), Inositol-tetrakisphosphate 1-kinase 1-like (PH), PYR-1 (PH), ferredoxin-like (SW) (Figure S17), Amino acid transporter AVT3B-like (SPAD), etc. (Figure S18).

Table 3 Candidate genes identified within mapped QTL regions with non-synonymous variations.
Full size table

Validation and development of indel markers

Arabidopsis is one of the most extensively studied model organisms in plant genetics, with a fully sequenced genome and a wealth of functional genomics data. Hence, well-characterized orthologs of Arabidopsis were identified and their pattern of expression was analyzed. For traits like plant height, two genes were identified through digital gene expression among the six identified candidate genes (Fig. 5). This digital gene expression can provide hypotheses for further experimental validation. For further use of these identified QTLs, InDel primers with more than 15 bp of insertions or deletions were developed (Table 4).

Fig. 5
figure 5

Digital gene expression of plant height candidate genes (a) Inositol-tetrakisphosphate 1-kinase 1-like (LOC106777318) and (b) Abscisic acid receptor PYR1 (LOC106776425).

Full size image
Table 4 QTL-Specific indel primers based on parental polymorphisms.
Full size table

Discussion

Despite its significant economic importance, mungbean is ignored at both national and international levels, in the field of genomics. Mungbean was one of the first legumes in which mapping for bruchid resistance using RFLP markers was reported in 199219. However, advances in mungbean genomics lag behind other legume crops like chickpea and soyabean20. Available literature shows mapping in mungbean for only a few traits like seed size, seed dormancy, bruchid, bean bug, MYMV, and powdery mildew resistance20. Detailed QTL identification using RIL population for yield contributing traits is still missing in mungbean.

Typically, the mean performance of RILs is expected to fall within the range defined by the parental means for quantitative traits. However, the mean values of RILs were recorded lower than the lower parent, despite the presence of transgressive segregants displaying higher performance than both parents. This could be because of highly contrasting nature of the parental lines having co-adapted gene complexes due to the selection21. In the process of RIL formation, it seems that the co-adapted gene complexes got disrupted upon crossing and recombination22. After a few generations of selfing, heterozygosity is lost, which might have disrupted the coadapted gene complexes in certain RILs, as also reported for Avena barbeta23.

Grain yield was found positively associated with the number of pods, number of leaves, and SPAD values; whereas, it was found not significantly affected by the seed weight and plant height. Pulses typically begin flowering after reaching a specific number of nodes, suggesting that plant height and leaf count indirectly affect pod number and seed yield24,25. Leaves and their chlorophyll content play crucial roles in photosynthesis, light interception, and also as a sign for disease and pest resistance26. More number of leaves per plant with higher chlorophyll content enhances photosynthetic activity, resulting in more assimilates for pod formation and grain filling. Previous studies also reported a positive correlation between single-leaf photosynthetic rates, biomass production, and yield in pulses like cowpea, soybean27,28, and mungbean14,29.

Number of leaves per plant (NL)

The NL in mungbean plays a pivotal role in determining the yield, as it directly influences the pod yield by increasing the number of nodes and photosynthate accumulation26. The number of nodes was mapped on LG3 and LG10 using InDel markers30; while number of primary branches was mapped using DArT markers on chr.7 and 11 in a RIL population with 26% and 27.7% PVE, respectively31. Recently, 3 SNPs were found associated with primary branches on chr. 4, 9, and 11 using GWAS5. This study identifeid 4 QTLs (qNL-7-1, qNL-7-2, qNL-9-1, and qNL-11-1) for leaf number. Both MIM and CIM detected a qNL-7-1 (9.9% PVE) located in 37 Mb region on chr. 7, as also reported previously31. Within the 4 QTL regions, 5 candidate genes were identified exhibiting non-synonymous mutations and having a role in leaf growth, regulation of flowering, and cellular development (Table 3). The LOC106777903 gene, a pentatricopeptide repeat-containing protein, functions by directly interacting with the protein GIGANTEA, which regulates various physiological processes like flowering time, leaf senescence, hypocotyl elongation, chlorophyll accumulation, etc. Moreover, GIGANTEA may negatively influence leaf number by promoting flowering32.

Plant height (PH)

Plant height influences the overall architecture of the mungbean plant. The varieties with appropriate plant height can optimize light interception and resource utilization, leading to higher yield30. Two major QTLs (qPH-11-1 and qPH-7-1) and one minor QTL (qPH-2-1) controlling plant height were identified on chromosome numbers 11, 7, and 2 respectively. The qPH-11-1 accounted for 21.5% PVE, with a negative additive effect (-2.3), suggesting its inheritance from PMR-1. Within the three QTL regions, 32 genes are identified with mutations, and of these, six genes [LOC106778804 (chr.2), LOC106756212 (chr.2), LOC106777391 (chr.11), LOC106777318 (chr.11), LOC106776425 (chr.11), LOC106777991 (chr.11)] were having role in growth, flower development, energy metabolism, and environmental adaptation.

The LOC106756212 protein contains a Calmodulin binding protein-like N-terminal ___domain, having a role in protein localization to ER exit site, cellulose biosynthesis, and salicylic acid metabolic process. The Inositol-tetrakisphosphate 1-kinase 1-like (LOC106777318) gene is involved in the biosynthesis of inositol phosphates (IPs), which act as secondary messengers and plays a crucial role in auxin-related processes33. The digital gene expression analysis showed that this gene was primarily expressed during pollen maturation and seed maturation stages (Fig. 5). Similarly, the PYR1 (LOC106776425) gene regulates plant responses to abscisic acid (ABA) and it is highly expressed in the second internode of stem and also in the imbibed seeds34. ABA at its basal level promotes plant growth, however, at elevated levels, it impairs plant growth35. When PYR1 binds to ABA, it inhibits ABI1 by forming a complex thereby preventing ABI1 from dephosphorylating (and thus inhibiting) the SnRK2 kinases. This activation of the SnRK2 kinases leads to the expression of ABA-responsive genes, promoting stress responses inhibiting shoot growth and increasing root growth34. LOC106777991 is involved in the synthesis of porphyrin and chlorophyll metabolism, which are crucial in determining plant height. Similarly, various reports of QTL identification for stem length on various chromosomes are known like six QTLs on six chromosomes (1, 2, 3, 6, 9, 10)36, three QTLs on chr. 4, 5, 8 (2.9–20% PVE)37, four QTLs on chr. 3, 4, 8, 11 (4.2–16.6% PVE)31, four QTLs on four chromosomes (1, 3, 7, 10)30, and 2 QTLs on two chromosomes (1, 7)5.

SPAD (soil-plant analysis development) values

The SPAD value is an indirect measure of chlorophyll content in the leaves. It holds considerable importance in determining plant growth, physiological status, and yield12. A consistent QTL (qSPAD25-7-1) on chr. 7 at 3 Mb position (15.8% PVE) was identified for both 25 and 45 DAS having a role in chlorophyll accumulation, which is inherited from PMR-1. Another major QTL (qSPAD25-5-1) was identified regulating SPAD value on chr. 5 (12.2% PVE). Similarly, a QTL on chr. 7 (5.5 Mb position, 12% PVE) was identified12, while another QTL on chr. 1 was identified through association mapping in mungbean5. In soybean, QTLs for SPAD value were also found located on chr. 738,39. Within these two QTL regions, 22 candidate genes were identified with mutations, and LOC106765948 was found coding for the amino acid transporter AVT3B-like having a role in the transport of amino acids, mediating nitrogen distribution, which is crucial for plant chlorophyll and growth40.

Seed weight (SW)

Seed weight is a critical trait which directly affects the amount of reserves present in the seeds to provide energy for early seedling growth and establishment. Seed weight is also associated with seed quality attributes and bigger seeds are often preferred by consumers and processors41. It is one of the most extensively mapped traits where researchers have identified a number of QTLs on different chromosomes42,43,44.

This study identified 3 QTLs for SW (qSW-10-1, qSW-11-1, qSW-6-1) on chromosomes 10, 11, and 6 with a PVE of 24.4%, 17.5%, and 10.8%, respectively. Similarly, 11 QTLs (PVE: 5.2–15.1%) have been identified for SW using RFLP markers45, while 5 QTLs were identified on chr.1, 2, 10 and 11 (PVE: 12.6%) using SSR markers15. SW is also mapped using an association mapping approach5,11,14. The most significant QTL (qSW-10-1) was identified from the parent Pusa Baisakhi at 20.6 Mb on chr. 10 which had 7 genes with variations between the parents. Among these only one gene, ferredoxin-like was found to have a role in the biological processes related to seed development. It is functional in chloroplasts of photosynthetic tissues, participates in carbon assimilation, plays a role in photosynthetic electron transport of PS1, and affects other FD-dependent genes46. The expression of its ortholog in Arabidopsis thaliana showed that it was highly expressed in the torpedo stage of embryo development. As seed formation progresses, the expression of this gene increases, reaching a peak around seed stage 6, after which its expression declines during the maturation phase. The involvement of ferredoxin in early and mid-stage seed formation seems linked to energy transfer and metabolic processes required for seed growth. It is also expressed in other plant tissues, like leaves and flowers, having a role in electron transport and redox reactions. However, its expression in siliques is low47.

In addition, three more candidate genes identified include, ATP-dependent zinc metalloprotease FTSH 4 (LOC106763742), L-type lectin-___domain containing receptor kinase (LOC106763733), and alpha-L-fucosidase 2-like (LOC106777853 respectively). Of these, LOC106763742 codes for a cytoplasmic-membrane protein that controls photosynthesis, chloroplast, and mitochondrial development. The ATPase and proteolytic activities were reported in the ATP-dependent zinc metalloprotease FTSH 4 protein which influences overall seed weight48. Similarly, L-type lectin-___domain-containing receptor kinases (LecRLKs) are plasma membrane-located proteins having a carbohydrate recognition ___domain (CRD), needed for the seed development49. Also, the LecRK-VIII.2 functions as a regulator of both seed size and number50. Similarly, an Alpha-L-fucosidase 2-like enzyme is likely to have its role in seed weight regulation by modulating the cell wall composition, especially xyloglucans51.

Number of pods per plant (NP)

A higher number of pods per plant generally results in an increased yield. This study identified 2 QTLs for NP on chr. 2 (qNP-2-1) and 4 (qNP-4-1), from the parent Pusa Baisakhi having substitution effects of 2.9 and 2.8, respectively. Both the QTLs collectively explained 30% PVE. Similarly, previous studies have also reported QTLs for NP on chromosomes 1, 3, and 915, and on chromosomes 2, 4, and 736, collectively explaining 6.7–12.7% PVE. The identified QTLs are having 4 candidate genes (MADS-box transcription factor, UDP-glycosyltransferase 92A1, ADP & ATP carrier protein 1, and E3 ubiquitin-protein ligase RDUF1). The MADS-box transcription factors regulate plant, seed, and pod development in soybean52. Other candidate genes need further confirmation for their role in the determination of NP in mungbean.

Grain yield per plant (GY)

The total yield of a genotype represents a complex quantitative trait influenced by various qualitative traits and a multitude of interconnected genetic, environmental, and management domains. This study identified 3 QTLs (qGY-7-1, qGY-3-1, and qGY-4-1) governing GY (30% collective PVE) on chr. 7, 3 and 1, respectively. Similarly, 11 QTLs are identified on various chromosomes (8.3–63% PVE) and qGY-4-1 was found linked with qNP-4-1, suggesting their pleiotropic effects on both GY and NP31. In the 2.7 Mb QTL region (qGY-7-1, qGY-3-1, and qGY-4-1), 67 genes are identified having mutations and are also part of several pathways like pollination, reproduction, embryo development, post-embryonic development, and fruit ripening. Of these, 6 candidate genes viz., endoglucanase 1 (LOC106757749), squamosa promoter-binding-like protein 13 A (SPL protein) (LOC106768860), ARF guanine-nucleotide exchange factor GNOM (LOC106757622), E3 ubiquitin-protein ligase RDUF1 (LOC106758074), glutamate–cysteine ligase chloroplastic-like (LOC106767912), and bidirectional sugar transporter (LOC106766939) seems having role in GY determination. SPL genes encode plant-specific TF, influencing phase transition from vegetative stage to flowering to fruit development. Under the influence of hormones, SPL genes may increase yield by efficient distribution of energy and nutrients53. Ubiquitin-ligase RDUF1 regulates yield by controlling the degradation of key regulatory proteins involved in stress responses, hormone signaling, and growth pathways54.

Conclusions

The knowledge of QTLs associated with yield and its component traits is crucial for the development of high-yielding mungbean varieties. This study revealed a significant influence on yield-related traits by both environmental and genetic factors. Grain yield was found positively influenced by parameters like number of pods per plant, number of leaves per plant, and chlorophyll content; while seed weight and plant height did not show any noticeable effect. Two QTLs with very high PVE (> 20%) for plant height (qPH-11-1) and seed weight (qSW-10-1) were identified as the ideal candidates for introgression into mungbean breeding programs. Chromosome 7 of mungbean harbored 5 novel QTLs for plant height, grain yield, chlorophyll content, and number of leaves. This makes chromosome 7 as the key chromosome for the development of markers and marker-assisted introgression of various yield-related traits in mungbean. QTLs for grain yield (qGY-7-1 at 54.43 Mb) and number of leaves (qNL-7-2 at 50.67 Mb) are found linked as evident from correlation analysis. The QTLs and the InDel markers identified could serve as a tools for improving the mungbean yield through molecular approaches. This study is the first to employ high-density SNP markers for QTL mapping of yield-related traits in mungbean, leading to the development of QTL-specific markers.

Materials and methods

Mapping population

A RIL population (F9:10; 166 No) was developed by crossing two contrasting genotypes (Pusa Baisakhi × PMR-1). Pusa Baisakhi - A released variety having high yield, high hundred seed weight, lower chlorophyll content with lesser height and number of leaves. PMR1 – An advanced breeding line is having low yield, less hundred seed weight, higher chlorophyll content, and taller plants with a greater number of leaves than Pusa Baisakhi. The RILs were advanced via single seed descent (SSD) method (F2 to F9:10) and evaluated for six yield contributing traits at the Indian Agricultural Research Institute (IARI), New Delhi (28.7°N, 77.1°E, 228 m AMSL). The experiment followed a randomized complete block design (RCBD) with two replications, adhering to the standard agricultural practices. Each genotype was sown in a 5.0 m row, maintaining a spacing of 20 cm between plants and 30 cm between rows.

Phenotypic evaluation

The parental genotypes and RILs were phenotyped for three consecutive seasons during 2022 (Sept-Nov) and 2023 (March-May) and 2023 (Sept-Nov) for six traits viz., plant height (PH, cm), number of leaves/plant (NL), number of pods/plant (NP), chlorophyll content as SPAD value at 25 and 45 days after sowing (DAS), 100 seed weight (SW, mg) and total grain yield (GY, g). The data from three plants per replicate (NL, PH, SPAD, and NP) were averaged, while 100 seed weight was measured thrice per replication and averaged. Each tri-foliate leaf was counted as a single leaf, the SPAD value was recorded between 09:00 to 11:00 AM on the matured third leaflet from the top of the plant. The analysis of variance (ANOVA) was measured using Agricole package in R55, the descriptive statistics and genetic parameters were estimated using Variability package56, combined best linear unbiased predictor values (BLUPs) were calculated using MetaR57 and correlation was measured using R base58,59. The studied traits were also tested for normality using the Shapiro-Wilks test.

Sequencing and variant calling

Genotyping by sequencing (GBS) was performed by isolating genomic DNA from the leaves of RILs using DNeasy Plant Mini Kit (Qiagen, Germany) (NCBI BioProject: PRJNA914895, https://www.ncbi.nlm.nih.gov/bioproject/?term=PRJNA914895). The quality and quantity of genomic DNA were assessed using Nanodrop (Thermo Fisher Scientific, Waltham, MA, USA). Libraries were prepared using TruSeq® DNA Nano LP Library Prep Kit (Illumina, San Diego, CA, USA) and were sequenced to produce raw 2 × 150-bp reads (Novaseq 6000, Illumina, San Diego, CA, USA). Trimmomatic software was used to process the raw reads by removing the adapters and low-quality bases60. Following trimming, ≥ 85 bp reads were retained having Q30 (base calling accuracy > 99.9%). Cleaned FASTQ files were then mapped to the reference genome of Vigna radiata (mungbean), assembly GCF_000741045.1 corresponding to cultivar VC1973 (https://www.ncbi.nlm.nih.gov/datasets/genome/GCF_000741045.1/)61. Variant calling was performed using the Unified Genotyper in the Genome Analysis Toolkit (GATK) version 3.6 for identification of SNPs from RILs; and whole-genome resequencing (WGRS) was employed to sequence the parents Pusa Baisakhi (NCBI BioProject: PRJNA1191361, https://www.ncbi.nlm.nih.gov/bioproject/?term=PRJNA1191361) and PMR-1 (NCBI BioProject: PRJNA1191361, https://www.ncbi.nlm.nih.gov/sra/PRJNA1191375).

High-density map construction and QTL mapping

A total of 38,931 SNPs were identified through GBS across 11 chromosomes of 166 RILs. The SNPs and individuals (taxa) were filtered based on a missing data threshold (> 30%). Then monomorphic, heterozygous (> 5%), and significantly distorted (p > 0.05) SNPs were removed using Tassel62. Between parental lines, 5494 filtered polymorphic SNPs were identified and those significantly deviating from the expected 1:1 ratio were eliminated based on a Chi-square test (p ≤ 0.05). Finally, the linkage map was constructed using 1374 SNPs (146 RIL sequence data) with 8% missing proportion and 0.46 as the minor allele frequency (MAF) (Table S4). The linkage map was generated using Kosambi function having fixed genomic physical positions.

QTL mapping was performed using both MIM and CIM through Rqtl package in R Studio63 and genetic distance between markers was converted into centiMorgan (cM). Combined BLUP values were generated for accurate mapping. Genotypic data were imputed using a hidden Markov model (HMM) (error probability: 0.001; step size: 0.5 cM; simulations: 64). MIM was conducted using the Haley-Knott regression method and penalties were calculated (p = 0.05). For CIM, up to 10 marker covariates were included (window size of 1.0 cM). To find the significant threshold for Logarithm of Odds (LOD) scores, 1000 permutations were conducted and a predefined LOD score threshold of 3.0 was used to identify the QTLs. Additionally, the confidence intervals were calculated using the LOD drop-off method with a decrease of 1.0 LOD unit on either side of the QTLs. PVE (%) was estimated using the following formula: \(\:PVE=1-{10}^{\left(\frac{-2\times\:LOD}{N}\right)}\). Where, n: sample size and LOD: Peak LOD value64. The QTLs were classified as major if they explained > 10% of phenotypic variance, while those with PVE < 10% were considered minor12.

Identification of candidate genes and indel marker development

The WGRS was used to sequence the parental genotypes (Pusa Baisakhi and PMR-1) and sequence information was used to identify the variants within the QTLs identified through mapping. The raw reads were trimmed, aligned and resulting variants were functionally annotated using SnpEff ver. 5.1d65. The annotated information was used to identify the candidate genes exhibiting variations between the parents and potentially affected expression. For the development of markers, WGRS data was aligned to the reference genome and variant calling was performed and those InDels that are monomorphic or of < 15 bp size were excluded. Then 200 bp sequence flanking each nearest InDel to the QTL (100 bp on either side) was used for primer designing via primer3 web tool66. Primers with a length of 18–24 bp, melting temperature (Tm) of 53–62 °C and GC content of 40–60% were idealised. While, primers forming secondary structures (e.g., hairpins) and primer-dimers were avoided.

Digital gene expression analysis of identified candidate genes

The digital gene expression was done to validate the candidate genes using well-characterized Arabidopsis orthologs, by analyzing their localized expression patterns through Expression Angler, an ePlant web browser (https://bar.utoronto.ca/eplant/)47,67. The Arabidopsis orthologs were retrieved using OrthoDB V12.068.