Abstract
Chickpea—the second most important grain legume worldwide—is cultivated mainly on marginal soils. Phosphorus (P) deficiency often restricts chickpea yields. Understanding the genetics of traits encoding P-acquisition efficiency and P-use efficiency will help develop strategies to reduce P-fertilizer application. A genome-wide association mapping approach was used to determine loci and genes associated with root architecture, root traits associated with P-acquisition efficiency and P-use efficiency, and any associated proxy traits. Using three statistical models—a generalized linear model (GLM), a mixed linear model (MLM), and a fixed and random model circulating probability unification (FarmCPU) —10, 51, and 40 marker-trait associations (MTAs), respectively were identified. A single nucleotide polymorphism (SNP) locus (Ca1_12310101) on Ca1 associated with three traits, i.e., physiological P-use efficiency, shoot dry weight, and shoot P content was identified. Genes related to shoot P concentration (NAD kinase 2, dynamin-related protein 1C), physiological P-use efficiency (fasciclin-like arabinogalactan protein), specific root length (4-coumarate–CoA ligase 1) and manganese concentration in mature leaves (ABC1 family protein) were identified. The MTAs and novel genes identified in this study can be used to improve P-use efficiency in chickpea.
Introduction
Phosphorus (P) is an essential nutrient for crop production. Using spatially explicit global maps for more than 100 crops, reported that 29% of the global cropland area is P deficient, while 71% has surplus P. For example, 42% of field soil in India is low in P, 38% is medium and 20% is high (). Excessive P fertilizer and manure application in industrialized countries have led to low P-use efficiency (PUE), with surplus P retained in soil as residual P () or lost to the environment where it is causing significant water quality problems. Globally 51–86% more P input will be needed for sustainable crop production by 2050 (), unless we work toward more P-efficient crops and cropping systems ().
Chickpea (Cicer arietinum L.), is one of the most important grain legumes cultivated by smallholder farmers in more than 50 countries. Advances in chickpea genomics in the last decade have made large-scale genomic resources available to the chickpea research community including molecular markers (; ), genetic maps (), genome sequences (), and resequencing of several germplasm lines (, ; ). These resources have improved our understanding of both abiotic () and biotic stress tolerance in chickpea () and enabled fine mapping of traits (; ). Furthermore, the resources have been successfully used to develop new varieties with enhanced tolerance or resistance (see ; ; ; ). Twelve chickpea genotypes with well-known responses under drought and irrigation were evaluated for profuse root length density (RLD) in surface soil and root dry weight (RDW) and root:shoot ratio (RSR) in deeper soil layers (). This study revealed that drought stress increased RLD below 300 mm soil depth, deep RDW, and RSR, but decreased root diameter. Moisture-conservation practices and optimum P levels to enhance PUE were reported recently in chickpea grown in vertisols in central India (). In soybean, root length is positively correlated with P accumulation in well-watered and water-stressed conditions (). Phenotypic plasticity and genetic variability in root architectural traits of chickpea and their role in drought tolerance using a novel semi-hydroponic system have been reported (). A recent study revealed that root system plasticity affects P acquisition efficiency, PUE, and photosynthetic PUE in 266 chickpea genotypes (). Root exudates like carboxylate enable P acquisition from the soils which are low in available P. Manganese concentration in mature chickpea leaves is positively correlated with the amount of rhizosheath carboxylates, offering an easily measurable proxy for assessing rhizosheath carboxylates in 100 chickpea genotypes grown under low P availability (). In the context of global climate change, summarized the factors affecting PUE and enhancing P-acquisition efficiency in legumes, and key areas for future research. Understanding the genetics of these traits, identification of genomic regions, molecular markers and or marker trait associations for PUE efficiency related traits will help improving these traits through marker-assisted selection or genomics-assisted breeding.
Various type of molecular markers have been used to establish marker-trait associations (MTA) in chickpea. For instance, 1,072 Diversity Arrays Technology (DArTs), 651 single nucleotide polymorphisms (SNPs), 113 gene-based SNPs, and 36 simple sequence repeats (SSRs) were used to establish 312 MTAs for drought- and heat-tolerance related traits (). SNPs were used extensively for association studies in chickpea due to their abundance and amenability for high-throughput genotyping (; ; ; ; ).
In legumes, genome-wide association studies (GWAS) are gaining momentum. A recent study reported two SNP markers tightly linked to seed iron (Fe) and one to seed zinc (Zn) concentration in lentils (Lens culinaris Medik.) (). Another study reported 159 quantitative trait nucleotides (QTNs) and 52 candidate genes associated with the photosynthetic response to low-P stress in soybean [Glycine max (L.) Merr.] (). In chickpea, association mapping between SNP markers and seed copper (Cu), P, and potassium (K) concentrations identified eight SNPs associated with variation in three nutrients in more than two environments (). Similarly, seed mineral concentration in pea (Pisum sativum L.; ), Fe chlorosis in soybean; ; ), and Fe bioavailability in cooked dry beans (Phaseolus vulgaris L.; ) have been mapped. Novel genes involved in the accumulation of P in Lotus japonicas have been reported using GWAS analysis ().
In view of above, this study was conducted to undertake GWAS analysis to identify MTAs for (i) root architectural traits evaluated in a high-throughput semi-hydroponic root phenotyping platform, and (ii) root morphological and physiological traits related to P-acquisition efficiency and P−use efficiency under low P supply. This is the first study that reports genomic regions associated with above mentioned traits by using three different models, namely fixed and random model circulating probability unification (FarmCPU), mixed linear model (MLM), and generalized linear model (GLM) in the GAPIT-R package. The reference genome has been used to identify the candidate genes in the identified MTAs associated with above traits.
Materials and Methods
Germplasm Lines and Phenotyping
The chickpea reference set () comprising of 300 diverse accessions (267 landraces, 13 advanced lines and cultivars, 7 wild Cicer accessions, and 13 accessions with unknown biological status) was used for phenotyping datasets as following: (i) 270 genotypes of the reference set were evaluated for 30 root architectural traits in a semi-hydroponic phenotyping system (); (ii) Two hundred and sixty-six genotypes (including 255 from the chickpea reference set along with 11 Australian chickpea cultivars (Ambar, Almaz, Neelam, Genesis 079, Genesis 090, Genesis 509, Genesis 836, Genesis Kalkee, PBA Boundary, PBA Slasher and PBA Striker) were evaluated for P-acquisition efficiency and P−use efficiency with P supplied as insoluble FePO4 (); and (iii) a selected subset of 100 chickpea genotypes of the reference set—showing visual differences in plant size or leaf symptoms of P deficiency—phenotyped for shoot/root morphological and physiological traits to understand the relative roles of root morphology and physiology in P-acquisition efficiency (, ).
Genotyping and Determining Population Structure
The SNP dataset, based on whole-genome resequencing data on the reference dataset was filtered for missing values (≥20%) and minor-allele frequency <5% using vcftools and imputed by BEAGLE-5.0 (; ). As the number of genotypes were different in data sets, Population structure was separately determined for the chickpea reference set and 91 genotypes (from the 100 selected subset of 199 genotypes) using ADMIXTURE v1.3.0 ().
Genome-Wide Association Analysis
GWAS for Root Traits and Seven P-Acquisition Efficiency Traits
To identify significant MTAs and avoid spurious associations of the 270 genotypes mentioned above, 233 genotypes (with phenotyping data and genotyping data) were considered for GWAS analysis. GWAS analysis was performed using 698,183 SNPs (SNP calls obtained based on aligning 233 genotypes to reference genome CDC Frontier) and phenotyping data for 37 traits (i.e., 30 root traits and seven P-acquisition efficiency traits). MTAs were determined using three models, namely fixed and random model circulating probability unification (FarmCPU), mixed linear model (MLM), and generalized linear model (GLM) in the GAPIT-R package ().
GWAS for Biochemical Traits/Proxy Traits for PUE
A total of 706,865 SNPs and phenotyping data generated on 91 chickpea germplasm lines was for analysis in GAPIT-R using FarmCPU and MLM models to determine MTAs for biochemical traits.
The Bonferroni correction threshold of 7.07E-08 was used to avoid spurious associations. The genes involving significant SNP markers were aligned against the NCBI non-redundant (nr) protein database taxon Viridiplantae using BLASTX, to obtain functional annotations. GO and KEGG (Kyoto Encyclopedia of Genes and Genomes) pathway identification were conducted on these sequences in the KEGG pathways in-built in BLAST2GO. The SNPEff- 4.3T open source program was used for variant annotation and prediction of significant SNP effects.
Results
Population Structure and Genome-Wide Association Study (GWAS)
Using genome-wide SNP data for 233 chickpea genotypes (with both genotypic and phenotypic data), three subpopulations were identified (Supplementary Figures 1A,B) using ADMIXTURE. Similarly, in the smaller subset of 91 genotypes, three subpopulations were identified (Supplementary Figures 1C,D). Three statistical models—GLM, MLM, and FarmCPU—enabled identification of 10, 51, and 40 MTAs, respectively (Table 1 and Supplementary Table 1) after applying Bonferroni and FDR corrections. Spurious MTAs were excluded by examining Q/Q plots. Forty-two of the MTAs identified by the GLM and MLM models were robust (>10% of the phenotypic variation explained); the phenotypic variation explained by the MTAs identified in the FarmCPU model was not computed. The GLM model identified MTAs for four traits (root growth rate, root mass density, specific root length, and shoot P concentration); this low number of traits could be due to the many spurious MTAs (Supplementary Figure 2). The FarmCPU model identified MTAs for the most traits, followed by the MLM model.
TABLE 1
| Trait | Code | GLM | MLM | FarmCPU | |||||
| No. of MTAs | p-value | R2 | No. of MTAs | p-value | R2 | No. of MTAs | p-value | ||
| Average branch length | (ABL; cm per branch) | – | – | – | 2 | 5.92 × 10–8 to 3.69 × 10–8 | 5.01–8.24 | 1 | 5.02 × 10–8 |
| Branch density | (BD; cm–1 taproot) | – | – | – | – | – | – | 2 | 6.41 × 10–10 to 2.67 × 10–10 |
| Branch intensity | (BI; cm–1 root) | – | – | – | – | – | – | 1 | 7.51 × 10–9 |
| Root diameter s2 | (RD_s2; mm) | – | – | – | – | – | – | 1 | 5.14 × 10–9 |
| Subsoil root diameter | (RD_sub; mm) | – | – | – | – | – | – | 1 | 7.2 × 10–14 |
| Topsoil root diameter | (RD_top; mm) | – | – | – | – | – | – | 1 | 7.46 × 10–12 |
| Root growth rate | (RGR; cm d–1) | 1 | 2.2 × 10–8 | 13.13 | 12 | 5.62 × 10–15 to 1.9 × 10–8 | 9.2–11.78 | – | – |
| Root length ratio | (RLR_top/sub) | – | – | – | – | – | – | 1 | 1.06 × 10–8 |
| Root length s2 | (RL_s2; cm) | – | – | – | – | – | – | 1 | 1.59 × 10–8 |
| Root mass | (RM; mg) | – | – | – | – | – | – | 1 | 3.86 × 10–8 |
| Root mass ratio | RMR | – | – | – | 6 | 4.4 × 10–25 to 1.73 × 10–8 | 20.86–30.12 | 1 | 3.8 × 10–8 |
| Root_2 | Root diameter (mm) | – | – | – | 1 | 1.77 × 10–12 | 7.52 | – | – |
| Root tissue density | (RTD; mg cm–3) | 2 | 2.88 × 10–9 to 6.73 × 10–8 | 14.23 to 16.95 | – | – | – | 2 | 7.18 × 10–15 to 4.58 × 10–10 |
| Shoot dry weight | SDW (mg) | – | – | – | – | – | – | 1 | 6.31 × 10–9 |
| Specific root length | (SRL; m g–1 dry mass) | 3 | 2.48 × 10–9 to 5.94 × 10–8 | 17.2 to 20.31 | 2 | 1.88 × 10–10 to 1.44 × 10–9 | 17.63 | 3 | 2.99 × 10–11 to 4.06 × 10–9 |
| Taproot length zone1 | (TRL_z1; cm) | – | – | – | 1 | 7.15 × 10–8 | 13.59 | 1 | 7.15 × 10–8 |
| Taproot length zone2 | (TRL_z2; cm) | – | – | – | 1 | 1.48 × 10–11 | 9.64 | 1 | 1.48 × 10–11 |
| WUE | WUE (A/gs) | – | – | – | – | – | – | 2 | 2.80 × 10–11 to 1.85 × 10–8 |
| Total rhizo dry soil (g plant-1) | Total rhizo (g plant–1) | – | – | – | 6 | 5.66 × 10–17 to 1.34 × 10–11 | 5–17.63 | – | – |
| Phosphorous-utilization efficiency | – | – | – | 13 | 4.45 × 10–18 to 8.41 × 10–15 | 20.23–26.16 | 1 | 2.13 × 10–9 | |
| Physiological P-use efficiency | PPUE (μmol g–1 P s–1) | 1 | 5.46519E-08 | 11.95 | – | – | – | – | – |
| Shoot phosphorous concentration | mg g–1 | 3 | 2.85 × 10–9 to 7.65 × 10–8 | 27.67 to 40.33 | 2 | 5.23 × 10–9 to 2.85 × 10–9 | 27.69–32.97 | – | – |
| Shoot P content | Shoot P content | – | – | – | 5 | 8.2 × 10–9 to 7.7 × 10–8 | 39.41.3 | 1 | 1.68 × 10–13 |
| Carboxylate_2 | Carboxylate conc (μmol g–1 root DW) | – | – | – | – | – | – | 1 | 5.93 × 10–8 |
| Ci | Ci | – | – | – | – | – | – | 1 | 8.52 × 10–10 |
| Citric_2 | Citric (μmol plant–1) | – | – | – | – | – | – | 1 | 1.42 × 10–9 |
| Malonic (μmol plant–1) | Malonic (μmol plant–1) | – | – | – | – | – | – | 1 | 4.82 × 10–10 |
| Mn concentration in mature leaves | Mn_ML | – | – | – | – | – | – | 7 | 6.08 × 10–21 |
| P_ML | P_ML (mg g–1) | – | – | – | – | – | – | 1 | 6.08 × 10–21 |
| P_Pn | Pn_area (μmol m–2 s–1) | – | – | – | – | – | – | 1 | 2.6 × 10–8 |
| Pn_mass | Pn_mass (μmol g–1 s–1) | – | – | – | – | – | – | 1 | 3.57 × 10–9 |
| Rhizo-pH | – | – | – | – | – | – | 2 | 1.72 × 10–15 to 2.24 × 10–8 | |
| Specific rhizosheath weight | (g g–1 root DW) | – | – | – | – | – | – | 1 | 2.97 × 10–8 |
| Total MTAs | 10 | 51 | 40 | ||||||
Summary of marker trait associations identified using three statistical models namely GLM, MLM and FarmCPU.
GWAS Signals for Root Architectural Traits
Fifty-seven MTAs were identified for 19 root architectural traits phenotyped in our earlier studies (; ; Table 1). The FarmCPU model identified one significant MTA each in root depth zone 1 (TRL_z1; cm) and root depth zone 2 (TRL_z1; cm) on Ca3 and Ca6 for taproot length (i.e., root depth, z1 + z2), one significant MTA (Ca3_26114159) for root mass (RM, mg), two significant MTAs (Ca4_29694614 and Ca6_57970784) for branch density (BD; cm–1 taproot), and one significant MTA (Ca2_6340118) for branch intensity (BI; cm–1 root). No significant MTAs were detected for the above root traits in the GLM and MLM methods. Nevertheless, the MLM model identified six significant MTAs and FarmCPU identified one MTA for root mass ratio (RMR), and the MLM and FarmCPU models identified two significant MTAs (one each) for root length ratio (RLR_top/sub) and three significant MTAs (two and one, respectively) for average branch length (ABL; cm per branch). The FarmCPU model identified one significant MTA each for root diameter s2 (RD_s2; mm), subsoil root diameter (RD_sub; mm), and topsoil root diameter (RD_top; mm) on Ca6, Ca1, and Ca2, respectively. The FarmCPU model also identified two significant MTAs for water-use efficiency (WUE; A/gs) on Ca6 and Ca8. Based on the physical position of the SNP loci (MTAs) associated with different traits on Ca4, none of the MTAs was mapped in “QTL-hotspot_a” or “QTL-hotspot_b” that harbors several drought-tolerance-related root traits (). The MTAs revealed key root traits for efficiently acquiring soil resources and adapting to drought and other abiotic stresses. Of the 13 significant MTAs (12 from MLM and one from GLM) identified for RGR (cm d–1), six (46.2%) were identified on Ca3 and three (23.1%) on Ca6. Six MTAs identified on Ca6 explained 3.6–10.0% of the phenotypic variation associated with five genes (Ca_08259, Ca_01151, Ca_01152, Ca_01156, and Ca_16553). Both statistical models (GLM and MLM) detected one significant MTA (Ca1_4716136), which explained higher phenotypic variation (11.8–13.1%) than the other MTAs (Supplementary Table 1). Eight significant MTAs were identified for specific root length (SRL, cm) in all three models; of these, one MTA on Ca7 (Ca7_3606123), explaining about 17.8% of the phenotypic variation, was consistent in all three models used for GWAS analysis. The MTA was associated with the Ca_03107 gene that encodes pectinesterase/pectinesterase inhibitor PPE8B. Six MTAs for total rhizosheath dry soil (g plant–1) were identified on Ca8 (3), Ca5 (1), Ca2 (1), and Ca1 (1).
GWAS Signals for P-Acquisition Efficiency and PUE-Related Traits
A total of 10 significant MTAs (three based on GLM and seven based on MLM) for shoot P content—six on Ca4, two on Ca6 and one each on Ca7 and Ca2 were identified (Figure 1A). Two significant MTAs on Ca4 (Ca4_38518152 and Ca4_8269508) were identified with the GLM and MLM models (Supplementary Table 1). The 15 MTAs identified—14 from MLM and one from FarmCPU—explained 5–21% of the phenotypic variation. Five of the MTAs were identified on Ca6 followed by Ca7 (4), Ca4 (2), Ca2 (2), Ca1 (1), and C3 (1) (Figure 1B). A SNP locus, Ca7_33808891, associated with physiological P-Use efficiency present in gene Ca_16189 on Ca7 explained 20.23% PVE (Figure 1C; Supplementary Table 1). A SNP locus, Ca4_38518152, associated with Shoot phosphorus content present in gene Ca_13110 on Ca4 explained 31.5% PVE (Figure 1D; Supplementary Table 1). In the case of physiological PUE one significant MTA (Ca1_12310101) on Ca1 explaining 12.0% of the phenotypic variation was identified. A single nucleotide polymorphism locus (Ca1_12310101) on Ca1 associated with three traits, i.e., physiological P-use efficiency, shoot dry weight, and shoot P content was identified (Figure 2).
FIGURE 1
FIGURE 2
GWAS Signals for Proxy Traits
One significant MTA for P concentration in mature leaves (P_ML; mg g–1) was identified on Ca4. Seven significant MTAs for Mn concentration in mature leaves were identified on Ca2 (3), Ca4 (2), and Ca7 (Figure 3). The SNP loci on Ca7 (Ca7_32383349) and Ca4 (Ca4_1791932) were associated with manganese concentration in mature leaves with two different p-values (Figure 3), while three SNPs associated with the Mn_ML on Ca2 (Ca2_7561143, Ca2_866639, and Ca2_359984). One MTA each for citric (μmol plant–1), Ci, carboxylate conc (μmol g–1 root dw), malonic (μmol plant–1), Pn_area (μmol m–2s–1), Pn_mass (μmol g–1 s–1), specific rhizosheath weight (g g–1root DW) were identified using FarmCPU. No significant MTAs were identified for these traits using GLM and MLM models.
FIGURE 3
Discussion
Understanding the genetics of traits associated with enhanced PUE, including P-acquisition efficiency and P-use efficiency, is essential for its manipulation. In our earlier studies on chickpea, we gained insight into root system architecture, shoot/root traits associated with P-acquisition efficiency and P-use efficiency, and associated proxy traits. Many SNPs associated with root-hair length in P-sufficient and P-deficient conditions have been reported recently (); here we report SNPs associated with PUE and P-acquisition efficiency related root traits as well as proxy traits, which can be deployed for breeding.
As population structure is important for avoiding spurious associations, we identified three subpopulations using ADMIXTURE in the reference set (233 genotypes) and a smaller subset (91 genotypes). Similarly, earlier studies reported three subpopulations using different marker systems () and genome-wide SNP markers (). In the smaller subset of 91 genotypes, we also identified three subpopulations (Supplementary Figures 1A–D). Four subpopulations were recently reported in a diverse set of 186 genotypes (including 20 Iranian landraces and 166 Kabuli advanced breeding lines from ICRISAT and ICARDA) using DArTseq markers (). Two subpopulations were reported in a set of 92 (77 landraces and five elite cultivars) chickpea germplasm lines that represent arid, semi-arid, and tropical climates using 8,113 genotyping-by-sequencing based SNPs ().
Three statistical models (GLM, MLM, and FarmCPU) were used to identify genome-wide association signals for root architectural traits, P-acquisition efficiency and PUE-related traits. Similarly, a recent study used different statistical models to compare selected traits with different heritabilities in soybean and maize—the FarmCPU model provided a closer number of QTL than those in the literature and known genomic regions (). In this study, identification of significant MTAs and determining their usefulness for chickpea improvement was a major focus, rather than comparative utilities of different statistical models. In the present study, significant MTAs for 17 of the 30 root traits analyzed was reported. Earlier, in chickpea, MTAs for drought-tolerance-related root traits were based on screening a reference set in polyvinyl chloride pipes or semi-automated root screening facility at ICRISAT (; ). In this study, MTAs for root traits studied under semi-hydroponic conditions and the role of root and shoot traits in P-acquisition efficiency and PUE are reported.
The MTA (Ca1_4716136) identified for RGR, is present in a gene (Ca_00555) that encodes a receptor-like cytosolic serine/threonine-protein kinase RBK2 involved in protein phosphorylation in Arabidopsis thaliana1. Rac-like GTP-binding protein (ARAC5), with 87.8% homology, is localized on the plasma membrane of Arabidopsis thaliana root tips. The Ca_01156 gene encodes Acyl-CoA-binding domain-containing protein 4; in Arabidopsis thaliana; overexpression of Acyl-CoA -binding protein 3 (ACBP3) results in leaf senescence (). In chickpea, the pectinesterase gene was downregulated in the roots of salt-tolerant genotypes (). Pectinestarases also play a key role in root-hair initiation and elongation (). Similarly, another MTA on Ca1 (Ca1_8712480) explained about 20% of the phenotypic variation detected in the GLM and MLM models. The MTA was present in the Ca_02941 gene that encodes 4-coumarate–CoA ligase 1. In tobacco, root length increased by 64% compared with the wildtype, on overexpression of Fm4CL-like 1 [4-coumarate:CoA ligase 4 (4CL-like 1) from Fraxinus mandshurica] under mannitol-simulated drought stress ().
Two MTAs (Ca4_38518152 and Ca4_8269508) for shoot P content were present in two genes, Ca_13110 and Ca_08315, which encode NAD kinase 2, chloroplastic-like isoform X2 and piezo-type mechanosensitive ion channel homolog isoform X2, respectively. NAD kinase 2 is involved in phosphorylation. In general, NADK genes show tissue specificity in expression. In Arabidopsis, NADK2 is expressed in leaves, while TaNADK2 is highly expressed in wheat pistils, caryopses, and endosperm during the reproductive stage (). All MTAs for shoot P content were robust and explained 27.7–41.4% of the phenotypic variation. QTL for shoot P content and PUE were located on chromosomes 3 and 4, respectively (). Of the five MTAs on Ca6, two were in the same gene (Ca_10411), encoding dynamin-related protein 1C. The MTA identified in the FarmCPU model, Ca1_16163105 on Ca1, was present in the Ca_06938 gene that encodes organic cation/carnitine transporter 4-like (OCT), which is involved in homeostasis in animals and has been well-studied in Arabidopsis thaliana—disruption of AtOCT1 affects root development (). Four major pathways in Lupinus albus that contribute to PUE are carbon fixation, cluster-root formation, soil P mobilization, and cellular P reuse (). A recent effort to understand the genetic basis of photosynthesis and PUE affecting yield reported that three major QTL (q14-2, q15-2, and q19-2) explained 6.6–58.9% of the phenotypic variation (). Furthermore, the gene that encodes purple acid phosphatase within the q19-2 region (Glyma.19G193900) is a potential candidate for regulating both soybean PUE and photosynthetic capacity. A total of 159 QTNs within 31 genomic regions and genes associated with photosynthesis-related traits under P stress conditions were genotyped in 2,019 soybean accessions using 292,035 high-quality SNPs and phenotyped under adequate- and low-P conditions for 2 years (). The MTA (Ca1_12310101) for physiological PUE is present in a gene that encodes fasciclin-like arabinogalactan protein 12 (Ca_02521); in Arabidopsis, this gene is involved in cell wall biogenesis (Figure 2C). Further, this SNP locus was associated with total shoot P content and shoot DW in different models (Supplementary Table 1 and Figure 2). Such shared associations were also reported using different models in mungbean PUE (). A recent study reported that low expression of selected fasciclin-like arabinogalactan protein genes led to kernel abortion in maize (Zea mays) and Arabidopsis thaliana seeds (). In cereals such as rice, SNP loci on chromosomes 1, 4, 11, and 12 are associated with PUE ().
reported that root foraging and root physiology, such as the exudation of carboxylates into the rhizosphere, are important strategies for plant P acquisition efficiency. A positive correlation was also identified between mature leaf Mn concentration and rhizosheath carboxylate amount relative to root DW, and hence the carboxylate-releasing P-mobilizing strategy was proxied by foliar Mn concentration in a large set of chickpea germplasm under low P supply (). The MTA for Mn concentration in mature leaves identified on Ca4 was in gene Ca_14893 that encodes zinc finger BED domain-containing protein RICESLEEPER 2-like protein. Carboxylate exudation is an important physiological root trait that enables plants to mine soil P (, ; ). Carboxylate concentrations in the rhizosheath are positively correlated with shoot P content (). A SNP locus associated with carboxylate amount was identified in the rhizosheath (μmol g–1 root DW) associated with a Ca_11019 gene that encodes for ABC1 family protein, and is an integral part of the membrane. The FarmCPU model identified two significant MTAs (Ca4_15651907 and Ca8_230753) in RhizoPH.
Conclusion
In summary, the SNP loci associated with more than one trait were identified. For instance, Ca1_12310101 on Ca1 is associated with three traits (i.e., physiological PUE, shoot DW, and shoot P content), Ca2_31290805 on Ca2 is associated with P utilization and total root length, Ca4_37796452 on Ca4 is associated with RMR and shoot P content, Ca3_6798755 on Ca3 is associated with RGR and TRLz1, and Ca7_5414752 on Ca7 is associated with leaf intracellular CO2 concentration and WUE. The MTA for Mn in mature leaves identified on Ca4 was in gene Ca_14893 that encodes zinc finger BED domain-containing protein RICESLEEPER 2-like protein. The MTAs reported in this study can be used in chickpea breeding programs to enhance PUE.
Statements
Data availability statement
Publicly available datasets were analyzed in this study. This data can be found here: CNSA (https://db.cngb.org/cnsa/) of CNGBdb with accession code CNP0000370.
Author contributions
KHMS and RKV conceived the study. MT and DK performed GWAS. YC and JP generated the phenotyping data. PB, MRo, AC, MRy, and HL contributed to the resources and the writing. MT, RKV, and KHMS prepared the manuscripts. All authors read and approved the final manuscript.
Funding
RKV greatly acknowledges financial support from the Department of Agriculture and Cooperation and Farmers Welfare, Government of India, UWA Research Support from the University of Western Australia, Australia and the Bill and Melinda Gates Foundation, United States.
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fpls.2021.636973/full#supplementary-material
Supplementary Figure 1Example of spurious associations detected for taproot length (TRL, mm) in the GLM model, with no associations detected in the MLM and FarmCPU models.
Supplementary Figure 2Population structure. (A,B) Three sub-populations among 233 genotypes phenotyped for root and phosphorus-related traits, and (C,D) three sub-populations among 91 genotypes phenotyped for shoot/root morphological and physiological traits. Each colored vertical line represents proportions of ancestral populations (K) for each individual. Optimum K value determined using ADMIXTURE’s cross-validation procedure.
Supplementary Table 1Details of significant MTAs with different traits and functional annotation of SNP loci.
Footnotes
References
1
AssefaT.ZhangJ.Chowda-ReddyR. V.LauterA. N. M.SinghA.O’RoukeJ. A.et al (2020). Deconstructing the genetic architecture of iron deficiency chlorosis in soybean using genome-wide approaches.BMC Plant Biol.20:42. 10.1186/s12870-020-2237-5
2
BharadwajC.TripathiS.SorenK. R.ThudiM.SinghR. K.SheoranS.et al (2020). Introgression of the “QTL-hotspot” region enhances drought tolerance and grain yield in multiple genetic backgrounds of elite chickpea cultivars in India.Plant Genome14:e20076. 10.1002/tpg2.20076
3
BrowningS. R.BrowningB. L. (2016). Genotype imputation with millions of reference samples.Am. J. Hum. Genet.98116–126. 10.1016/j.ajhg.2015.11.020
4
CagnolaJ. I.Dumont de ChassartG. J.IbarraS. E.ChimentiC.RicardiM. M.DelzerB.et al (2018). Reduced expression of selected fasciclin-like arabinogalactan protein genes associates with the abortion of kernels in field crops of Zea mays (maize) and of Arabidopsis seeds.Plant Cell Environ.41661–674. 10.1111/pce.13136
5
ChaudharyR. S.SomasundaramJ.MandalK. G.HatiK. M. (2018). Enhancing water and phosphorus use efficiency through moisture conservation practices and optimum phosphorus application in rainfed maize–chickpea system in vertisols of Central India.Agric. Res.7176–186. 10.1007/s40003-018-0316-8
6
ChenX.WangH.LiX.MaK.ZhanY.ZengF. (2019). Molecular cloning and functional analysis of 4-Coumarate:CoA ligase 4 (4CL-like 1) from Fraxinus mandshurica and its role in abiotic stress tolerance and cell wall synthesis.BMC Plant Biol.19:231. 10.1186/s12870-019-1812-0
7
ChenY.GhanemM. E.SiddiqueK. H. M. (2017). Characterising root trait variability in chickpea (Cicer arietinum L.) germplasm.J. Exp. Bot.681987–1999.
8
CongW.-F.SuriyagodaL. D. B.LambersH. (2020). Tightening the phosphorus cycle through phosphorus-efficient crop genotypes.Trends Plant Sci.25967–975. 10.1016/j.tplants.2020.04.013
9
CosgroveD. J. (2016). Plant cell wall extensibility: connecting plant cell growth with cell wall structure, mechanics, and the action of wall-modifying enzymes.J. Exp. Bot.67463–476. 10.1093/jxb/erv511
10
DiapariM.SindhuA.BettK.DeokarA.WarkentinT. D.Tar’anB. (2014). Genetic diversity and association mapping of iron and zinc concentrations in chickpea (Cicer arietinum L.).Genome57459–468. 10.1139/gen-2014-0108
11
FarahaniS.MalekiM.MehrabiR.KanouniH.SchebenA.BatleyJ.et al (2019). Whole genome diversity, population structure, and linkage disequilibrium analysis of chickpea (Cicer arietinum L.) genotypes using genome-wide DArTseq-based SNP markers.Genes10:676. 10.3390/genes10090676
12
GaliK. K.SackvilleA.TafesseE. G.LachagariV. B. R.McPheeK.HyblM.et al (2019). Genome-wide association mapping for agronomic and seed quality traits of field pea (Pisum sativum L.).Front. Plant Sci.10:1538. 10.3389/fpls.2019.01538
13
GiovannettiM.GöschlC.DietzenC.AndersenS. U.KoprivaS.BuschW. (2019). Identification of novel genes involved in phosphate accumulation in Lotus japonicus through genome wide association mapping of root system architecture and anion content.PLoS Genet15:e1008126. 10.1371/journal.pgen.1008126
14
HammondJ. P.BroadleyM. R.WhiteP. J.KingG. J.BowenH. C.HaydenR.et al (2009). Shoot yield drives phosphorus use efficiency in Brassica oleracea and correlates with root architecture traits.J. Exp. Bot.601953–1968. 10.1093/jxb/erp083
15
HeJ.JinY.DuY.-L.WangT.TurnerN. C.YangR.-P.et al (2017). Genotypic variation in yield, yield components, root morphology and architecture, in soybean in relation to water and phosphorus supply.Front. Plant Sci.8:1499. 10.3389/fpls.2017.01499
16
HiremathP. J.KumarA.PenmetsaR. V.FarmerA.SchlueterJ. A.ChamarthiS. K.et al (2012). Large-scale development of cost-effective SNP marker assays for diversity assessment and genetic mapping in chickpea and comparative mapping in legumes.Plant Biotechnol. J.10716–732. 10.1111/j.1467-7652.2012.00710.x
17
JadhavA. A.RayateS. J.MhaseL. B.ThudiM.ChitikineniA.HarerP. N.et al (2015). Marker-trait association study for protein content in chickpea (Cicer arietinum L.).J. Genet.94279–286. 10.1007/s12041-015-0529-6
18
JaganathanD.ThudiM.KaleS.AzamS.RoorkiwalR.GaurP. M.et al (2015). Genotyping-by-sequencing based intra-specific genetic map refines a “QTL-hotspot” region for drought tolerance in chickpea.Mol. Genet. Genom.290559–571. 10.1007/s00438-014-0932-3
19
KaashyapM.FordR.KudapaH.JainM.EdwardsD.VarshneyR. K.et al (2018). Differential regulation of genes involved in root morphogenesis and cell wall modification is associated with salinity tolerance in chickpea.Sci. Rep.8:4855.
20
KaleS. M.JaganathanD.RuperaoP.ChenC.PunnaR.KudapaH.et al (2015). Prioritization of candidate genes in “QTL-hotspot” region for drought tolerance in chickpea (Cicer arietinum L.).Sci. Rep.5:15296.
21
KalerA. S.GillmanJ. D.BeissingerT.PurcellL. C. (2020). Comparing different statistical models and multiple testing corrections for association mapping in soybean and maize.Front. Plant Sci.10:1794. 10.3389/fpls.2019.01794
22
KatuuramuD. N.HartJ. P.PorchT. G.GrusakM. A.GlahnR. P.CichyK. A. (2018). Genome−wide association analysis of nutritional composition−related traits and iron bioavailability in cooked dry beans (Phaseolus vulgaris L.).Mol. Breed.38:44.
23
KhazaeiH.PodderR.CaronC. T.KunduS. S.DiapariM.VandenbergA.et al (2017). Marker–trait association analysis of iron and zinc concentration in lentil (Lens culinaris Medik.) seeds.Plant Genome10:lantgenome2017.02.0007. 10.3835/plantgenome2017.02.0007
24
KohliP. S.Kumar VermaP.VermaR.ParidaS. K.ThakurJ. K.GiriJ. (2020). Genome-wide association study for phosphate deficiency responsive root hair elongation in chickpea.Funct. Integr. Genomics20775–786. 10.1007/s10142-020-00749-6
25
LambersH.HayesP. E.LalibertéE.OliveiraR. S.TurnerB. L. (2015). Leaf manganese accumulation and phosphorus-acquisition efficiency.Trends Plant Sci.2083–90. 10.1016/j.tplants.2014.10.007
26
LambersH.RavenJ. A.ShaverG. R.SmithS. E. (2008). Plant nutrient-acquisition strategies change with soil age.Trends Ecol. Evol.2395–103. 10.1016/j.tree.2007.10.008
27
Lelandais-BrièreC.JovanovicM.TorresG. A.PerrinY.LemoineR.Corre-MenguyF.et al (2007). Disruption of AtOCT1, an organic cation transporter gene, affects root development and carnitine-related responses in Arabidopsis.Plant J.51154–164. 10.1111/j.1365-313x.2007.03131.x
28
LiB.-B.WangX.TaiL.MaT.-T.ShalmaniA.LiuW.-T.et al (2018). NAD kinases: metabolic targets controlling redox co-enzymes and reducing power partitioning in plant stress and development.Front. Plant Sci.9:379. 10.3389/fpls.2018.00379
29
LiH.YangY.ZhangH.ChuS.ZhangX.YinD.et al (2016). A genetic relationship between phosphorus efficiency and photosynthetic traits in soybean as revealed by QTL analysis using a high-density genetic map.Front. Plant Sci.7:924. 10.3389/fpls.2016.00924
30
LipkaA. E.TianF.WangQ.PeifferJ.LiM.BradburyP. J.et al (2012). GAPIT: genome association and prediction integrated tool.Bioinformatics282397–2399. 10.1093/bioinformatics/bts444
31
LüH.YangY.LiH.LiuQ.ZhangJ.YinJ.et al (2018). Genome-wide association studies of photosynthetic traits related to phosphorus efficiency in soybean.Front. Plant Sci.9:1226. 10.3389/fpls.2018.01226
32
MacDonaldG. K.BennettE. M.PotterP. A.RamankuttyN. (2011). Agronomic phosphorus imbalances across the world’s croplands.Proc. Natl. Acad. Sci.1083086–3091. 10.1073/pnas.1010808108
33
MamidiS.LeeR. K.GoosJ. R.McCleanP. E. (2014). Genome−wide association studies identify seven major regions responsible for iron deficiency chlorosis in soybean (Glycine max).PLoS One9:e107469. 10.1371/journal.pone.0107469
34
MannurD. M.BabbarA.ThudiM.SabbavarapuM. M.RoorkiwalM.SharanabasappaY.et al (2019). Super Annigeri 1 and improved JG 74: two Fusarium wilt-resistant introgression lines developed using marker-assisted backcrossing approach in chickpea (Cicer arietinum L.).Mol. Breed.39:2.
35
MogollónJ. M.BeusenA. H. W.van GrinsvenH. J. M.WesthoekH.BouwmanA. F. (2018). Future agricultural phosphorus demand according to the shared socioeconomic pathways.Global Env. Change50149–163. 10.1016/j.gloenvcha.2018.03.007
36
MotsaraM. R. (2002). Fertility status of Indian soils.Fertilizer News4715–22.
37
NayakS. N.ZhuH.VargheseN.DattaS.ChoiH.-K.HorresR.et al (2010). Integration of novel SSR and gene-based SNP marker loci in the chickpea genetic map and establishment of new anchor points with Medicago truncatula genome.Theor. Appl. Genet.1201415–1441. 10.1007/s00122-010-1265-1
38
OzkuruE.AtesD.NemliS.ErdogmusS.KaracaN.YilmazH.et al (2018). Association mapping of loci linked to copper, phosphorus, and potassium concentrations in the seeds of C. arietinum and C. reticulatum.Genomics1111873–1881. 10.1016/j.ygeno.2018.12.010
39
PangJ.ZhaoH.BansalR.BohuonE.LambersH.RyanM. H.et al (2018c). Leaf transpiration plays a role in phosphorus acquisition among a large set of chickpea genotypes.Plant Cell Environ.412069–2079.
40
PangJ.BansalR.ZhaoH.BohuonE.LambersH.RyanM. H.et al (2018a). The carboxylate-releasing phosphorus-mobilizing strategy can be proxied by foliar manganese concentration in a large set of chickpea germplasm under low phosphorus supply.New Phytologist219518–529. 10.1111/nph.15200
41
PangJ.RyanM. H.LambersH.SidddiqueK. H. M. (2018b). Phosphorus acquisition and utilisation in crop legumes under global change.Curr. Opin. Plant Biol.45248–254. 10.1016/j.pbi.2018.05.012
42
PurushothamanR.KrishnamurthyL.UpadhyayaH. D.VadezV.VarshneyR. K. (2017). Root traits confer grain yield advantages under terminal drought in chickpea (Cicer arietinum L.).Field Crops Res.201146–161. 10.1016/j.fcr.2016.11.004
43
ReddyV. R. P.DasS.DikshitH. K.MishraG. P.AskiM.MeenaS. K.et al (2020). Genome-wide association analysis for phosphorus use efficiency traits in mungbean (Vigna radiata L. Wilczek) using genotyping by sequencing approach.Front. Plant Sci.11:537766. 10.3389/fpls.2020.537766
44
RichardsonA. E.SimpsonR. J. (2011). Soil microorganisms mediating phosphorus availability.Plant Physiol.156989–996. 10.1104/pp.111.175448
45
RoorkiwalM.BharadwajC.BarmukhR.DixitG. P.ThudiM.GaurP. M.et al (2020). Integrating genomics for chickpea improvement: achievements and opportunities.Theor. Appl. Genet.1331703–1720. 10.1007/s00122-020-03584-2
46
SabS.LokeshaR.MannurD. M.Somasekhar, JadhavK.MallikarjunaB. P.et al (2020). Genome-wide SNP discovery and mapping QTLs for seed iron and zinc concentrations in chickpea (Cicer arietinum L.).Front. Nutr.7:559120. 10.3389/fnut.2020.559120
47
SabbavarapuM. M.SharmaM.ChamarthiS. K.SwapnaN.RathoreA.ThudiM.et al (2013). Molecular mapping of QTLs for resistance to Fusarium wilt (race 1) and Ascochyta blight in chickpea (Cicer arietinum L.).Euphytica193121–133. 10.1007/s10681-013-0959-2
48
SaniS. G. A. S.ChangP. L.ZubairA.Carrasquilla-GarciaN.CordeiroM.PenmetsaR. V.et al (2017). Genetic diversity, population structure, and genetic correlation with climatic variation in chickpea (Cicer arietinum) landraces from Pakistan.Plant Genome11:170067. 10.3835/plantgenome2017.08.0067
49
SyersJ. K.JohnstonA. A.CurtinD. (2008). Efficiency of Soil and Fertilizer Phosphorus Use. Reconciling Changing Concepts of Soil Phosphorus Behaviour with Agronomic Information.Rome: Food and Agriculture Organization of the United Nations (FAO), 110.
50
ThudiM.BohraA.NayakS. N.VargheseN.ShahT. M.PenmetsaR. V.et al (2011). Novel SSR markers from BAC-end sequences, DArT arrays and a comprehensive genetic map with 1,291 marker loci for chickpea (Cicer arietinum L.).PLoS One6:e27275. 10.1371/journal.pone.0027275
51
ThudiM.KhanA. W.KumarV.GaurP. M.KattaA. V. S. K.GargV.et al (2016b). Whole genome re-sequencing reveals genome wide variations among parental lines of mapping populations in chickpea (Cicer arietinum).BMC Plant Biol.16:10. 10.1186/s12870-015-0690-3
52
ThudiM.ChitikineniA.LiuX.HeW.RoorkiwalM.YangW.et al (2016a). Recent breeding programs enhanced genetic diversity in both desi and kabuli varieties of chickpea (Cicer arietinum L.).Sci. Rep.6:38636.
53
ThudiM.UpadhyayaH. D.RathoreA.GaurP. M.KrishnamurthyL.RoorkiwalM.et al (2014). Genetic dissection of drought and heat tolerance in chickpea through genome-wide and candidate gene-based association mapping approaches.PLoS One9:e96758. 10.1371/journal.pone.0096758
54
UpadhyayaH. D.BajajD.NarnoliyaL.DasS.KumarV.GowdaC. L. L.et al (2016). Genome-wide scans for delineation of candidate genes regulating seed-protein content in chickpea.Front. Plant Sci.7:302. 10.3389/fpls.2016.00302
55
UpadhyayaH. D.DwivediS. L.BaumM.VarshneyR. K.UdupaS. M.GowdaC. L. L.et al (2008). Genetic structure, diversity, and allelic richness in composite collection and reference set in chickpea (Cicer arietinumL.).BMC Plant Biol.8:106. 10.1186/1471-2229-8-106
56
VarshneyR. K.SongC.SaxenaR. K.AzamS.YuS.SharpeA. G.et al (2013b). Draft genome sequence of chickpea (Cicer arietinum) provides a resource for trait improvement.Nat. Biotechnol.31240–246.
57
VarshneyR. K.GaurP. M.ChamarthiS. K.KrishnamurthyL.TripathiS.KashiwagiJ.et al (2013a). Fast-track introgression of “QTL-hotspot” for root traits and other drought tolerance traits in JG 11, an elite and leading variety of chickpea.Plant Genome6:lantgenome2013.07.0022. 10.3835/plantgenome2013.07.0022
58
VarshneyR. K.ThudiM.NayakS. N.GaurP. M.KashiwagiJ.KrishnamurthyL.et al (2014). Genetic dissection of drought tolerance in chickpea (Cicer arietinum L.).Theor. Appl. Genet.127445–462.
59
VarshneyR. K.ThudiM.RoorkiwalM.HeW.UpadhyayaH. D.YangW.et al (2019). Resequencing of 429 chickpea accessions from 45 countries provides insights into genome diversity, domestication and agronomic traits.Nat. Genet.51857–864. 10.1038/s41588-019-0401-3
60
WissuwaM.KondoK.FukudaT.MoriA.RoseM. T.Pariasca-TanakaJ.et al (2015). Unmasking novel loci for internal phosphorus utilization efficiency in rice germplasm through genome-wide association analysis.PLoS One10:e0124215. 10.1371/journal.pone.0124215
61
XiaoS.GaoW.ChenQ. F.ChanS. W.ZhengS. X.MaJ.et al (2010). Overexpression of Arabidopsis acyl-CoA binding protein ACBP3 promotes starvation-induced and age-dependent leaf senescence.Plant Cell221463–1482. 10.1105/tpc.110.075333
62
XuW.ZhangQ.YuanW.XuF.AslamM. M.MiaoR.et al (2020). The genome evolution and low-phosphorus adaptation in white lupin.Nat. Commun.11:1069.
63
ZhouH.AlexanderD.LangeK. (2011). A quasi-Newton acceleration for high-dimensional optimization algorithms.Stat. Comput.21261–273. 10.1007/s11222-009-9166-3
Summary
Keywords
chickpea, genome-wide association study, phosphorus-acquisition efficiency, phosphorus-use efficiency, root traits, genetic mapping
Citation
Thudi M, Chen Y, Pang J, Kalavikatte D, Bajaj P, Roorkiwal M, Chitikineni A, Ryan MH, Lambers H, Siddique KHM and Varshney RK (2021) Novel Genes and Genetic Loci Associated With Root Morphological Traits, Phosphorus-Acquisition Efficiency and Phosphorus-Use Efficiency in Chickpea. Front. Plant Sci. 12:636973. doi: 10.3389/fpls.2021.636973
Received
02 December 2020
Accepted
04 May 2021
Published
28 May 2021
Volume
12 - 2021
Edited by
Aditya Pratap, Indian Institute of Pulses Research (ICAR), India
Reviewed by
Kun Lu, Southwest University, China; Mubshar Hussain, Bahauddin Zakariya University, Pakistan
Updates
Copyright
© 2021 Thudi, Chen, Pang, Kalavikatte, Bajaj, Roorkiwal, Chitikineni, Ryan, Lambers, Siddique and Varshney.
This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Kadambot H. M. Siddique, kadambot.siddique@uwa.edu.auRajeev K. Varshney, r.k.varshney@cgiar.org; rajeev.varshney@murdoch.edu.au
This article was submitted to Plant Breeding, a section of the journal Frontiers in Plant Science
Disclaimer
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.