Genome-Wide Association Analysis for Phosphorus Use Efficiency Traits in Mungbean (Vigna radiata L. Wilczek) Using Genotyping by Sequencing Approach

Mungbean (Vigna radiata L. Wilczek) is an annual grain legume crop affected by low availability of phosphorus. Phosphorus deficiency mainly affects the growth and development of plants along with changes in root morphology and increase in root-to-shoot ratio. Deciphering the genetic basis of phosphorus use efficiency (PUE) traits can benefit our understanding of mungbean tolerance to low-phosphorus condition. To address this issue, 144 diverse mungbean genotypes were evaluated for 12 PUE traits under hydroponics with optimum- and low-phosphorus levels. The broad sense heritability of traits ranged from 0.63 to 0.92 and 0.58 to 0.92 under optimum- and low-phosphorus conditions, respectively. This study, reports for the first time such a large number of genome wide Single nucleotide polymorphisms (SNPs) (76,160) in mungbean. Further, genome wide association study was conducted using 55,634 SNPs obtained by genotyping-by-sequencing method. The results indicated that total 136 SNPs shared by both GLM and MLM models were associated with tested PUE traits under different phosphorus regimes. We have identified SNPs with highest p value (–log10(p)) for some traits like, TLA and RDW with p value (–log10(p)) of more than 6.0 at LP/OP and OP condition. We have identified nine SNPs (three for TLA and six for RDW trait) which was found to be present in chromosomes 8, 4, and 7. One SNP present in Vradi07g06230 gene contains zinc finger CCCH domain. In total, 71 protein coding genes were identified, of which 13 genes were found to be putative candidate genes controlling PUE by regulating nutrient uptake and root architectural development pathways in mungbean. Moreover, we identified three potential candidate genes VRADI11G08340, VRADI01G05520, and VRADI04G10750 with missense SNPs in coding sequence region, which results in significant variation in protein structure at tertiary level. The identified SNPs and candidate genes provide the essential information for genetic studies and marker-assisted breeding program for improving low-phosphorus tolerance in mungbean.

Mungbean (Vigna radiata L. Wilczek) is an annual grain legume crop affected by low availability of phosphorus. Phosphorus deficiency mainly affects the growth and development of plants along with changes in root morphology and increase in rootto-shoot ratio. Deciphering the genetic basis of phosphorus use efficiency (PUE) traits can benefit our understanding of mungbean tolerance to low-phosphorus condition. To address this issue, 144 diverse mungbean genotypes were evaluated for 12 PUE traits under hydroponics with optimum-and low-phosphorus levels. The broad sense heritability of traits ranged from 0.63 to 0.92 and 0.58 to 0.92 under optimum-and low-phosphorus conditions, respectively. This study, reports for the first time such a large number of genome wide Single nucleotide polymorphisms (SNPs) (76,160) in mungbean. Further, genome wide association study was conducted using 55,634 SNPs obtained by genotyping-by-sequencing method. The results indicated that total 136 SNPs shared by both GLM and MLM models were associated with tested PUE traits under different phosphorus regimes. We have identified SNPs with highest p value (-log 10 (p)) for some traits like, TLA and RDW with p value (-log 10 (p)) of more than 6.0 at LP/OP and OP condition. We have identified nine SNPs (three for TLA and six for RDW trait) which was found to be present in chromosomes 8, 4, and 7. One SNP present in Vradi07g06230 gene contains zinc finger CCCH domain. In total, 71 protein coding genes were identified, of which 13 genes were found to be putative candidate genes controlling PUE by regulating nutrient uptake and root architectural development pathways in mungbean. Moreover, we identified three

