Genetic Dissection of Phosphorous Uptake and Utilization Efﬁciency Traits Using GWAS in Mungbean

: Mungbean ( Vigna radiata L. Wilczek) is an early maturing legume grown predominantly in Asia for its protein-rich seeds. P deﬁciency can lead to several physiological disorders which ultimately result in a low grain yield in mungbean. The genetic dissection of PUpE (P Uptake efﬁciency) and PUtE (P utilization efﬁciency) traits are essential for breeding mungbean varieties with a high P uptake and utilization efﬁciency. The study involves an association mapping panel consisting of 120 mungbean genotypes which were phenotyped for total dry weight, P concentration, total P uptake, and P utilization efﬁciency under low P (LP) and normal P (NP) conditions in a hydroponic system. A genotyping-by-sequencing (GBS) based genome-wide association study (GWAS) approach was employed to dissect the complexity of PUpE and PUtE traits at the genetic level in mungbean. This has identiﬁed 116 SNPs in 61 protein-coding genes and of these, 16 have been found to enhance phosphorous uptake and utilization efﬁciency in mungbeans. We identiﬁed six genes with a high expression (VRADI01G04370, VRADI05G20860, VRADI06G12490, VRADI08G20910, VRADI08G00070 and VRADI09G09030) in root, shoot apical meristem and leaf, indicating their role in the regulation of P uptake and utilization efﬁciency in mungbean. The SNPs present in three genes have also been validated using a Sanger sequencing approach. the hotspot region several major QTLs/genes regulating P uptake and P utilization efﬁciency traits in the mungbean. Therefore, chromosome six can be enriched and targeted for genomics assisted breeding to improve the phosphorous deﬁciency response in mungbean. The LD analysis revealed a comparatively short LD decay in mungbean. This study leads towards deciphering the genetic complexity for PUpE and PUtE in mungbean. The genetic and genomic resources generated in this study can be exploited to develop phosphorous stress responsive cultivars which can overcome P deﬁciency in mungbean.


