A meta-QTL analysis highlights genomic hotspots associated with phosphorus use efficiency in rice (Oryza sativa L.)

Phosphorus use efficiency (PUE) is a complex trait, governed by many minor quantitative trait loci (QTLs) with small effects. Advances in molecular marker technology have led to the identification of QTLs underlying PUE. However, their practical use in breeding programs remains challenging due to the unstable effects in different genetic backgrounds and environments, interaction with soil status, and linkage drag. Here, we compiled PUE QTL information from 16 independent studies. A total of 192 QTLs were subjected to meta-QTL (MQTL) analysis and were projected into a high-density SNP consensus map. A total of 60 MQTLs, with significantly reduced number of initial QTLs and confidence intervals (CI), were identified across the rice genome. Candidate gene (CG) mining was carried out for the 38 MQTLs supported by multiple QTLs from at least two independent studies. Genes related to amino and organic acid transport and auxin response were found to be abundant in the MQTLs linked to PUE. CGs were cross validated using a root transcriptome database (RiceXPro) and haplotype analysis. This led to the identification of the eight CGs (OsARF8, OsSPX-MFS3, OsRING141, OsMIOX, HsfC2b, OsFER2, OsWRKY64, and OsYUCCA11) modulating PUE. Potential donors for superior PUE CG haplotypes were identified through haplotype analysis. The distribution of superior haplotypes varied among subspecies being mostly found in indica but were largely scarce in japonica. Our study offers an insight on the complex genetic networks that modulate PUE in rice. The MQTLs, CGs, and superior CG haplotypes identified in our study are useful in the combination of beneficial alleles for PUE in rice.


Introduction
Phosphorus (P) is one of the most important macronutrients in plants and is required in large quantities.Inorganic phosphate (Pi) is a crucial component of phospholipids and plays a significant role in deoxyribonucleic acid (DNA) and ribonucleic acid (RNA) synthesis and other nucleotide-containing molecules.Moreover, P present in adenosine triphosphate (ATP), adenosine diphosphate (ADP), and nicotinamide adenine dinucleotide phosphate (NADPH), provides plant cells with energy required for metabolic and catabolic cellular processes (Heuer et al., 2017).In rice (Oryza sativa), P deficiency substantially decreases overall productivity by reducing plant height, tiller number, panicle length (Irfan et al., 2020) and suppresses the root growth (Liu, 2021).Delayed maturity, high sterility, and poor grain quality (Ismail et al., 2007) have also been reported.
Approximately 5.7 billion hectares of land lack plant-available P due to its immobility in soil (Dobermann and Fairhurst, 2000).P tends to get fixed in soils with extreme levels of pH due to the complexation of P by aluminum (Al) or iron (Fe) in acidic conditions and by calcium (Ca) in alkaline soils (Haefele et al., 2014).Unlike N, which has an unlimited source due to its abundance in the atmosphere thanks to the Haber-Bosch process (Lynch, 2011), P fertilizer sources are finite, mainly consisting of phosphate rocks that are estimated to last for only 300~400 years (Shimizu et al., 2004;Van Kauwenbergh, 2010).
Plants have evolved to cope with P-starvation by undergoing physiological changes in their root morphology, such as the promotion of lateral root and root hair growth, and the inhibition of primary root development.Additionally, plants foster symbiosis between their roots and arbuscular mycorrhizal fungi to scavenge Pi from the soil (Peŕet et al., 2011).The uptake of P in plants is determined by three primary factors: (1) root morphology; (2) P uptake efficiency; and (3) internal P-use efficiency.Among these mechanisms, genetic variations in root morphology are the primary causal factor of P uptake, whereas P uptake efficiency and internal P-use contribute less (Ismail et al., 2007).Consequently, enhancing root morphology through breeding efforts may lead to the development of rice varieties with high phosphorus use efficiency (PUE).
Development of P-efficient rice varieties can be achieved through improved uptake of phosphate from soil (P-acquisition efficiency) (Mori et al., 2016) and improved biomass and/or yield per unit P taken up [internal P-utilization efficiency/P-use efficiency] (PUE) (Rose and Wissuwa, 2012).Developing Pefficient rice genotypes has gained significant attention to breeders.The major PUE quantitative trait locus (QTL), Pup1, has been extensively utilized due to its large additive effect.The QTL was mapped from a backcross population derived from a cross between the P-deficiency intolerant Nipponbare and P-deficiency tolerant Kasalath and is located on the long arm of chromosome 12.The use of Pup1 in molecular marker-assisted backcrossing (MABC) has proven successful (Chin et al., 2011).The Pup1 QTL harbors the Kasalath-derived OsPSTOL1 gene, encoding a protein kinase that enhances early-stage root morphology in rice (Gamuyao et al., 2012).Similarly, a minor QTL on the long arm of chromosome 6 (Shimizu et al., 2004) where three of the known rice Pi-responsive regulatory genes OsERF3, OsTHS1 (Wasaki et al., 2003), and OsPTF1 (Yi et al., 2005) collocate, has been mapped but is yet to be implemented in large-scale breeding programs.In addition, numerous promising genes that are involved in PUE have been identified through overexpression and knockout studies in rice and are now undergoing advanced testing in cereal crops (Heuer et al., 2017).These genes include AVP1 (Yang et al., 2014), PHO1 (Khan et al., 2014), OsPHT1;6 (Ai et al., 2009), and SPX-MFS (Wang et al., 2015).Moreover, the next-generation sequencing (NGS) has facilitated the development of high-density rice linkage maps, which have been utilized to identify QTLs associated with various rice traits, including PUE.Several studies have used SNPbased linkage maps to identify PUE QTLs (Koide et al., 2013;Ogawa et al., 2014;Navea et al., 2017;Fu et al., 2019;Ranaivo et al., 2022;Wang K et al., 2014).
Despite the progress in ascertaining QTLs and genes underlying PUE, their practical utility in breeding programs remain elusive due to their unstable effects across different genetic backgrounds and environments (Arcade et al., 2004;Daware et al., 2017;Navea et al., 2022) as well as their complex relationship with soil status (Heuer et al., 2017) and linkage drag (Kumar et al., 2020).It is therefore necessary to identify genomic regions conferring robust and stable effects across a wide range of genetic background and environments, as well as QTLs with small CI to increase the efficiency and precision of genomics-assisted breeding (Collard and Mackill, 2008).
Meta-QTL analysis is a powerful tool that can help achieve breeding precision by compiling QTLs identified from various mapping populations used in independent studies.It can also provide target genomic regions with considerably small CI (Goffinet and Gerber, 2000;Khahani et al., 2021).Previous studies successfully identified the meta-QTL's (MQTL) underlying yield (Khahani et al., 2021;Aloryi et al., 2022), nitrogen-use efficiency (Sandhu et al., 2021), salinity tolerance (Islam et al., 2019), grain zinc content (Joshi et al., 2023), drought (Selamat and Nadarajah, 2021), and grain traits (Selamat and Nadarajah, 2021) in rice.However, to date, MQTLs associated with PUE remain to be explored in rice.The present study aims to fill this gap by identifying promising MQTLs through an extensive literature search on previous PUE QTLs mapped in various independent studies and identify potential rice donors for pyramiding of the beneficial alleles of the candidate genes (CGs) identified within the MQTLs.