INTRODUCTION
The genus Vigna includes about 100 species dispersed in subtropical and tropical region globally (Tomooka et al., 2002;Maxted et al., 2004). Among these green gram, black gram (Vigna mungo), cowpea (Vigna unguiculata), adzuki bean (Vigna angularis), moth bean (Vigna aconitifolia), creole bean (Vigna glabrescens), and tuber cowpea (Vigna vexillata) have been domesticated (Takahashi et al., 2016). The wild species are adapted to acid soils, deserts, limestone rocks, sandy beaches, and wetlands (Tomooka et al., 2011). Mungbean is native of India and Central Asia. The closest relative of mungbean is Vigna radiata var. sublobata (Roxb.) Verdc., occurring in wild state in Western Ghats and Sub-Himalayan tract. From India, mungbean spread to Southeast Asia and through silk road (Tomooka et al., 1992). Mungbean is cultivated in over 7 mha with production of about 5 mt in the world (Nair et al., 2019). The crop is generally grown in East, South, and Southeast Asia. Mungbean is a rich source of protein, minerals, and vitamins (Mubarak, 2005). Mungbean is consumed as dhal (whole grain or split) in South Asia and as sprouts, noodles, and paste in East Asia.
Phosphorus (P) is required for photosynthesis, redox reactions, energy production, and phosphorylation/ dephosphorylation. The uptake of inorganic P by plants ranges from 10 to 25% (Johnston et al., 2014). Despite being abundant in agricultural soils, the majority of P occurs in an insoluble form (Menzies and Lucia, 2009). P is fixed in acidic soils by metal oxides and by carbonate compounds in alkaline soil or it exists as organic P (Heuer et al., 2017). Low bioavailability of P in soil is due to slow diffusion and high fixation. It is reported that 50% of the world's arable land is P deficient (Lynch, 2011). Sedimentary deposits are main source of P reserves and account for 80-90% P production. The phosphate rocks are generally located in China, Western Sahara, Morocco, and Russia (Kauwenbergh et al., 2013). These nonrenewable, irreplaceable, and finite global phosphate reserves are depleting and are expected to exhaust by the end of this century (Manschadi et al., 2014). The demand for P is increasing to meet the food security of growing population. Environmental problems like water eutrophication and soil degradation result from overuse of phosphorus.
The natural P fertilizer is a limited resource. It is imperative to enhance P fertilizer use efficiency to develop plants with higher uptake and utilization of phosphorous. Phosphorus use efficiency (PUE) measures ability of plant for phosphorus uptake and its use for production of biomass or yield (Wiel et al., 2016). Under P starvation condition, plants activate a range of mechanisms which results in higher uptake of P from soil (Ramaekers et al., 2010). Plant traits increasing P uptake ability includes change in root architecture such as increased lateral root number, root hair length and density, and high root-toshoot ratio (Lin et al., 2009;Peret et al., 2011). The modification of root architecture along with the exudation of P mobilizing organic acids which results in enhancement of P uptake by plants (Lambers et al., 2006). It is established that vigorous root system with increased nutrient acquisition ability can lead to yield increase under optimized nutrient management system (Vinod and Heuer, 2012). The Pup1-specific protein kinase gene and phosphorus-starvation tolerance 1 is reported to be enhancer of early root growth which leads to more acquisition of P in plants (Gamuyao et al., 2012).
The quantitative trait loci (QTLs) for PUE traits in response to P starvation have been reported in crop plants such as wheat (Su et al., 2009;Yuan et al., 2017), rice (Luo et al., 2017;Mahender et al., 2018), corn (Zhu et al., 2005;Chen et al., 2008), common beans (Yan et al., 2004), and soybean (Liang et al., 2010). Conventional linkage mapping utilizes bi-parental mapping populations developed through controlled crosses. Identification of QTLs with small effect requires large population size in biparental mapping (Collard et al., 2005). QTL mapping covers small portion of the genetic architecture for a trait as only the allele differing between the parents segregate. Furthermore, biparental mapping has less generation and offspring to extensive shuffling of genome and QTLs are localized to a large part of chromosomal regions (10 to 20 cM) (Flint-Garcia et al., 2003).
The association mapping (AM) utilizes panel of diverse genotypes for genome-wide association studies for identification of trait-linked molecular tags. QTLs are identified based on association of marker with trait elucidated by LD between makers and polymorphisms across the AM panel. AM has advantage of the populations undergoing many generations of recombination. Thus, AM provides much finer resolution than genetic mapping. Genome-wide association study (GWAS) have evolved over the last few years into a powerful tool which exploits hundreds of markers across all chromosomes to find genetic variations associated with particular trait (Hirschhorn and Daly, 2005). GWAS and high-throughput next-generation sequencing (NGS) platform make AM as best approach for QTL mapping in plants (Atwell et al., 2010). The genotyping by sequencing (GBS) is a simple and attractive high-throughput multiplexed sequencing method that involves discovery and genotyping of marker (Poland and Rife, 2012). The robustness and low cost of GBS makes this as an outstanding tool for many applications in plants. Since its inception, GBS assay has been successfully applied to develop genetic linkage maps, GWAS, and fine mapping/map-based cloning of QTLs in diverse crops, including soybean (Bastien et al., 2014;Jarquin et al., 2014) and chickpea Upadhyaya et al., 2016).
Single nucleotide polymorphism (SNP) markers are highly desirable for marker-assisted selection (MAS) due to their abundance, locus specificity, codominant inheritance, and highthroughput genotyping (Terakami et al., 2011). Significant SNPs associated with PUE-related traits have been reported in rice (Wissuwa et al., 2015), maize (Wang et al., 2019), cowpea (Ravelombola et al., 2017), oilseed rape (Wang et al., 2017), and soybean (Ning et al., 2016). To our knowledge, to date, a SNP marker associated with PUE traits in mungbean has not been reported. Keeping the importance of P uptake and utilization efficiency traits in mind, the present study aimed to decipher the complex genetic architecture of PUE in mungbean. The study will also lead to advancement of mungbean genomics by enriching the mungbean genome with the development of SNP markers.