Introduction
Mungbean (Vigna radiata L. Wilczek), a nutrient-dense grain legume is predominantly cultivated in India, Myanmar, Pakistan, Bangladesh, Thailand, China, Indonesia, and Vietnam. Globally in 2020, this crop was grown in 7.3 million hectares with a production of 5.3 million tones [1]. Mungbean seeds are rich in protein, carbohydrates, minerals, and vitamins. Mungbean grains are easy to cook and possess digestible protein [2]. The mungbean grains are consumed as split grains/dhal (with or without seed coat), whole grains, and sprouts [3] and are an important constituent of several traditional preparations [4]. This crop is grown as an intercrop with corn or millets in East Asia and between wheat and rice in South and Southeast Asia [5]. Worldwide, mungbean is the preferred crop for intercropping due to its nitrogen-fixing ability and short maturity duration. The crop residue of mungbean improves soil fertility and can also be used as animal feed [6].
Phosphorus (P) is an important macronutrient essential for energy metabolism, nucleic acid synthesis, membrane stability, photosynthesis, and disintegration of carbohydrates and is involved in biological and key physiological processes of the plant [7]. P exhibits slow diffusion in soil and a high fixation by soil minerals [8]. In most soil conditions, P exhibits low mobility and poor availability to plants compared to other nutrients [9]. Nearly 70% of the global arable land is reported to be deficient in P [10,11]. The availability of P in deficient soils is around 1.0 µmol/L and the optimum requirement of P for plants is nearly 30 µmol/L [12]. The occurrence of P is in both organic and inorganic forms. In acidic soil, P reacts with iron, aluminum, and manganese oxides, whereas in alkaline soils it makes a complex with carbonate of calcium. Plants counter P scarcity by adaptive changes in root morphology, distribution, and topology. P deficiency reduces the growth and yield of plants. P-deficient plants exhibit increased root elongation, root hairs, root branching, and root/shoot ratios [13]. P is a major yield-limiting factor in subtropical and tropical environments [14].
To overcome P deficiency in the plant, phosphatic fertilizers are applied to the soil. The raw material for phosphatic fertilizer is derived from rock phosphate, which is only minable in a few countries of the world [15]. The single country, Morocco possesses nearly 85% of the remaining resources of rock phosphate [16]. Further, the demand for P is expected to increase annually by 2.2% during 2015-2020 [17]. In addition, the uptake of applied P is only about 15-30% [18]. However, the development of cultivars with improved phosphorous use efficiency (PUE) is considered as the sustainable approach to address this problem. PUE is defined as the total biomass production per unit of P uptake [19]. PUE is reported as a quantitative trait that is governed by root and shoot architectural traits [20]. PUE can be differentiated into PUpE and PUtE [21]. The phosphorous uptake efficiency (PUpE) refers to the amount of total P taken up by the plant and is expressed as mg P per plant. The phosphorous utilization efficiency (PUtE) refers to the mobilization of P within the plant for sustainable biomass production and is expressed as the ratio of plant biomass produced per unit of P taken up [21]. Development of mungbean cultivars with high phosphorus uptake and utilization efficiency is aimed to reduce the overall application of phosphatic fertilizer in the soil and also to improve the utilization of soil P. Further, various traits governing PUE are polygenic and are controlled by quantitative trait loci [22,23].
A combinatorial genomics-assisted breeding strategy is used in the application of findings to improve varieties in a targeted and efficient manner in crop plants. For the identification of QTLs, biparental mapping (BPM) and association mapping (AM) are the two most commonly used approaches that are being followed in several crop plants. The AM provides a higher mapping resolution over the BPM approach due to the occurrence of more recombination events over the generations in any given population [24]. In recent years, the evolution of next-generation sequencing (NGS) technology using a high throughput sequencing approach has drastically reduced the sequencing cost [25]. Single nucleotide polymorphisms (SNPs) are bountiful in the plant genome and known to influence the phenotype if located within the exon [26]. SNPs have the potential for exploitation in QTL mapping, increasing marker density, and high throughput marker-assisted selection [27].
In mungbean, the publication of a reference genome sequence of the cultivar VC1973A has provided a promising route for extensive genomic research in mungbean [28]. Also, genotyping-by-sequencing (GBS) is considered as a NGS-based genotyping methodology that has been used efficiently for discovering genome-wide SNPs [27,29]. The genome-wide association study (GWAS) is an established tool to scan markers across the genome and to identify the markers associated with the trait of interest [30]. GWAS combined with a high throughput genotyping platform enables association mapping as an excellent approach for detecting significant SNPs and candidate genes associated with a particular trait [31]. By employing GWAS, SNPs have been associated with PUE regulating traits in different grain legume crops like soybean [32], cowpea [33], and mungbean [20]. With this backdrop, this study aimed to (i) identify the genome-wide SNPs; (ii) study the genetic diversity and population structure; and (iii) GWAS for various PUpE and total PUtE-related traits in mungbean.

Plant Materials and Experimental Conditions
The investigation was conducted using 120 diverse genotypes (Supplementary Table S1) having various ABLs (advanced breeding lines), released varieties, and exotic germplasm lines as obtained from World Vegetable Centre (Taiwan). The experiment was performed under a hydroponic system to study the P uptake and P utilization efficiency in the selected mungbean genotypes. The plants were grown in the controlled greenhouse at the Indian Agricultural Research Institute, New Delhi having a day/night temperature of 30/18 • C, 90% RH (relative humidity), and a 12 h photoperiod. The seedlings were evaluated in a completely randomized design with three replications. In each replication, eight plants per genotype were evaluated. The seeds of all genotypes were first sterilized using HgCl 2 (0.1% w/v) and then kept for germination. On the fifth day, when cotyledonary leaves emerged, the seedlings were moved to hydroponic trays with a modified Hoagland solution. The nutrient solution was composed of K 2 SO 4 (0.92 mM), MgSO 4 (1.0 mM), Urea (5.0 mM), Fe-EDTA (0.04 mM), ZnSO 4 (0.6 µM), CaCl 2 .2H 2 O (0.75 mM), and micronutrients (CuSO 4 (0.62 µM), H 3 BO 3 (2.4 µM), MnSO 4 (0.9 µM) and Na 2 MoO 4 (0.6 µM)) [34]. The two levels of P were maintained as control and treatment i.e., normal P (NP) (250 µM) and low P (3.0 µM) using a KH 2 PO 4 salt solution [35]. The freshly prepared nutrient solution was used every alternate day and the pH was kept at 6.0 using 1.0 M HCL or 1.0 M KOH.