Compilation of PUE QTLs from various studies
An exhaustive literature search on rice QTLs linked to PUE from studies published between 1998 and 2022 was performed using the publicly available QTL databases, such as PubMed (https://pubmed.ncbi.nlm.nih.gov),Google scholar (https:// scholar.google.com),and Gramene QTL database (https:// archive.gramene.org).QTLs and traits associated with PUE were defined as traits contributing to biomass and yield (Rose and Wissuwa, 2012) under low P condition, such as shoot dry weight (SDW); absolute value of root traits (RT); root-shoot ratio (RSR); total dry biomass (BM); seed P content (SPC); relative response to P (RRP); internal P translocation (IPT); and yield component (YLD).QTLs identified under low P input were chosen for the meta-QTL (MQTL) anlaysis.Additionally, only QTL mapping studies carried out under field conditions were selected.Sixteen independent studies were used to perform the MQTL analysis.A total of 192 PUE QTLs mapped from 15 non-overlapping bi-parental populations were utilized in the analysis.Information on the type of parental lines, population types, size of the mapping population, type of molecular markers, number of QTLs, and trait types were recorded.Individual QTLs were assessed for QTL names, type of mapping population, trait class, year of study, linkage group, logarithm of odds (LOD) scores, phenotypic variance explained (PVE), peak position (cM), and CI.These attributes were used to generate the "QTL file".The peak positions were inferred from the midpoint between two flanking markers whenever information on the peak position was lacking.An LOD threshold of 3.0 was set when LOD scores were not supplied.Missing CI values were estimated using the following formulas proposed in a previous study (Darvasi and Seller, 1997;Weller and Soller, 2004): Where N = population size and PVE = phenotypic variance explained.Equation I was used to estimate CI in recombinant inbred lines (RILs) whereas equation II was used for populations derived from F 2 , backcross inbred lines (BILs), and chromosome substitution lines (CSSLs).Missing PVE values were calculated using the following formula (Nagelkerke, 1991):

Construction of consensus map and QTL projection
A high-density consensus map was constructed by integrating SNP markers from the 6K Infinium SNP array (Thomson et al., 2017) and the flanking markers of 192 QTLs.The physical locations of the flanking markers were determined by aligning their sequences to the Nipponbare reference genome (IRGSP v. 1.0) using the Basic Local Alignment Tool (https://www.ncbi.nlm.nih.gov).Then, the closest markers from the Infinium SNP array were used to project QTLs on the consensus map.Markers were arranged according to their physical positions.The consensus map was created by converting the physical position into centimorgan units (cM) using the conversion factor of 1 cM = 250 kb (Raghavan et al., 2017).An individual "map file" for each chromosome was generated, containing information on linkage group, marker name, and the genetic position in centimorgan (cM).The QTL and map files were used as input files to project the consolidated QTLs on the consensus map.

Meta-QTL analysis
After QTL projection, MQTL analysis was performed using Biomercator v4.2.3 software (Arcade et al., 2004;Sosnowski et al., 2012).Two different approaches were applied based on the number of QTLs on each chromosome, namely the Goffinet and Gerber method when the number of QTLs was ≤ 10, and the Veyrieras approach when the number of QTLs was > 10 ( Goffinet and Gerber, 2000;Veyrieras et al., 2007).The algorithms and statistical procedures for both methods were previously described in detail (Sosnowski et al., 2012).The "true" number of MQTLs per chromosome was defined from the model with the lowest Akaike information criterion (AIC) value, which was selected as it contained the least amount of information loss (Akaike, 1987).The PVE value for each MQTL were estimated using that of its "initial QTL" members.The best model, along with its corresponding AIC values, MQTL peak positions, 95% CI, and physical positions, are presented in (Table 1).Additionally, only MQTLs with an average PVE value of ≥ 5% and supported by QTLs identified in at least two independent studies were selected for further analysis.The QTLs identified under P-deficient or P-nonsupplied conditions were categorized into the following traits: 1) shoot dry weight (SDW); 2) absolute values of root traits (RT); 3) root-shoot ratio (RSR); 4) total dry biomass (BM); 5) seed P content (SPC); 6) relative response to P (RRP); 7) internal P translocation (IPT); and 8) yield components (YLD).