Constitution of Association Mapping Panel and Growing Condition in Hydroponics System
The AM panel for GWAS comprised 144 mungbean genotypes including 40 released varieties, 42 advanced breeding lines (ABL), and 62 exotic germplasm lines (Supplementary Table S1). The experiment was conducted under controlled environment facility at the Indian Agricultural Research Institute, New Delhi, India using hydroponics approach. The seedlings were evaluated in completely randomized design with three replications. In each replication, eight plants per genotype were evaluated. Air temperature in the greenhouse was maintained at 30/18 • C (day/night), photoperiod for 12 h, and the relative humidity at 90%. In order to screen genotypes, surfaces of seeds were sterilized with 0.1% (w/v) HgCl 2 for 3 min and washed thrice with distilled water to remove detergent. The seeds were kept for germination in germination paper. After emergence of cotyledonary leaves, uniformly germinated seedlings were transferred to modified Hoagland solution with two levels of P. Two P levels were balanced using KH 2 PO 4 as optimum P (250 µM) and low P (3 µM). The composition comprised CaCl 2 ·2H 2 O (0.75 mM), K 2 SO 4 (0.92 mM), Fe-EDTA (0.04 mM), urea (5 mM), MgSO 4 (1 mM), and micronutrients [H 3 BO 3 (2.4 µM), ZnSO 4 (0.6 µM), MnSO 4 (0.9 µM), CuSO 4 (0.62 µM), and Na 2 MoO 4 (0.6 µM)] (Sivasakthi et al., 2017). Nutrient solution was renewed at 2 days interval, and pH was maintained around 6.0 using 1 M HCL or 1 M KOH.

Phenotypic Evaluation for Phosphorous Use Efficiency Trait in Mungbean
Twenty-one-day-old seedlings grown under optimum-and low-P conditions were carefully separated into leaf, stem, and root to record observations on PUE traits. Isolated complete root system of each plant was placed on a tray with no overlapping of roots for recording observations on root traits. Root traits recorded include total root volume (TRV), total root length (TRL), root average diameter (RAD), total root surface area (TSA), total root tips (TRT), and root forks (RF). These were measured using a root scanner (EPSON XL 1000) with WinRhizo software (Pro version 2016a; Regent Instrument Inc., Quebec, QC, Canada). Primary root length (PRL) was determined manually using meter scale. Chlorophyll concentration (CHL) was measured by using MC-100 chlorophyll concentration meter (Apogee Instruments, Inc., United States). Total leaf area (TLA) was measured using leaf area meter (LI-COR 3000, Lincoln, NE, United States). Separated shoots and roots were dried at 60 • C until constant mass of root dry weight (RDW) and shoot dry weight (SDW) was obtained and then root-to-shoot ration (RSR) was calculated. All the traits have been measured in three replications. The PRL and TRL have been measured in centimeters per plant. TSA, TRV, and RAD have been measured in square centimeters per plant, cubic centimeters per plant, and centimeters, respectively. CHL and LA have been measured in micromoles per square meter and square meters per plant, respectively. RDW and SDW have been measured in grams per plant ( Table 1). Descriptive statistics, Pearson correlation coefficients, and analysis of variance were performed using the Statistical Tool for Agricultural Research (STAR) 2.1.0 software (Gulles et al., 2014). The broad sense heritability is defined by the formula Vg/Vg + Ve, where Vg and Ve are the genotype and environmental variances, respectively.

Discovery, Genotyping, Structural, and Functional Annotation of GBS-Based SNPs
The genomic DNA was isolated from the AM panel genotypes. The isolated genotypes were used to constitute two 144-plex GBS library and sequenced (2 × 150-bp single end) using Illumina HiSeq 4000 NGS platform according to Liu et al. (2012) and Bastien et al. (2014). Three accessions were taken as biological replicates. The demultiplexing and mapping of FASTQ reads onto the reference mungbean genome (Kang et al., 2014) had been carried out. The high-quality SNPs were detected using Bowtie v2.1.0 (Langmead and Salzberg, 2012) and reference-based GBS pipeline/genotyping approach of STACKS v1.0 1 according to Kujur et al. (2015). The unmapped sequence reads were analyzed through de novo-based GBS approach of STACKS. Subsequently, the putative SNPs were discovered among genotypes following Kujur et al. (2015). In order to validate SNPs identified through GBS assay, the SNPs present in three candidate genes have been selected. The candidate genebased SNPs have been validated through Sanger sequencing. Therefore, the primers were designed from 300-bp sequences which were flanking either side of selected SNPs. Next, the PCR amplification carried out using the genomic DNA of 25 mungbean accessions. The 25 mungbean accessions have been selected on the basis of response of root to available phosphorous in soil (Reddy et al., 2020). The purification and sequencing of PCR products were carried out. The alignment of sequences is performed to detect SNPs among accessions according to Kujur et al. (2013). The reference genome GBS-based SNPs identified in genic and intergenic component of genomes were annotated according to mungbean genome annotation (Kang et al., 2014).

Construction of Phylogenetic Tree, Determination of Population Structure, and LD Patterns
The genotyping data among mungbean genotypes were analyzed with PowerMarker v3.51 (Liu and Muse, 2005) and MEGA v6.0 (Tamura et al., 2013). An unrooted neighbor-joining (NJ)based phylogenetic tree (with 1000 bootstrap replicates) had been constructed. The population genetic structure among genotypes was determined by a model-based program STRUCTURE v2.3.4 according to Kujur et al. (2015). The genotyping data of chromosome-based SNPs was analyzed by PLINK and the fullmatrix approach of TASSELv5.0 (Saxena et al., 2014). The LD patterns (r 2 ) and LD decay (by plotting average r 2 against 50 kb uniform physical intervals across chromosomes) had been determined. The LD had been measured as frequency correlation among pair of alleles across a pair of SNP loci.