Trait Measurement
The 21-day old mungbean plants which were grown under NP and LP conditions were removed and dried in the oven at 60 • C until constant biomass was obtained for the measurement of various parameters. The total dry weight (TDW) was estimated as g/plant using a precision balance and the dried samples were ground to obtain a fine powder. The powder (0.1 g) was then digested using a diacid mixture (HNO 3 :HCLO 4 in 9:4 ratio) until a clear solution was obtained [36]. After filtering the solution, the P concentration (PC) (mg P/g dry weight) was estimated using a colorimetric method [37] and the total P uptake (TPU) (mg P/g dry weight) was calculated by multiplying the TDW and PC of the sample [38]. The P utilization efficiency (PUtE) (g dry weight/mg P) was calculated using the following formula under both NP and LP conditions [38].
PUtE (g dry weight/mg P) = Total dry weight/total P uptake by plant

Large Scale Genotyping of Association Panel Genotypes
A very large number of chromosome-based SNPs (55,634) were identified by our group in a previous study [20] using a GBS-based NGS assay. These chromosome-based SNPs were used for genotyping the P uptake and utilization-specific association mapping panel, which was constituted from 120 diverse mungbean genotypes. The structural and functional annotation of the data has been used by Reddy and co-workers [20].

Depiction of Linkage Disequilibrium, Phylogenetic Details, and Population Structure
To find the genome-wide linkage disequilibrium (LD) patterns (r 2 ) and LD decay, the genotyping data of 120 genotypes were analyzed using PLINK and TASSELv5.0 software [39]. For the preparation of an un-rooted neighbor-joining (NJ)-based phylogenetic tree, SNP genotyping data were used and analyzed using MEGA v6.0 [40] and Power-Marker v3.51 [41]. The principal component analysis (PCA) was performed using GAPIT, and population genetic structure was determined using STRUCTURE v2.3.4 following the details of Upadhyaya and co-workers [42]. The population structure analysis was performed and the population number or the K varied from 1 to 10, while the number of replications was kept at 20. The population number was then determined using the delta K value which is embedded in the second-order statistics of the STRUCTUREv2.3.4 software.

GWAS for P Use and P Utilization Efficiency-Associated Traits in Mungbean
GAPIT was used to integrate the phenotyping and genotyping (SNP) data, population structure coefficient (Q), kinship matrix (K), and PCA (P) using a mixed model (P+K, K, and Q+K) based MLM (mixed linear model) and CMLM (compressed mixed linear model) [42]. Afterward, p-value threshold corrections were performed using the false discovery rate (FDR cutoff ≤ 0.05) to improve the overall marker-trait association accuracy.