Identification of CGs underlying PUE in rice
After identifying MQTLs with an average PVE value of ≥ 5% and support from at least two independent studies using the model with the lowest AIC value, we aimed to mine CGs associated with PUE in rice.To this end, we utilized the physical position of markers flanking the MQTLs as query terms in the Rice Annotation Project Database (IRGSP v1.0 and MSU 7) to batch download functionally annotated gene models within the MQTLs.All genes within the MQTLs were initially considered as CGs, which were further filtered based on gene ontology (GO) terms and/or keywords related to PUE such as P homeostasis, phosphate, inorganic phosphate (Pi) transporter, crown root, root hairs, phosphorus translocation, phosphorus uptake, abiotic stress, and/ or secondary traits associated with PUE, as described in a review article by Heuer et al. (2017).
To identify P-responsive genes, we conducted in silico analysis on the CGs identified in the previous step.We used a microarray dataset (RXP_5002, available at https://ricexpro.dna.affrc.go.jp) containing root gene expression data of 7-day-old seedlings grown under P-deficient and control conditions.Details of the plant growth conditions and treatments can be found on the RiceXPro website (https://ricexpro.dna.affrc.go.jp/RXP_5002/ details-of-methods.html).Briefly, 7-day-old Nipponbare seedlings were exposed to P-non-supplied and control (P-supplied) conditions.Root samples were collected at 6-and 24-hour (h) post-treatment for RNA extraction.The RNA was labeled with Cy3 and subjected to hybridization using the Agilent one-color microarray analysis system.The resulting gene expression profile was presented in terms of raw signal intensity.We carried out gene ontology (GO) enrichment analysis using the Shiny GO 0.76 database (Ge et al., 2020) to determine the most enriched biological pathways in the P-responsive genes with at least 1.5fold change in expression (P-non-supplied/control condition).We used all the unfiltered genes within the MQTLs as the background against the P-responsive genes in the GO enrichment analysis.In the next step, we further filtered the CGs to those that were responsive to P at both 6-and 24-h post-treatment, with a consistent direction of regulation at both time points, and at least a 2-fold change in expression levels at 24-h post-treatment (P-nonsupplied relative to control).

Statistical analysis
Statistical analysis was performed using the two-tailed t-test to determine the significance of differences in the root gene expression levels at different P application at P<0.10, P<0.05, and P<0.01 using Minitab release v.14.

Identification of potential donors for the pyramiding of beneficial CG alleles
Haplotype analysis carried out for the PUE CGs with significant responses to P treatment was performed using the SNP Seek database's built-in tool (Mansueto et al., 2016) in the 3K Rice Genome Project (RGP) database.The database included all rice subpopulations with Nipponbare as the reference genome.PUE CGs with non-synonymous SNPs, at least a 2-fold change in gene expression under P-deficient vs. control conditions, and with at least two haplotypes were utilized in the analysis.The number of haplotypes were determined using the default parameters with the Calinski criterion for determining the optimal number of groups (Caliñski and Harabasz, 1974).Haplotypes with fewer than three genotypes were excluded.Genotypes with more than 20% missing SNPs or heterozygous loci were filtered out from the data set.Haplotypes identical to those of the beneficial allele donors (whenever available in the 3K RGP dataset) in the QTL studies used in the MQTL analysis were regarded as superior haplotypes.The abundance of superior haplotypes across the subpopulations was calculated, and rice accessions with superior haplotypes in indica, japonica, and aus backgrounds from the 3K RGP were identified.These lines were considered as potential donors for pyramiding superior CG haplotypes for improving PUE in rice.

The rice consensus map and QTL projection
The consensus map included a total of 5,694 markers, comprising both SNPs from the 6K Infinium SNP array and the flanking markers of the QTLs used in the MQTL analysis.They were well-distributed across the twelve chromosomes.The SNP density per 500 kb window is presented in (Figure S1).In addition, the positions of known PUE genes, such as OsPSTOL1, OsPTF1, OsERF3, and OsTHS1, were incorporated onto the consensus map.The cumulative length of the consensus map was 1,637 cM.The average distance between the adjacent markers was 0.29 cM.Chromosomes 1 (171.9cM) and 6 (160 cM) were the longest, while chromosomes 9 (108.5 cM) and 10 (92.2 cM) were the shortest.The number of markers varied from 332 to 653 per chromosome.Chromosomes 1 and 2 had the highest counts of 653 and 574 markers, respectively, while chromosomes 9 (n=354) and 10 (n=332) were the least saturated.We projected the PUE QTLs onto the consensus map using both the physical position of the QTLs and the SNP markers (Figure 3).The number of QTLs ranged from 5 to 33.Chromosomes 6 (n=33), 2 (n=32), and 1 (n=25) were the most saturated regions, while chromosomes 9 (n=5), 7 (n=5), and 3 (n=7) had the fewest.

MQTL analysis
MQTL analysis of the 192 initial PUE QTLs resulted in the identification of 60 MQTLs (Table 2; Figure 3).The distribution of MQTLs across chromosomes was uneven.The highest number of MQTLs were detected on chromosomes 1 (n=9) and 5 (n=8), while chromosomes 3 (n=2), 7 (n=2), and 9 (n=2) had the least.Reductions in MQTL CI were observed, with fold reductions ranging from 2.8 to 11.4 with an overall average reduction of 9.7 cM (Figure 4).The CI for MQTLs located on chromosomes 3, 5, 7, 9, and 10 were at least five times smaller than their corresponding initial QTLs, while those on chromosomes 1, 2, 6, and 11 were at least twice as small.Chromosomes 4 and 12 exhibited the least reduction in CI, with fold change values of 1.7 and 1.02, respectively.The MQTLs had PVE values from 2.3 to 20.3.
Subsequently, we selected MQTLs that had an average PVE value of at least five and were supported by QTLs from at least two independent studies (Table 3).This resulted in 38 MQTLs, with the number of QTL members ranging from 2 to 11. MQTL1.1,MQTL8.1, and MQTL8.2 had the most QTLs (n=11), while MQTL3.1,MQTL4.1, and MQTL10.1 had only two initial QTLs.The MQTLs had CI ranging from 0.2 to 9 cM, with an average of 3.5 cM per MQTL.The smallest CI (0.2~2 cM) were observed on
To understand the biological pathways involved, we conducted gene ontology (GO) enrichment analysis using the 103 PUE CGs.Amino acid transmembrane transport, organic transport, and response to auxin were the most enriched biological pathways with five, five, and six gene members, respectively.Transcription, DNA templated, nucleic acid-templated transcription, and regulation of nucleobase-containing compound metabolic process  had the most gene members (20, 21, and 21, respectively) (Figure 5).
To identify genes that are responsive to P treatment, we set a threshold of a 1.5-fold change in expression (P non-supplied condition vs. control), which narrowed down the list of CGs to 103 (Table 4), of which 86 were upregulated and 17 were downregulated (Figure 6).

Haplotype analysis and identification of potential donors
We performed haplotype analysis on the selected PUE CGs and identified potential donors for pyramiding the beneficial PUE CG alleles.We found 25 genes with at least a 2-fold change in expression levels at 24h post-treatment (Table S6).Among these, the eight genes (OsARF8, OsSPX MFS3, OsRING141, OsMIOX, HsfC2b, OsFER2, OsWRKY64, and OsYUCCA11) (Figure 7) had at least two haplotypes in the 3K RGP (Table 5) and were therefore used for further analysis.The number of synonymous SNPs ranged from three to nine (Table S7).OsARF8, OsSPX-MFS3, and HSfC2b had five SNPs, while OsYUCCA11 had the most with nine SNPs.OsRING141, OsMIOX, OsFER2, OsWRKY64 had four SNPs.The number of haplotypes ranged from two to six (Table S8).Four CGs had four haplotypes each (OsARF8, OsMIOX, HSfC2b, OsFER2, and OsWRKY64).Two CGs had six haplotypes each (OsSPX-MF3 and OsYUCCA11).OsRING141 had only two haplotypes.
We inferred the superior haplotypes from the beneficial allele donors in the original QTLs used in MQTL analysis.Kasalath and IR20, which were donors of the beneficial alleles in MQTLs from which the eight PUE CGs were located, had haplotype information in the 3K RGP.Thus, we used them to infer superior haplotypes.The Kasalath-types were regarded as the superior haplotype for OsARF8 (Haplotype 2), OsRING141 (Haplotype 4), OsSPX-MFS3 (Haplotype 4), and HSfC2b (Haplotype 2).IR20-types were selected for OsMIOX (Haplotype 1) OsFER2 (Haplotype 2), OsWRKY64 (Haplotype 2), and OsYUCCA11 (Haplotype 2) (Table S8).The abundance of the superior haplotypes was evaluated in the 3K RGP (Table S9; Figure 8).Superior haplotypes were found most abundant in the indica variety group.The frequency of superior alleles in japonica was considerably lower than that of the other subpopulations; in fact, superior haplotype for OsSPX-MFS3 was completely absent.We identified rice accessions from the 3K RGP that cover the beneficial haplotypes for the eight CGs (Table 6).Four accessions were identified as the potential donors for the indica breeding programs, whereas two and three accessions were selected for the PUE breeding programs for aus and japonica, respectively.The population development using only Kasalath and IR20 can be suggested to pyramid the superior haplotypes of the eight PUE CGs.

Discussion
The development of rice varieties requiring less agricultural inputs is essential in the face of a changing climate and the rising cost of agricultural inputs such as fertilizers.Unlike N fertilizers, P fertilizer sources are estimated to last for only another four to five generations (Shimizu et al., 2004;Van Kauwenbergh, 2010).This is further complicated by the lack of plant-available P, despite fertilizer application due to the inherently low PUE in rice (Rose et al., 2011).It is therefore vital to utilize genetic variation in rice to address these issues.PUE-related QTLs and genes have been previously studied and mapped.However, their practical breeding utility has been hampered by their unstable effect across environments, genetic backgrounds, QTL to QTL or gene to gene interactions, and linkage drag caused by the large introgression size of QTLs that have not undergone fine-mapping.One example is the instability of the Kasalath-derived major PUE QTL Pup1, which confers tolerance to low P application as well as to mild drought.Pup1 improves yield under low P in several indica varieties such as IR64, IR74, Situ Bagendit, Batur, and Dodokan under tropical conditions (Chin et al., 2011).However, Pup1 introgression into Dasanbyeo, a Tongil-type indica, neither improved phosphorus uptake nor yield under low P application in temperate regions (Navea et al., 2022).In addition, a study has revealed the inability of Pup1 to confer tolerance to low P in early shoot growth when introgressed simultaneously with the submergence tolerance major QTL Sub1 into IR64 (Shin et al., 2021).The inconsistencies in QTL  We conducted an MQTL analysis on 192 PUE QTLs identified in 16 independent studies, in 15 non-overlapping bi-parental mapping populations.MQTL analysis offers the benefit of identifying reliable and robust genomic regions with reduced numbers of QTLs and narrower introgression sizes (Goffinet and Gerber, 2000;Kumari et al., 2021;Aloryi et al., 2022;Anilkumar et al., 2022;Joshi et al., 2023).To the best of our knowledge, our study is the first attempt to identify MQTLs underlying PUE in rice.Initially, we identified 60 PUE MQTLs and further selected 38 MQTLs that were supported by at least two independent studies and had PVE values of at least 5%.The majority of QTLs underlying these MQTLs had small PVE values ranging from 3% to 6%, implying that PUE in rice is controlled by many minor loci.This observation is consistent with results from several PUE QTL mapping studies (Navea et al., 2017;Wang et al., 2014).The distribution of the initial QTLs largely varied throughout the whole chromosome (Figure 1).However, it is not largely attributable to the SNP density (Figure S2).This is in contrast to the results from other MQTL studies in rice (Khahani et al., 2021;Selamat and Nadarajah, 2021;Aloryi et al., 2022;Joshi et al., 2023), wherein the uneven distribution of the initial QTLs mostly depended on the number of markers used in constructing the consensus map.Chromosome 6 was the chromosome with the largest number of initial PUE QTLs, as reported in the various previous PUE studies in rice (Ismail et al., 2007;Wissuwa et al., 2015;Heuer et al., 2017).
Linkage drag, due to the large CI of the QTLs utilized in breeding programs, is one of the hindrances in achieving successful in marker-assisted breeding (MAB) and genomic selection (GS) programs (Collard and Mackill, 2008;Kumar et al., 2020).Here, the average CI of the PUE MQTLs was reduced compared to their initial QTL members (Figure 4), except in case of chromosomes 4 and 12. MQTL analysis was not able to narrow down the CI on chromosomes 4 and 12 partly due to the small number of the initial QTL members as well as the already small CI of the QTLs.With the exceptions of the MQTLs detected on chromosomes 4 and 12, the average CI was reduced from 14.04 cM to 2.70 cM, suggesting the power of MQTL analysis in identifying precise genomic segments modulating a trait of The most enriched biological pathways in the PUE CGs.

LOC_Os12g09300 Os12g0194900 OsAAP10B
Amino acid permease, A member of the amino acid transporter (AAT) family, Regulation of tillering and grain yield, Regulation of neutral amino acid transport LOC_Os12g09590 Os12g0197700 -Region of unknown function, putative Zinc finger, XS and XH domain containing protein.
interest.Our results are consistent with the observations in previous MQTL studies on various traits in cereal crops (Chen et al., 2017;Venske et al., 2019;Soriano et al., 2021;Miao et al., 2022).
The successful integration of pleiotropic QTLs into breeding programs has been documented in rice (Ookawa et al., 2010;Yan et al., 2011;Vishnukiran et al., 2020).In this study, we categorized QTLs identified under field conditions, based on traits welldistributed throughout the rice growth stages.These included seedling stage (RT, RSR, RRP), vegetative stage (SDW, BM, IPT), and reproductive stage traits (SPC and YLD).This enabled us to identify PUE MQTLs that may have pleiotropic effects, as observed in MQTL1.5, MQTL1.6, and MQTL1.7, on multiple traits including BM, SCP, RSR, IPT, YLD, and RRP.This can aid in a more efficient improvement of PUE in rice through the accumulation of beneficial PUE MQTL alleles.Among the QTLs used in this study, genomic regions regulating root traits under P deficient conditions were the least abundant across rice chromosomes, which can be attributed to the challenges in establishing reliable and high-throughput root phenotyping techniques under field conditions (Heuer et al., 2017), as well as the lack of genetic diversity (Ismail et al., 2007).It is worth noting that the number of PUE CGs underlying MQTL2.7,MQTL4.3, and MQTL11.4 are relatively large despite their narrow CI.This is unlike the observations in previous MQTL studies in rice, wherein the number of CGs underlying an MQTL had a strong positive correlation with the size of the CI (Daware et al., 2017;Islam et al., 2019;Khahani et al., 2021;Selamat and Nadarajah, 2021;Aloryi et al., 2022;Anilkumar et al., 2022;Joshi et al., 2023).These MQTLs have the potential to be effective genomic targets for use in the MAB after being validated in a wide range of genetic backgrounds and environments by utilizing the aforementioned linked flanking markers (Table 3).
We identified CGs underlying the PUE MQTLs using GO terms directly related to PUE traits, as well as for the secondary PUE traits.We further selected genes that could be validated through the root expression data under control and P deficient/non-supplied conditions.Most of the PUE CGs (83%) were upregulated (Figure 6F), implying that genes underlying the PUE MQTLs confer an active response to nutrient deficiency, rather than conservation of resources.This agrees with the pattern of gene regulation under nutrient deficiency in previous studies in various plants (Loṕez-Bucio et al., 2003;Chen et al., 2022;Wang M et al., 2019).The results of our GO analysis give an insight into the complexity of PUE in rice.CGs underlying PUE MQTLs were heavily enriched with genes involved in amino acid transmembrane transport, organic acid transport, and response to auxin.Amino acid transporters in rice, although not reported to be directly involved in PUE, have been associated with the regulation of flowering time and defense against abiotic stresses and pathogen attack (Guo et al., 2021).In particular, the genes involved in the amino acid transport pathways function in various plant species defending against both biotic and abiotic stresses through the regulation of the salicylic acid pathway (Kan et al., 2017), including drought, salinity, UV radiation, heavy metals, and pathogens (Szabados and Savoure, 2010).In plants, organic acid transporters are upregulated by Al and/or P deficiency (Yu et al., 2016).This is caused by enhanced phosphorylation levels of the plasma membrane H+-ATPase, which creates an electrochemical potential across the cell membrane.This leads to an increase in the activity of the organic acid transporters and the passive release of organic anions from root tips.Ultimately, this process causes organic  acids to exude from the roots and form a stable complex with Al, which allows P to become soluble for plant assimilation (Ligaba et al., 2004).Plant root architecture undergoes adaptive changes, including the inhibition of primary root growth and the increase in the number and length of lateral roots (Peŕet et al., 2011) modulated by the change in sensitivity of auxin receptors such as TIR1, under a P deficient condition (Wu et al., 2023).Upregulation of auxin receptors degrades auxin repressors, releasing auxin response factor, ARF19, which then leads to the activation of genes related to lateral root morphogenesis (Peŕez-Torres et al., 2008).We further narrowed down the PUE CGs to eight genes and investigated their natural variation in a diverse set of rice (3K RGP).These genes included OsARF8, OsSPX-MFS3, OsRING141, OsMIOX, HsfC2b, OsFER2, OsWRKY64, and OsYUCCA11 (Figure 7).OsSPX-MFS3 is a member of the rice SPX-MFS family which mediates Pi transport between the cytosol and vacuole (Yang et al., 2017).OsSPX-MFS3 plays a major role in the transport of Pi from the cytosol to vacuole (Guo et al., 2023) and is downregulated under P deficiency.The regulation pattern of OsSPX-MFS3 is consistent with the observation in the microarray data utilized in this study.It is worth noting that, except for OsSPX-MFS3, the rest of the CGs did not have GO terms directly related to PUE but were nevertheless associated to secondary PUE traits and other abiotic stresses.Auxin is one of the phytohormones that regulate root architecture modifications in response to Pi deficiency (Nacry et al., 2005).In this study, we identified two CGs, namely OsARF8 and OsYUCCA11, implicated in auxin response and which were significantly downregulated under P non-supplied condition.OsARF8 is a member of the ARF family implicated in the crosstalk between auxin signaling and P status (Wang S et al., 2014).
Genes under the zinc-finger family may play an important role in the regulation tolerance to multiple stresses in rice (Deng et al., 2018).The present study identified a P-deficiency-upregulated zincfinger gene, OsRING141, among the genes subtending MQTL6.5.Similarly, previous studies found two C2H2-type zinc finger protein genes, ZOS3-12 and ZOS5-08, to be responsive to P and N deficiencies (Huang et al., 2012).MQTL12.1 harbors two genes modulating iron (Fe) homeostasis (OsFER2) and iron stress tolerance (OsWRKY64).Both genes were significantly upregulated in the P non-supplied condition.Previous studies have shown that the reduction of Fe concentration can lead to the recovery of primary root elongation under low P conditions, suggesting that Fe may play a role in the Pi deficiency-induced reduction of primary root growth (PRG) (Ward et al., 2008).Moreover, research on Arabidopsis has demonstrated that when Fe is deposited in the root tip meristem, it can trigger the accumulation of reactive oxygen species (ROS) and callose, likely through LPR1-dependent redox signaling (Müller et al., 2015).The accumulation of ROS and callose can interfere with cell-to-cell communication that is essential for maintaining the stem cell niche and inhibiting PRG (Müller et al., 2015).
MQTL6.7 harbors two abiotic stress-related genes, namely OsMIOX (drought) and HSfC2b (heat).A previous study reported that the overexpression of OsMIOX in transgenic rice greatly improved growth performance under drought conditions by decreasing oxidative damage (Duan et al., 2012).Severe phosphorus deficiency can lead to changes in the photosynthetic apparatus, such as decreased rates of carbon dioxide assimilation, reduced expression of photosynthesis-related genes, and photoinhibition at the photosystem II level.These changes can potentially cause photo-oxidative stress, which can damage the plant's cells.By preventing oxidative damage, the plant may be able to better tolerate the effects of phosphorus deficiency on its photosynthetic apparatus (Hernańdez and Munne-Bosch, 2015).Therefore, as OsMIOX may prevent oxidative damage, it might be beneficial to phosphorus-starved plants.A heat shock factor (HSF) member, HSfC2b, harbored within MQTL6.7, was upregulated under P deficient condition.The first attempt to elucidate the genome wide expression of these HSFs was conducted by Chauhan et al. (2011).HSF encoding genes showed significant upregulation in abiotic stresses such as cold, drought, and salinity.Our study is the first to show the link between an HSF gene and PUE in rice.In Arabidopsis, the overexpression of HSFs in transgenic plants confer simultaneous tolerances to multiple abiotic stresses such as heat and anoxia (Banti et al., 2010).Nevertheless, further study is needed to elucidate the relationship between HSF and PUE in rice.
The use of NGS-based genotyping methods have aided in the identification of SNPs associated with agronomic traits in rice.However, the applicability of SNPs in breeding programs is constrained by their bi-allelic nature by cross breeding, the presence of uncommon alleles, and the abundance of linkage drag Abundance of superior PUE CG haplotypes in the rice 3K rice genome panel.Values indicate the number of accessions.(Annicchiarico et al., 2017).Considering the gene haplotypes for genome-wide analysis will help overcome SNP marker's limitation.Haplotypes are specific combinations of jointly inherited nucleotides or DNA markers from polymorphic sites in the same chromosomal segment (Stephens et al., 2001;Lu et al., 2010).Haplotype-based breeding holds a promise in accumulating beneficial alleles for a trait of interest in living organisms.This can be achieved through the identification of haplotypes for different genes and utilizing them through techniques such as allele mining, pyramiding, or GS.This approach has resulted in genetic gains in crops (Abbai et al., 2019;Anandan et al., 2022;Du et al., 2022).We identified subspecific-wise potential donors of the superior haplotypes for PUE CGs: OsARF8, OsSPX-MFS3, OsRING141, OsMIOX, HsfC2b, OsFER2, OsWRKY64, and OsYUCCA11 (Table 6).We inferred the superior haplotypes for these genes based on the haplotype of beneficial allele donors (Kasalath and IR20) in the initial QTL studies used for the meta-QTL analysis.The frequencies of superior haplotypes differed among sub-populations, with particular groups exhibiting a higher prevalence (Figure 8; Table S9).Similarly, previous studies have suggested that the distribution of haplotypes is influenced by evolutionary and population genetic factors, such as rates of mutation and recombination, as well as selection pressures (Magwa et al., 2016;Sinha et al., 2020).Superior haplotypes for almost all the CGs were predominant in indica subspecies.The frequency of the superior haplotypes in japonica was considerably smaller, compared to that of other subspecies.As mentioned, OsSPX-MFS3 was completely absent in japonica varieties in 3K RGP.This implies a big opportunity to improve PUE in the japonica varieties using the donors identified in this study.Our results suggest that indica varieties are a rich source of PUE CG superior haplotypes and could be utilized in the genomics-assisted breeding programs for PUE in rice.In the case of the japonica, it is necessary to generate prebreeding lines that encompass the superior PUE CG haplotypes.
We identified potential CGs associated with PUE as well as superior haplotypes that could be used to accumulate beneficial alleles to improve PUE in rice.However, it should be noted that the PUE CG mining pipeline used in the present study was limited by the availability of P-supplied and P-non-supplied root transcriptome data.Initially, we identified 273 PUE genes (Table S1), but we could only analyze 238 CGs (Table S2) that had gene expression data in the RiceXPro database.Therefore, a complete set of transcriptome data would provide offer higher precision in the identification of PUE CGs in rice.Another limitation of our pipeline is the limited scope of the reference genome used in both the annotation of genes and the root expression analysis.We utilized gene annotation from the IRGSP v. 1.0 (Nipponbare) reference genome, a widely used reference genome for rice research, to identify PUE CGs.Here, it was not possible for us to identify genes that were absent in the reference genome and therefore were possibly neglected in our analysis.For instance, the P uptake gene OsPSTOL1 harbored in the P uptake major QTL Pup1 is absent in the reference genome Nipponbare, as well as in the genomes of several other commonly used rice varieties (Chin et al., 2011;Gamuyao et al., 2012).Similarly, Sub1A, a gene conferring submergence tolerance in rice, was not found in the fully sequenced genome of Nipponbare (Xu et al., 2006).The limited  availability of gene information from the reference genome may fail to precisely capture the genetic variability linked to PUE.Consequently, it is important to generate supplementary reference genomes, or a "rice meta-genome," that encompass diverse rice varieties and species to address this limitation.Nevertheless, our approach provides breeders with valuable information regarding the selection of optimal donors for their desired traits from the gene bank.However, the assumed genetic gains resulting from the accumulation of superior haplotypes of PUE CGs require validation in practical breeding programs.

Conclusion
We identified 38 meta-QTLs (MQTLs) for phosphorus use efficiency (PUE) that were supported by multiple QTLs from independent studies, which had a phenotypic variation explained (PVE) value of at least 5%.We subjected the 38 PUE MQTLs to candidate gene (CG) mining.The genomic regions associated with PUE MQTLs were found to be enriched with genes involved in the transmembrane transport of amino acids and organic acids, as well as genes involved in the response to auxin.Some superior haplotypes containing eight CGs for PUE could be considered for the genomics-assisted breeding in rice.

FIGURE 1
FIGURE 1Phenotypic trait classes and chromosome-wise distribution of QTLs utilized in the MQTL analysis for PUE in rice.

FIGURE 4
FIGURE 4    Comparison of average CI between initial the QTLs and MQTLs on the 12 rice chromosomes.Values on top of the lines represent the reduction in the average CI.

FIGURE 6
FIGURE 6 Schematic representation of the distribution patterns of PUE QTLs, MQTLs, and PUE CGs on rice chromosomes.The outermost circle represents the rice karyotype.The second circle outlines the consensus map SNP density (500kb window).The third circle display the initial PUE QTLs used in the MQTL analysis.The fourth inner circle represents the MQTLs supported by QTLs from at least 2 independent studies and with a PVE value of ≥ 5%.The fifth and sixth inner circle represent the P-responsive CGs and the fold change in gene expression (P non-supplied vs. control treatment), respectively.

TABLE 1
QTL studies used in the QTL meta-analysis for PUE in rice.
a CSSL, chromosome substitution lines; c SDW, shoot dry weight; RT, absolute value of root traits; RSR, root-shoot ratio; BM, total dry biomass; SPC, seed P content; RRP, relative response to P; IPT, internal P translocation; YLD, yield component.d LOD,logarithm of odds; PVE, phenotypic variation explained.

TABLE 2
Number of initial QTLs and MQTLs identified on rice chromosomes.
a Initial QTLs: QTLs used for MQTL analysis.

TABLE 3
Summary of rice PUE MQTLs supported by multiple QTLs from at least two independent studies and with at least 5% average PVE values.
effects necessitate the identification of fine-mapped genomic regions related to PUE that are effective across diverse genetic backgrounds and environments.

TABLE 3 Continued
a For the model number in the parentheses of AIC value, see the materials and methods.b SDW, shoot dry weight; RT, absolute value of root traits; RSR, root-shoot ratio; BM, total dry biomass; SPC, seed P content; RRP, relative response to P; IPT, internal P translocation; YLD, yield component.c CI, Confidence interval.

TABLE 4 P
-responsive candidate genes underlying PUE MQTLs supported by multiple QTLs from at least two independent studies and with at least 5% average PVE values.Mini zinc finger protein, A member of the ZF-HD (zinc finger-homeodomain) family, Negative regulation of deep sowing tolerance, Mesocotyl elongation.
S-Domain receptor like kinase-5, Response to drought in a tolerant genotype, Response to submergence.

TABLE 5
PUE-related candidate genes used for the haplotype analysis in the rice 3K genome panel.

TABLE 6
Suggested donors based on a haplotype-based selection for the pyramiding of favorable PUE CG alleles.
'+' represents the presence of superior haplotype of each gene Figure1.Phenotypic trait classes and chromosome-wise distribution of QTLs utilized in the MQTL analysis for PUE in rice.