Marker Trait Association Analysis
The genome-wide SNP genotyping information of mungbean AM panel was integrated with PUE trait phenotyping data, ancestry coefficient (Q matrix), and relative kinship matrix (K) data using general linear model (GLM) (Q model) and mixed linear model (MLM) approach. The quantile-quantile plot had been generated to determine relative distribution of observed and expected log 10 (p) value for each associated SNP. The adjusted p value threshold of significance was corrected for multiple comparisons according to false-discovery rate (FDR cutoff ≤ 0.05) (Benjamini and Hochberg, 1995). The SNP loci having significant association with PUE trait at highest R 2 (degree of associated SNP) and lowest FDR adjusted p values (threshold p < 1 × 10 −5 ) were identified by integrating all the model-based outputs of TASSELv5.0.

Digital Gene Expression Analysis for the Identified Candidate Gene
The digital gene expression analysis has been performed for the candidate genes which have been identified in our study. We have taken the Arabidopsis ortholog of candidate gene from mungbean genome for analysis. The gene expression pattern has been searched using Expression Angler, an online search tools (Austin et al., 2016). The experiment had been carried out using different developmental stages like root, shoot, xylem etc.

Phenotypic Variation for 12 PUE Traits
Twelve PUE traits were investigated in 144 genotypes of the AM panel under two levels of P conditions. Analysis of variance revealed highly significant variation among 144 genotypes for all tested traits under two P regimes (Supplementary Table S2). We recorded highly significant phenotypic difference between two P regimes (p < 0.01) for all traits except PRL, TRV, and TSA, while RDW showed general significance (p < 0.05) ( Table 2). The mean values of TRL, TSA, TRT, RF, TLA, and SDW were significantly higher under optimum-P condition compared with low-P condition. Whereas, mean values of RAD, CHL, RDW, and RSR were significantly higher in low-P compared with optimum-P condition. The coefficient of variation varied from 3.59 to 15.47% and from 3.32 to 24.31% under optimum-and low-P conditions, respectively. This study showed a maximum heritability (0.92) for TLA under optimum-P condition. At low-P condition, the heritability of CHL in leaf was found to be highest (0.92). The variation in all tested traits under optimum-P condition was significantly correlated with the corresponding traits in low-P condition with Pearson correlation coefficients ranging from 0.39 to 0.76. The trait ratios between low-P and optimum-P conditions ranged from 0.58 (TLA) to 1.78 (RSR) with six values being lower than 1.00 ( Table 3). The coefficient of variation for trait ratios ranged from 4.72 to 28.22%.

Large-Scale Discovery, High-Throughput Genotyping and Annotation of Genome-Wide GBS-and Candidate Gene-Derived SNPs
A total of 264.40 million high-quality sequence reads had been generated. The reads are evenly distributed (mean, 1.83 million reads) across 144 genotypes of AM panel. Notably, on an average, 75% of these reads were mapped StDev, standard deviation; CV, coefficient of variation; h 2 , broad sense heritability; PRL, primary root length; TRL, total root length; TSA, total root surface area; TRV, total root volume; RAD, root average diameter; TRT, total root tips; RF, root forks; CHL, chlorophyll concentration; TLA, total leaf area; RDW, root dry weight; SDW, shoot dry weight; RSR, root-to-shoot ratio; OP, traits tested under optimum-phosphorus conditions; LP, traits tested under low-phosphorus conditions. *p < 0.05; **p < 0.01; ***p < 0.001, levels of significance. a The traits tested in the study. b Independent t test was performed at 1% level of significance to test the significant difference between optimum-and low-P conditions for all tested traits.
on V. radiata reference genome. In total, 76,160 high-quality SNPs (with read-depth 10, <5% missing data, and 8% minor allelic frequency) differentiating 144 samples were discovered using mungbean reference genome-based GBS assay. The genome-wide GBS-based SNP genotyping overall identified 55,634 chromosome-based SNPs showing polymorphism among mungbean genotypes of AM panel. The distribution of all SNPs was found to be in different components of genome including intergenic (27%), intragenic (23%), regulatory (31%), and scaffold (19%) variants ( Figure 1A). A total of 20,526 SNPs identified in unanchored scaffold region ( Table 4). A maximum of 8,805 SNPs (15.8%) were mapped on chromosome 1, whereas, a minimum of 2,264 SNPs (4.01%) were mapped on chromosome 3. The SNP density was found to be low on chromosome 10 (1.27 SNPs per 100 bp). The structural annotation of 55,634 SNPs identified a number of 25,663 (46.12%) and 29,923 (53.78%) SNPs in the 11,068 protein coding genes (intragenic region) and intergenic regions of mungbean genome, respectively. The abundance of gene-based SNPs was found in the regulatory region (33,724 SNPs, 60.1%), followed by CDS (13,030 SNPs, 23.4%) and least (10,423 SNPs, 18.73%) in the intron region ( Figure 1B). In total, 5,643 (43.4% of coding SNP) and 7,387 (56.6% of coding SNP) missense and synonymous SNPs have been found. The relative distribution of GBSbased SNPs mapped on 11 mungbean chromosomes showing polymorphism among AM panel genotypes are presented in a Circos circular ideogram ( Figure 1C). The functional annotation of 55,634 SNP-carrying genes revealed that their correspondence lies to 2,481 growth, 7,741 development, 4,936 metabolismrelated proteins followed by 764 signal transduction protein.
We performed Sanger sequencing-based validation of SNPs present in three candidate genes. Two SNPs present in two candidate genes have been mined through Sanger sequencing (Supplementary Table S4). The structurally and functionally annotated genome-wide SNPs physically mapped on mungbean genome can be used for genotyping applications for quantitative trait dissection of mungbean.