Digital Gene Expression Analysis and Validation for the Candidate Gene
In order to identify putative candidate genes, the CDS sequences of 65 associated protein coding genes were retrieved (https://legumeinfo.org/genomes/gbrowse/Vr1.0 (accessed on 12 March 2021)). A BLAST (nucleotide BLAST) search was performed against the Arabidopsis genome database (with a default parameter p value ≤ 1.0) (https://blast. ncbi.nlm.nih.gov/Blast.cgi (accessed on 12 March 2021)).To study the expression pattern of the candidate genes, digital gene expression analysis was performed. For this, the Arabidopsis orthologue of the identified candidate genes from the mungbean genome was used for the analysis. An online search tool viz., Expression Angler was then used to decode the expression pattern of the genes [43]. This platform identifies Arabidopsis genes with similar expression patterns. It calculates the correlation coefficients for all gene expression vectors as compared to an expression or to an expression pattern associated with an AGI ID or gene name that is specified [43]. The experiments were performed using several tissues, such as shoot apical meristem, shoot, root, xylem, etc. and at varying developmental stages. The BLAST search of the candidate genes was also performed against the Phaseolus vulgaris genome. The expression of orthologous genes were also checked using the Phaseolus vulgaris Gene Expression Atlas (Pv GEA: http://plantgrn.noble.org/PvGEA/ (accessed on 16 June 2021)) which facilitates functional genomic studies in the common bean (Table S4). The expression atlas has been developed using RNA-seq data [44].The SNPs present within three genes (VRADI01G04370, VRADI05G20860, and VRADI08G00070) were also validated using the Sanger sequencing approach (Table S3).

Phenotypic Variation for P Uptake and P Utilization Efficiency-Related Traits
An association mapping panel consisting of 120 diverse mungbean genotypes was evaluated for four traits viz., TDW, PC, TPU, and PUtE under both NP and LP conditions. Through analysis of variance studies (Table 1), highly significant variations could be identified among the studied genotypes for all four traits when studied under both P levels. The mean values of traits such as TDW, PC, and TPU were found to be significantly higher under NP when compared to LP conditions. However, a trait like PUtE showed significantly higher values under the LP condition. The coefficient of variation was recorded ranging from 9.73 to 18.50 and 13.10 to 36.47 under NP and LP conditions, respectively. The highest broad-sense heritability was recorded for TDW under both NP (0.84) and LP (0.80) conditions. The relative values (trait ratio between LP and NP conditions) of three traits viz., TDW, PC, and TPU, were recorded as less than 1.00, whereas the relative value for PUtE was recorded as more than 1.00. The frequency distribution pattern of these four traits showed a normal distribution pattern which means that they are complex and quantitative ( Figure 1). (0.84) and LP (0.80) conditions. The relative values (trait ratio between LP and NP condi tions) of three traits viz., TDW, PC, and TPU, were recorded as less than 1.00, whereas the relative value for PUtE was recorded as more than 1.00. The frequency distribution pat tern of these four traits showed a normal distribution pattern which means that they are complex and quantitative ( Figure 1).  Frequency distribution of variation for total dry weight (TDW), P concentration (PC), total P uptake (TPU) and P utilization efficiency (PUtE) in 120 diverse mungbean genotypes under normal P (NP) and low P (LP) conditions. P, phosphorus; DW, dry weight.

Linkage Disequilibrium, Phylogenetic Tree and Population Structure of AM Panel Genotypes
A total of 55634 genome-wide and chromosome-based SNPs were used to construc the unrooted NJ phylogenetic tree; to find the LD estimates (r 2 ), and also to study the decay among the genotypes of the association panel. The LD decay was estimated by pooling the r 2 value across the eleven mungbean chromosomes. Subsequently, the aver age r 2 was plotted against equal physical intervals (50 kb). The analysis showed a higher Figure 1. Frequency distribution of variation for total dry weight (TDW), P concentration (PC), total P uptake (TPU) and P utilization efficiency (PUtE) in 120 diverse mungbean genotypes under normal P (NP) and low P (LP) conditions. P, phosphorus; DW, dry weight.

Linkage Disequilibrium, Phylogenetic Tree and Population Structure of AM Panel Genotypes
A total of 55,634 genome-wide and chromosome-based SNPs were used to construct the un-rooted NJ phylogenetic tree; to find the LD estimates (r 2 ), and also to study the decay among the genotypes of the association panel. The LD decay was estimated by pooling the r 2 value across the eleven mungbean chromosomes. Subsequently, the average r 2 was plotted against equal physical intervals (50 kb). The analysis showed a higher LD estimate (r 2 :0.72) and a comparatively less extensive decay (r 2 decreased to half of its maximum value: 0.31) at around a 70-100 kb distance in mungbean chromosomes ( Figure 2). The LD pattern recorded an increase and then an LD decay (r 2 ≥ 0.3) which followed a consistent pattern with increasing physical distance (Kb) of the SNPs [20,45]. LD estimate (r 2 :0.72) and a comparatively less extensive decay (r 2 decreased to half of its maximum value: 0.31) at around a 70-100 kb distance in mungbean chromosomes (Figure 2). The LD pattern recorded an increase and then an LD decay (r 2 ≥0.3) which followed a consistent pattern with increasing physical distance (Kb) of the SNPs [20,45]. The unrooted phylogenetic tree constructed using the NJ method depicted a clear grouping among the mungbean genotypes of the AM panel ( Figure 3). The genotypes were grouped into five major clusters. Furthermore, the population genetic structure of the mungbean genotypes was depicted by employing STRUCTUREv2.3.4 software. The delta K value showed its peak at five ( Figure 4A), confirming the grouping of 120 genotypes into five genetically distinct population groups (POP I-V) ( Figure 4B. This grouping was further confirmed by PCA, along with population structure analysis and construction of a phylogenetic tree ( Figure 4C).  The un-rooted phylogenetic tree constructed using the NJ method depicted a clear grouping among the mungbean genotypes of the AM panel ( Figure 3). The genotypes were grouped into five major clusters. Furthermore, the population genetic structure of the mungbean genotypes was depicted by employing STRUCTUREv2.3.4 software. The delta K value showed its peak at five ( Figure 4A), confirming the grouping of 120 genotypes into five genetically distinct population groups (POP I-V) ( Figure 4B). This grouping was further confirmed by PCA, along with population structure analysis and construction of a phylogenetic tree ( Figure 4C). LD estimate (r 2 :0.72) and a comparatively less extensive decay (r 2 decreased to half of its maximum value: 0.31) at around a 70-100 kb distance in mungbean chromosomes (Figure 2). The LD pattern recorded an increase and then an LD decay (r 2 ≥0.3) which followed a consistent pattern with increasing physical distance (Kb) of the SNPs [20,45]. The unrooted phylogenetic tree constructed using the NJ method depicted a clear grouping among the mungbean genotypes of the AM panel ( Figure 3). The genotypes were grouped into five major clusters. Furthermore, the population genetic structure of the mungbean genotypes was depicted by employing STRUCTUREv2.3.4 software. The delta K value showed its peak at five ( Figure 4A), confirming the grouping of 120 genotypes into five genetically distinct population groups (POP I-V) ( Figure 4B. This grouping was further confirmed by PCA, along with population structure analysis and construction of a phylogenetic tree ( Figure 4C).

Candidate Gene Identification for PUse and PUtilization Traits in Mungbean Using Genome-Wide GBS
The genome-wide genotyping and phenotyping data of the AM panel consisting of 120 mungbean genotypes were integrated for the identification of a population structure ancestry coefficient (Q) and kinship matrix using MLM and CMLM models. The threshold Bonferroni correction value of −log (p) = 3.5 was used as a cutoff to identify the significant SNPs associated with four traits viz., TDW, PC, TPU, and PUtE under NP, LP, and LP/NP conditions. The association study could identify a nearly similar number of SNPs by both MLM (125) and CMLM (123) models (Table 2, Figure 5A,B). Further, 116 common SNPs were found in both MLM and CMLM models (Supplementary Table S2). TPU showed a maximum number of associated SNPs under the NP condition (10) followed by LP (7) and LP/NP (5) conditions. For PUtE, maximum SNPs (23) were observed under the LP condition followed by the NP (10) and LP/NP (9) conditions. Similarly, for TDW and PC, maximum SNPs were found associated under LP (12) and NP (18) Table S2). Among these, 62 are in the intergenic region, 23 are in the coding region, 24 are in the upstream region and 7 are found in the downstream region. These genes play diverse roles which include ion transport, metabolism, stress, and development. Interestingly, a large number of genes were also found to have a role in nutrient transport activity.  Figure 5A,B). Further, 116 common SNPs were found in both MLM and CMLM models (Supplementary Table S2). TPU showed a maximum number of associated SNPs under the NP condition (10) followed by LP (7) and LP/NP (5) conditions. For PUtE, maximum SNPs (23) were observed under the LP condition followed by the NP (10) and LP/NP (9) conditions. Similarly, for TDW and PC, maximum SNPs were found associated under LP (12) and NP (18) Table S2). Among these, 62 are in the intergenic region, 23 are in the coding region, 24 are in the upstream region and 7 are found in the downstream region. These genes play diverse roles which include ion transport, metabolism, stress, and development. Interestingly, a large number of genes were also found to have a role in nutrient transport activity.