Molecular Diversity, Linkage Disequilibrium, and Population Structure
The identified 55,634 genome-wide GBS-based SNPs revealed polymorphism among AM panel mungbean genotypes with higher nucleotide diversity (θπ, 0.35547; θω, 0.18039; and Tajima's D, 3.21848) potential. Any significant difference for nucleotide diversity value have not been found within the different set of genotype group such as, advanced varieties, landrace, etc. However, we have compared the diversity (θπ) value among the related species, like soybean and chickpea. It has been found that the diversity in soybean was θπ: 0.33 which was comparable to mungbean (θπ, 0.35). In chickpea, θπ was found to be 0.99, which is higher as compared with mungbean. The unrooted neighbor-joining phylogenetic tree construction using 55,634 genome-wide SNPs depicted a distinct differentiation among the AM panel of mungbean (Figure 2). These genotypes were clustered into three major groups. The LD estimates (r 2 ) and extent of LD decay were determined using 55,634 SNPs mapped on 11 mungbean chromosomes. We determined the LD decay of 55,634 SNP pairs by pooling the r 2  estimates across 11 chromosomes of mungbean and plotting their average r 2 against the 50-kb equal intervals of physical distance (maximum up to 1,000 kb). The determination of LD patterns in a population of AM panel genotypes using 55,634 SNPs (physically mapped across eleven mungbean chromosomes) showed a higher LD estimate (r 2 , 0.62) and less-extensive LD decay (r 2 decreased half of its maximum value) at 50-100 kb physical distance in mungbean chromosomes (Supplementary Figure S3A). A slight increasing and then consistent pattern of LD decay (r 2 ≥ 0.3) was observed with increase in the physical distance (cM) of SNP markers mapped on the 11 chromosomes. The population genetic structure among AM panel genotypes was determined using SNPs (mapped on 11 mungbean chromosomes) employing STRUCTUREv2.3.4 software. The admixture model-based STRUCTURE simulations were performed at different levels of possible population numbers (K = 1 to 10) with 20 replications. The population numbers (K) were determined using the second-order statistics of STRUCTUREv2.3.4 by delta K estimation. The results showed a maximum value of delta K at 3 (Supplementary Figure S3B), which confirms the classification of 144 mungbean genotypes into three genetically distinct population groups (POP I-III) with high-resolution population structure (Supplementary Figure S3C). The population I predominantly consists of ABL (23) followed by exotic (17) and released cultivars (9). In population II, exotic (25) genotypes were found to be abundant as compared with ABL (4) and released cultivars (7). In population III, the numbers of released cultivars (24) were found to be higher as compared with exotic (20) and ABL (15). Populations I, II, and III consist of 34, 25, and 41% of total accessions, respectively. In population I, ABL is dominant (47%), population II is dominated by exotic (69%), and population III is found to be dominated by released cultivars (41%). The grouping of all AM panel genotypes into three groups was further confirmed by principal component analysis (PCA) along with phylogenetic tree and population structure analysis (Supplementary Figure S4).

GWAS for Phosphorous Use Efficiency Traits in Mungbean
For GWAS, the genotyping data of 55,634 GBS-SNPs developed from genotypes of AM panel were integrated with phenotyping data, high-resolution population genetic structure, and PCA data. The Bonferroni-corrected threshold (-log(p) = 4.0) was used as a cutoff to identify the significant marker trait associations. A total of 1,171 and 139 significant SNPs were identified by GLM and MLM, respectively. Of these, 136 SNPs were identified by both models (Table 5 and Supplementary Table S3). PRL, primary root length; TRL, total root length; TSA, total root surface area; TRV, total root volume; RAD, root average diameter; TRT, total root tips; RF, root forks; CHL, chlorophyll concentration; TLA, total leaf area; RDW, root dry weight; SDW, shoot dry weight; RSR, root-to-shoot ratio. a The traits tested in the study. b The total number of significant SNPs detected at the threshold of −log(p) = 4.0. c The phenotypic variation revealed by ANOVA with total significant SNPs detected using the two models. d The number of significant SNPs detected by both models.
Under optimum-P condition (OP), a total of 597 SNPs identified by GLM were found to be associated with all tested traits in the study. The GLM detected 39 and 18 significant SNPs for RAD and CHL, respectively; for which no significant SNP was detected by MLM. The MLM-detected maximum number of SNP loci (13) was found to be associated with TRT with 13.93% contribution of phenotypic variation. The MLM detected 3, 4, 7, 9, 6, 10, 10, 6, and 1 significant SNPs for PRL, TRL, TSA, TRV, RF, TLA, RDW, SDW, and RSR, respectively (Figure 3 and Supplementary Figure S1). The phenotypic variation (R 2 ) explained by the SNPs ranged from 13.28 to 17.25%.
At low-P and optimum-P conditions, MLM detected a maximum number of SNP loci (10) associated with RDW trait with around 14.59 and 17.25% of average phenotypic variation, respectively. Whereas, at LP/OP condition, a number of five SNP loci were found to be associated with RDW trait. At LP/OP conditions, MLM detected a maximum number of SNP loci (18) found to be associated with TLA trait and contributed around 15% of phenotypic variation. Whereas, at OP and LP conditions, a number of 10 and 5 SNP loci were found to be associated with TLA traits, respectively. The association results were further confirmed by analyzing quantile-quantile plot (Supplementary Figure S2). In total, 136 SNP loci were found to be present in a number of 71 protein coding genes. These genes were found to have diverse functions, including in stress, metal ion transport, metabolism, and development. Interestingly, most of the genes identified were found to be involved in nutrient transporting activity. These results suggest the utility of implicating GWAS in mungbean to decipher potential candidate genes for PUE traits in mungbean.

Delineation of Putative Candidate Genes for PUE Traits in Mungbean
In order to delineate the putative candidate genes regulating PUE trait in mungbean, we retrieved the CDS sequences of 71 associated protein coding genes and performed BLAST 2 search against Arabidopsis genome database (with default parameter p value ≤ 1.0). The study identified a total of 13 genes with highest similarity and identity (≥80%) to the Arabidopsis genes. These genes were well characterized and found to be involved in different nutrient uptake and root architecture development pathway (Table 6). Therefore, we concluded that these genes may be putative candidate genes regulating PUE traits in mungbean. We reported a total of 1,171 and 139 significant SNPs associated with PUE traits by GLM and MLM models, respectively. Of these, a total of 136 SNP loci shared by both models were found to be associated with 71 protein-coding genes (Figure 5). Most of these genes are found to be involved in diverse role, including nutrient uptake, root architecture development, and abiotic stress in mungbean. The Arabidopsis orthologous genes are found to be involved mostly in root development and nutrient uptake. Among these 71 protein coding genes, a total of 13 Chr, chromosome; RF, root forks; TRT, total root tips; TLA, total leaf area; RAD, root average diameter; RDW, root dry weight; SDW, shoot dry weight; RSR, root-toshoot ratio; TSA, total root surface area; TRV, total root volume; LP, low-phosphorus condition; OP, optimum-phosphorus condition; LP/OP, ratio of low phosphorus to optimum-phosphorus condition. a The key candidate genes associated with phosphorus use efficient traits under different phosphorus conditions.
genes are found to be strongly involved in root development and nutrient uptake (Figure 5). The Arabidopsis orthologs of these genes are plastidal glycolate/glycerate translocator 1 (PLGG1) (At1G32080), At5PTASE9 (At2G01900), Ca 2+ -sensor protein kinases-NIN-like protein 7 (CPK-NLP7) (At2G37590), cellulose synthase-like D gene (CSLD3) (At5G16910), K + channel Arabidopsis thaliana (KAT1) (At4G11480), ubiquitin-conjugating enzyme E2 13 (UBC13) (At3G46460), fibroblast growth factor receptor substrate 3 (FRS3) (At2G27110), GTPase root hair defective 3 (RHD3) (At4G39250), iron-regulated transport er1 (IRT1) (At5G58440), PHO80 (At5G64260), K + uptake 7 (KUP7) (At5G09400), speechless (SPCH) (At5G53210), and KAN protein (KANADI) (At4G25440). We found that three genes VRADI11G08340, VRADI01G05520, and VRADI04G10750 containing missense SNP at CDS region. These three genes may regulate PUE traits by changing at protein structure level. Therefore, we determined the protein structure by I-TASSER 3 for these three genes in order to decipher the functional complexity at protein structure level. However, we were able to predict the protein structure for two genes (VRADI11G08340; VRADI01G05520) (Figure 4). The structural analysis of these genes at protein level revealed a significant variation for native and mutated version of the genes. The missense SNPs present in the CDS of VRADI11G08340, VRADI01G05520, and VRADI04G10750 genes changes the amino acids from A/G; methionine to valine, T/C; tyrosine to histidine and T/C; and leucine to serine, respectively. Therefore, it may be concluded that these genes are potential candidate genes regulating PUE in mungbean by changing protein structure at tertiary level. Furthermore, we found that 20 potential candidate genes tended to have pleiotropic effects (Supplementary Table S3).
For example, candidate gene VRADI06G09990 was significantly associated with RDW and RSR under low-P condition, with RDW under optimum-P condition. Similarly, candidate gene VRADI07G06230 was significantly associated with TLA, RDW, SDW, TRV, and TSA under optimum-P condition. These genes may be involved in many metabolic pathways and affect the PUE traits and therefore worthy of future attention.

Phenotypic P Response of Mungbean Genotypes
The present study showed significant greater root biomass accumulation resulting in higher root-to-shoot ratio under low-P condition. Under P-deficient condition, P-efficient genotype maintained higher root-to-shoot ratio than P-inefficient genotype in maize (Mollier and Pellerin, 1999) and common bean (Nielsen et al., 2001). RDW has been indicated as key trait for determining P starvation response in soybean (Li et al., 2005). Higher chlorophyll concentration in leaves causes dark green foliage under low-P condition. P deficiency causes stunted growth and foliage turns dark green color because of accumulation of starch and sugars in the leaves (Pandey, 2015).

GBS Assay Accelerate Genomic Resources Development in Mungbean
In the present study, we employed GBS assay for discovery and genotyping of SNPs in mungbean. The length of pseudomolecule assemblies of mungbean genomes was 333.4 Mb representing 57% of total genome. The SNPs mined with GBS assay were found to be dispersed in different portions of mungbean genomes, which may lead to differences in physical map density of SNPs within the genomes. This large number of SNPs (76,160) reported in our study can be utilized for genotyping and QTL mapping for several traits as well as comparative genome mapping involving mungbean and legumes. Additionally, 55,634 chromosomebased SNPs can be used for marker-based applications of mungbean breeding.

Functional Utility of GBS-Based SNPs in Mungbean
The discovered SNPs were annotated in different components of protein-coding genes of mungbean genomes. The structural and functional annotation revealed that a total of 29,923 SNPs present in intergenic region and 25,663 SNPs in intragenic region. Among the intragenic variants, 10,423 and 13,030 SNPs were found to be present in intron and CDS regions, respectively. Further, the coding variants represented 5,643 and 7,387 missense and synonymous variants. A large number of SNPs were found to be present in regulatory region. The upstream and downstream regulatory region consisted of likely equal number of SNPs (16,260 upstream SNPs,17,464 downstream SNPs). The structural and functional annotation of genes will accelerate the pace of AM to identify a number of genes/QTLs governing qualitative and quantitative agronomic traits in mungbean. Tajima's D is a statistical test which is performed to test the null hypothesis of mutation-drift equilibrium and constant population size. A significantly higher Tajima's D (3.21848) value were observed in our used AM panel. The results indicate that mungbean may undergo for a population reduction (Maruyama and Fuerst, 1985), population subdivision, or a recent bottleneck.
In our study, high-resolution LD pattern were determined in a structured population of 144 mungbean accessions. In soybean and chickpea, several studies have been carried out using GBS and WGRS (Whole genome re-sequence) approach to find the LD pattern and dissect complex trait Das et al., 2015;Sonah et al., 2015;Zhou et al., 2015;Varshney et al., 2019). In a study, the LD extent was found to be ∼27 kb in wild soybean. In the soybean landraces and improved cultivars, LD was found to be 83 and 133 kb, respectively, which are comparable with mungbean (Zhou et al., 2015). The LD pattern in mungbean is found to be different from chickpea with larger decay at physical position (150-200 kb) (Saxena et al., 2014;Kujur et al., 2015;Upadhyaya et al., 2016;Varshney et al., 2019).
The study revealed a higher LD estimate (r 2 , 0.62) and less extensive LD decay with approximately 50-100 kb physical distance in mungbean chromosomes. In a previous study, the LD decay was found to be extensive with approximately 60-100 kb physical distance in a panel of cultivated and wild mungbean accessions (Noble et al., 2018). The LD pattern in mungbean was found to be likely the same to other selfpollinated crop species, like soybean (Patil et al., 2016) but different from chickpea (Saxena et al., 2014;Kujur et al., 2015) with larger decay at physical position. The population structure analysis of the composed panel revealed three population groups which was further confirmed by phylogenetic and principal component analysis. Additionally, this is the first high-resolution determination of LD decay in mungbean with such a large number of genome-wide markers (55,634 SNP markers).

GWAS to Dissect Complex PUE Traits in Mungbean
The PUE is an important trait which determines yield potential in several crops, including mungbean. Crops with enhanced PUE lead to higher growth for the same amount of P taken up at a time (Rengel and Marschner, 2005). However, the trait is complex and quantitative in nature (Heuer et al., 2017). Therefore, it is desirable to dissect complex genetic inheritance patterns of PUE traits in mungbean. Efforts have been made toward molecular dissection of this complex trait. Several QTLs/alleles for root, shoot, and grain yield-related traits underlying low-phosphorus tolerance have been reported in different legume crop plants including common bean (Beebe et al., 2006;Ochoa et al., 2006;Diaz et al., 2017), Soybean (Song et al., 2014;Zhang et al., 2017) and cowpea (Ravelombola et al., 2017). However, AM has not been explored in mungbean genome to identify functionally relevant molecular tags regulating important traits. In a previous study, a total of 22,000 polymorphic genome-wide SNPs were developed and utilized for AM for seed coat color in mungbean (Noble et al., 2018). In this study, for the first time, such a large number of high-quality SNPs has been developed (76,160 SNPs) in mungbean. Furthermore, these genome-wide welldistributed SNPs have been utilized to perform AM for PUE traits in mungbean. We reported a total of 1,171 and 139 significant SNPs associated with PUE traits by GLM and MLM models, respectively. Of these, a total of 136 SNP loci shared by both models were found to be associated with 71 protein coding genes ( Figure 5). Most of these genes are found to be involved in diverse roles, including nutrient uptake, root architecture development, and abiotic stress in mungbean. The Arabidopsis orthologous genes are found to be involved mostly in root development and nutrient uptake.
In our study, we identified bHLH, MYB, and F-box domaincontaining and cellulose synthase-like genes. The genes having basic helix-loop-helix DNA binding domains are reported to be involved in root development and biomass accumulation. These genes may help to increase P content and root hair formation (Chen et al., 2007;Li et al., 2011;Giehl and von Wiren, 2014). The genes with MYB DNA binding domain is reported to act in controlling architecture of root as well as adaptation to P starvation (Devaiah et al., 2009;Dai et al., 2012;Baek et al., 2013). The F-box domain-containing genes are reported to influence response to low-P stress in crop plants (Baldwin et al., 2008;Perez-Torres et al., 2008;Nakamura et al., 2009). The cellulose synthase-like genes are reported to be involved in cell wall remodeling under low-P condition (Xu et al., 2018). Interestingly, most of the genes exhibited root-specific expression which clearly indicates the possible role of these genes in phosphorous nutrient uptake.
Among these 13 putative candidate genes, we identified missense SNPs in CDS region of three genes, VRADI11G08340, VRADI01G05520, and VRADI04G10750. The Arabidopsis ortholog genes of these are plastidal glycolate/glycerate translocator 1 (PLGG1), inositol polyphosphate 5-phosphatase (At5PTASE9), and CPK-NLP7. The PLGG1 encodes a chloroplast protein which regulates lateral root development in response to ABA (Dong et al., 2018). Trull et al. (1997) described the role of ABA in root development under P deficiency with ABAdeficient mutants in Arabidopsis. Interestingly, the AT5PTASE9, inositol polyphosphate 5-phosphatase is previously reported to be involved in phosphate homeostasis (Kuo et al., 2018;Zhu et al., 2019). The AT5PTASE9 is a member of the AT5PTASE family, which regulates several responses, such as the calcium influx, ROS production, bulk endocytosis, and induction of stress-responsive genes (Golani et al., 2013). The CPK-NLP7 is found to be involved in a transcriptional regulatory network that regulates the root and shoot architecture in response to changes in nitrogen availability (Gaudinier et al., 2018). The missense SNPs present in the CDS region resulted in structural changes at protein level. Allelic variation of VRADI11G08340, VRADI01G05520, and VRADI04G10750 between native and mutated sequences suggests that the respective protein differ in domains and structure for post-transcriptional modification (Figure 6). The differences in protein-protein interactions and signal integration results in differential transcriptional modulations that are causal for observed difference in PUE (Wissuwa et al., 2015). Therefore, these genes may be the potential candidate genes regulating PUE trait in mungbean by changing the protein structure at tertiary level.