Candidate Gene Identification for PUse and PUtilization Traits in Mungbean Using Genome-Wide GBS
The genome-wide genotyping and phenotyping data of the AM panel consisting of 121 mungbean genotypes were integrated for the identification of a population structure ancestry coefficient (Q) and kinship matrix using MLM and CMLM models. The threshold Bonferroni correction value of -log (p) = 3.5 was used as a cutoff to identify the significant SNPs associated with four traits viz., TDW, PC, TPU, and PUtE under NP, LP, and LP/NP conditions. The association study could identify a nearly similar number of SNPs by both MLM (125) and CMLM (123) models (Table 2, Figure 5A,B). Further, 116 common SNPs were found in both MLM and CMLM models (Supplementary Table S2). TPU showed a maximum number of associated SNPs under the NP condition (10) fol-   Table 2. Genome-wide association study for total dry weight (TDW), P concentration (PC), total P uptake (TPU), and P utilization efficiency (PUtE) studied using a mixed linear model (MLM) and compressed mixed linear model (CMLM).  Table S2). Among these, 62 are in the intergenic region, 23 are in the coding region, 24 are in the upstream region and 7 are found in the downstream region. These genes play diverse roles which include ion transport, metabolism, stress, and development. Interestingly, a large number of genes were also found to have a role in nutrient transport activity.