CONCLUSION
We performed GWAS for PUE traits in a population containing 144 diverse mungbean genotypes with 55,634 SNP markers. The genetic diversity analysis, LD analysis, and population structure of this study will accelerate the pace of mungbean genomics research which can be directed toward genetic improvement of mungbean. The identified SNPs and candidate genes associated with PUE traits can be used in mungbean breeding program to enhance phosphorus deficiency tolerance. These results demonstrated the efficacy and reliability of implementing highresolution AM in mungbean which has a narrow genetic base and least availability of genomic resources. However, these genes need to be validated further to shed light on their functional role in regulation of the PUE traits in mungbean. Overall, the large number of genomic resources developed in our study will lead to enrich mungbean genomics. These resources can also be utilized to advance genomic-assisted breeding in mungbean.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are publicly available. This data can be found here: https://www.ncbi. nlm.nih.gov/sra/PRJNA609409.

AUTHOR CONTRIBUTIONS
HD, GM, RP, MS, and VR planned and designed the research. VR, MA, SM, and AS performed the hydroponics experiment. VR, P, TB, KT, and PG helped in data recording. VR, MA, and GM conducted genotyping of the association mapping panel. VR and SD conducted the STRUCTURE, TASSEL, and Bioinformatic analysis. HD, VR, and SD prepared the manuscript. RN, TS, and SK edited the manuscript for publication. All authors have read and approved the final manuscript.

FUNDING
This research work was carried out from the grants received from ICAR-IARI (research grant from the NAHEP-CAAST project) and DST (CRG/2018/002642). The funds from DST were utilized for phenotyping work and ICAR-IARI grants were utilized for GBS. The funding support for open access publication fees is being contributed by AVRDC.

ACKNOWLEDGMENTS
The authors are thankful to ICAR-IARI (research grant from the NAHEP-CAAST project) and DST for financial support to this research study. The authors acknowledge the facilities and guidance provided by the Head, Division of Genetics, Joint Director Research and Director, ICAR-IARI. The authors are grateful to the Head, Division of Plant Physiology and Professor, Division of Genetics for help and guidance during the investigations. VR is also thankful to Acharya N. G. Ranga Agricultural University, Andhra Pradesh for according deputation to pursue doctoral degree program at IARI, New Delhi.