Delineation of Putative Candidate Genes for P Uptake and P Utilization Efficiency Traits in Mungbean
The BLAST search resulted in the identification of 16 genes with more than an 80% identity index to the Arabidopsis genes, having varied roles in nutrient uptake and stressrelated pathways (Table 3). These genes were concluded to be putative candidate genes regulating the P uptake and P utilization efficiency in mungbean crops. The digital expression pattern of 16 genes revealed that six genes, namely VRADI01G04370 (AT3G18480), VRADI05G20860 (AT1G51700), VRADI06G12490 (AT3G13560), VRADI08G20910 (AT1G49975), VRADI08G00070 (AT4G02650) and VRADI09G09030 (AT3G44700) have high expression in different plant parts including root, shoot apical meristem, leaf, etc. The gene VRADI08G20910 is a ubiquitin-conjugating enzyme E3 NLA gene and is involved in PHT1 (P transporter) ubiquitination. The gene VRADI06G12490 has been reported to be from the SPX gene family and is involved in P signaling and homeostasis by negatively regulating the activity of PHR (Phosphate starvation response regulator). The gene VRADI09G09030 is a kinase gene that is underlying in the PUP1(Phosphorus uptake 1) QTL (Phosphorus-starvation tolerance 1 (PSTOL1)) region and is involved in promoting P uptake by enhancing early root growth in rice. The VRADI08G00070 gene is a phosphatidylinositol binding clathrin assembly protein with overlapping functions in recycling ANX (ANXUR receptor-like kinases) proteins to the pollen tube membrane. Therefore, we conclude that these genes may play a crucial role as potential candidate genes in regulating the P uptake and P utilization efficiency in mungbean crops.

Discussion
To improve the PUE, both TPU and PUtE are required to be exploited in crop plants [16,62]. The significant contribution of total biomass, P concentration, and total P uptake towards PUE was reported in rice [63] and wheat [64]. The mean values of TDW, PC, and TPU were found to be lower under LP conditions, whereas the mean value of PUtE was found to be higher under LP conditions, which is attributed to the higher TPU of the genotypes. The results are in good agreement with the earlier reports forrice [65] and wheat [66]. The highly significant interaction between genotype and P treatment indicated the significant effect of two P regimes on the studied genotypes for the tested traits. The studied genotypes showed the presence of a significant amount of variability for the investigated traits. Silva et al. [67] in common bean and Wang et al. [68] in maize, reported the presence of significantly higher genetic variability with high heritability as reliable selection criteria for PUE improvement in crop plants.
It is imperative to decode the LD pattern in order to dissect the genomic landscape of a crop species. In mungbean, limited effort has been made to decipher the high-resolution LD pattern [20,45]. In a previous study, the LD decay pattern was found to be between 50-100 kb in mungbean [20]. In this study of 120 mungbean accessions, a less extensive LD decay was observed which was around 80 kb (r 2 :0.31) for mungbean. In soybean, the extent of LD decay was found to be at a physical distance of around 27 kb, 83 kb, and 133 kb for wild soybean, landraces, and improved cultivars, respectively [69]. In chickpea, a larger decay was reported which was nearly 150-200 kb [39,[70][71][72][73]; whereas, in common bean and pigeonpea, the LD decay extent was reported as 50-70 kb [74] and 70 kb [75], respectively, which is comparable to mungbean [20]. In rice, the extent of LD decay was found to be nearly 75-150 kb, which is slightly higher than that of mungbean [76]. However, in wheat, the extent of LD was 2.1 Mb for subgenome A, 4.2 Mb for subgenome D, and 5.9 Mb for subgenome B, which is very much higher when compared to mungbean [77]. Population structure analysis revealed the presence of five subpopulations, which were further confirmed by PCA and un-rooted phylogenetic tree formation.
The GWAS revealed 116 novel SNPs that were shared by both MLM and CMLM models and are found to be associated with four traits viz., TDW, PC, TPU, and PUtE. Further, 61 associated protein-coding genes were also found to have diverse roles in stress and other metabolic pathways. This study could identify genes containing zinc finger domain, bHLH domain, metallophos domain, and SPX domains as also identified in our previous study [20]. These domain-containing genes have been reported to play a diverse role in controlling the root architecture as well as adaptation to P starvation stress [72][73][74][75][76]. Of these 61 genes, 16 were found to have a crucial role in the imposition of phosphorous stress tolerance in plants. Thus, it may be concluded that these genes might be the putative candidate genes having their role in the imposition of P stress tolerance in plants.
The VRADI07G06240 is a zinc finger CCCH domain-containing protein that is found to be expressed under short-term P-deprivation conditions in soybean leaves [77]. Whereas the VRADI08G11310 gene encodes a bHLH domain-containing protein, that may be involved in root hair formation and thereby increasing the P uptake under P starvation conditions [78]. Also, VRADI05G06040 is a metallophos domain-containing protein that seems to play a role in P uptake and utilization by releasing inorganic P [79]. Furthermore, VRADI08G10870, VRADI07G24790, and VRADI03G02620 are found to be associated with the root development [80,81]. The VRADI08G20910 is a ubiquitin-conjugating enzyme coding gene that is known to cause degradation of PHT1 gene function. PHT1 is a membrane protein that helps in the translocation of P in and out of the cell under both P-deficient and sufficient conditions. It is speculated that when PHT1 gets mutated, the normal process of transport or degradation of P gets hampered and ultimately the process of influx and efflux of P gets interrupted ( Figure 6).
VRADI06G12490 is an SPX domain-containing gene family that is involved in the negative regulation of the PHR (Phosphate starvation response regulator) gene by regulating the expression of the PSI (Phosphate starvation-induced) gene [53]. When this gene gets mutated, then its binding to the SPX domain gets hampered, which in turn causes interrupted expression of the PSI gene and ultimately disruption of P uptake (Figure 7). VRADI06G14930 gene is involved in the reduction of protein catabolism under P-deficient conditions; while the VRADI06G03940 gene is involved in the transport of heavy metals and detoxification in plant cells. Two more genes viz., VRADI05G21880 and VRADI06G13450 are known to function under P-deficient conditions. Besides, the gene VRADI08G00070 is known to recycle ANXUR proteins to the pollen tube membrane and is expressed under nutrient-deficient conditions. VRADI01G04370 and VRADI05G20860 genes are involved in the abscisic acid-deficient stress response [54][55][56][57].
gronomy 2021, 11, x FOR PEER REVIEW 13 of 1 P-deficient and sufficient conditions. It is speculated that when PHT1 gets mutated, th normal process of transport or degradation of P gets hampered and ultimately the pro cess of influx and efflux of P gets interrupted ( Figure 6). VRADI06G12490 is an SPX domain-containing gene family that is involved in th negative regulation of the PHR (Phosphate starvation response regulator) gene by regu lating the expression of the PSI (Phosphate starvation-induced) gene [53]. When this gen gets mutated, then its binding to the SPX domain gets hampered, which in turn cause interrupted expression of the PSI gene and ultimately disruption of Puptake (Figure 7) VRADI06G14930 gene is involved in the reduction of protein catabolism unde P-deficient conditions; while the VRADI06G03940 gene is involved in the transport o heavy metals and detoxification in plant cells. Two more genes viz., VRADI05G21880and VRADI06G13450are known to function under P-deficient conditions. Besides, the gen VRADI08G00070 is known to recycle ANXUR proteins to the pollen tube membrane and is expressed under nutrient-deficient conditions. VRADI01G04370 and VRA DI05G20860genes are involved in the abscisic acid-deficient stress response [54][55][56][57].  VRADI06G12490 is an SPX domain-containing gene family that is involved in th negative regulation of the PHR (Phosphate starvation response regulator) gene by regu lating the expression of the PSI (Phosphate starvation-induced) gene [53]. When this gen gets mutated, then its binding to the SPX domain gets hampered, which in turn cause interrupted expression of the PSI gene and ultimately disruption of Puptake ( Figure 7 VRADI06G14930 gene is involved in the reduction of protein catabolism unde P-deficient conditions; while the VRADI06G03940 gene is involved in the transport o heavy metals and detoxification in plant cells. Two more genes viz., VRADI05G21880an VRADI06G13450are known to function under P-deficient conditions. Besides, the gen VRADI08G00070 is known to recycle ANXUR proteins to the pollen tube membrane an is expressed under nutrient-deficient conditions. VRADI01G04370 and VRA DI05G20860genes are involved in the abscisic acid-deficient stress response [54][55][56][57]. . CMLM based Manhattan plots depicting the significant association of SNP markers with total dry weight (TDW), P concentration (PC), total P uptake (TPU) and P utilization efficiency (PUtE) under NP, LP and LP/NP conditions. Furthermore, digital expression analysis using 16 genes revealed significantly higher expression of six genes namely, VRADI01G04370, VRADI05G20860, VRADI06G12490, VRADI08G20910, VRADI08G00070 and VRADI09G09030 in the root, shoot apical meristem, and leaf tissues (Figures S1-S12). The digital gene expression analysis was also carried out for the orthologous gene of Phaseolus vulgaris. The result showed that most of the gene was highly expressed in leaf, pod and root (Table S3). Therefore, it can be concluded that these genes might act as potential candidate genes in regulating the P uptake and P utilization efficiency in mungbean plants. We found a smaller degree of expression in the pollen tube as well as different kinds of tissues in different developmental stages, viz. ovary and hypocotyls. Phosphorous is an important component of biological macromolecules and after uptake, it is sequestered in different components of plant tissues throughout the plant body. Therefore, it may be possible these genes also play an important role in its sequestration in those tissues. In our previous study [20] we have reported 13 candidate genes regulating PUE traits in mungbean, whereas this study identified 16 novel genes controlling P uptake and P utilization efficiency traits in mungbean. Chromosome-based analysis revealed that out of 29 genes, 7 (~25%) were found located on chromosome 6, while there were 4 genes on chromosome 8 and 3 genes on chromosome 7 (Figure 8). This result suggests that chromosome six may be the hotspot region harboring several major QTLs/genes regulating P uptake and P utilization efficiency traits in mungbean. However, further in-depth studies are required to reconfirm and strengthen our results.

Conclusions
This study employs a genome-wide GBS approach to identify the significant SNPs and candidate genes associated with P uptake and P utilization efficiency traits in mungbean. It resulted in the identification of 116 SNPs associated with the 61 protein-coding genes, of which 16 genes were found to play a crucial role in enhancing the P uptake and P utilization efficiency in several crop plants, including Arabidopsis. A total of seven genes out of twenty-nine(including our previous study) were found to be located on chromosome six. This result suggests that chromosome six may be the hotspot region harboring several major QTLs/genes regulating P uptake and P utilization efficiency traits in the mungbean. Therefore, chromosome six can be enriched and targeted for genomics assisted breeding to improve the phosphorous deficiency response in mungbean. The LD analysis revealed a comparatively short LD decay in mungbean. This study leads towards deciphering the genetic complexity for PUpE and PUtE in mungbean. The genetic and genomic resources generated in this study can be exploited to develop phosphorous stress responsive cultivars which can overcome P deficiency in mungbean.