A total of 132 tropical maize inbred lines were evaluated, composed of 77 field corn and 55 popcorn lines from the Maize Breeding Program of the State University of Maringá (UEM), Maringá, Brazil. These inbred lines represent part of the genetic variability present in the main hybrids and varieties of field corn and popcorn cultivated in Brazil. Seed samples were obtained from the UEM Germplasm Bank and later multiplied in the 2017-2018 summer season to standardize seed germination. Information on the origin of the maize lines is presented in Supplementary Table 1 .

Inbred lines were evaluated under controlled conditions in a greenhouse at the Agronomy Department of the State University of Londrina (UEL) in Londrina, Paraná, Brazil (23°17'34”S, 51°10'24”W, and 550 m altitude) in September 2019. Two independent hydroponic systems were developed for low (2.5 μM) and normal P (250 μM) conditions using Magnavaca nutrient solution (Magnavaca et al., 1987). Each hydroponic system consisted of eight 28-L polyethylene boxes (58 × 40 × 16.5 cm). The hydroponic boxes were interconnected and connected to a 300-L reservoir, that is, the same nutrient solution ran throughout the entire system. The pH of the nutrient solution was maintained at 5.8 ± 0.2 and permanently aerated by a pressurizing water pump (Model Pl400P 40 MCA, Lorenzetti Ltda, Campinas, Brazil).

Seeds were disinfested by immersion in 95% ethanol solution (v/v) for 30 s followed by immersion in 5% H₂O₂ solution (v/v) for 10 min and six washes with sterile deionized water. The seeds were then placed on filter paper (Germtest^TM, 24 × 33 × 0.02 cm) (Cienlab Ltda., Campinas, Brazil), moistened with distilled water, and placed in a growth chamber at a temperature of 27 ± 2°C and 70% relative humidity. Seven days after sowing (at development stage V1), uniform seedlings were selected and transplanted to hydroponic systems after removing the endosperm to eliminate seed reserves. The experiment was carried out using a completely randomized design with three repetitions. After 12 days of growth in nutrient solution (stage V2), the inbred lines root system was digitized in a 600 dpi image using an Epson L3110 scanner (Seiko Epson Corp., California, USA), and the root architecture traits were assessed using the GiaRoots software (Galkovskyi et al., 2012).

Inbred lines were evaluated at the School Farm of UEL (23°17'34”S, 51°10'24”W, and 550 m altitude) during the 2018–2019 summer and 2019 fall/winter seasons. This region has a humid subtropical climate and soil classified as Dystroferric Red Latosol. Two distinct areas were selected to conduct the experiments under low and normal P conditions, carried out in the same areas during the two seasons. P (Mehlich 1) contents in the areas ranged from 4.87–5.11 mg dm³ (low P) to 17.32–19.32 mg dm³ (normal P). The soil physicochemical analyses and other traits of the experimental areas are presented in Supplementary Table 2 .

The experiment was set in a randomized complete block design with three repetitions. Plots consisted of a 4 m long line, 0.45 m spacing between lines, and 0.20 m between plants. Two levels of phosphate fertilization were used in each environment, totaling four independent experiments. Before planting, normal P environments were fertilized with 120 kg P₂O₅ ha^–1, 40 kg K₂O ha^–1, and 40 kg N ha^–1, while low P environments only received 40 kg K₂O ha^–1 and 40 kg N ha^–1. The topdressing nitrogen (N) fertilization was carried out using 180 kg N ha^–1 applied at stage V6, while the other cultural practices were performed according to crop demand. At physiological maturation (stage R6), three uniform and representative plants from each plot were collected for phenotypic evaluations.

The detailed description of the 24 traits evaluated in the greenhouse and field experiments is presented in Table 1 . To determine the P content in plant tissues, samples were oven-dried at 60°C for 72 h and milled in Willey-type knife mill MA340 (Piracicaba, São Paulo, Brazil). Then 0.1 g aliquots were digested in nitroperchloric solution (HNO₃:HClO₄) according to Malavolta et al. (1989). P content was determined by the molybdenum blue spectrophotometric method (Pradhan and Pokhrel, 2013), reading the samples in an Agilent 8453 spectrophotometer (Agilent Technologies, California, USA) at 660 nm.

Data were analyzed using the best linear unbiased predictor (BLUP) and restricted maximum likelihood (REML) methods by the software Selegen–REML/BLUP (Resende, 2016) and R version 3.6.0 (https://www.r-project.org) via the 'lme4' package. The deviance analyses (ANADEV) of the traits obtained from the greenhouse and field trials were carried out using the following mathematical models, respectively:

Where µ is the overall mean, G _i is the random effect of the i-th genotype, P _m is the fixed effect of the m-th level of P, GP _im is the random effect of the genotype × P level interaction, and ϵ _im ~ N(0, σ²) is the random effect of the error associated with each experimental unit.

Where µ is the overall mean, G _i is the random effect of the i-th genotype, B _j/k/m is the random effect of the j-th block within the k-th season and within the m-th P level, S _k is the fixed effect of the k-th season, P _m is the fixed effect of the m-th P level, GS _ik is the random effect of the genotype × season interaction, GP _im is the random effect of the genotype × P level interaction, SP _km is the fixed effect of the yield × P level interaction, GSP _ikm is the random effect of the genotype × yield × P level interaction, and ϵ _ijkn ~ N(0, σ²) is the random effect of the error associated with each experimental unit.

The significance of the ANADEV random effects was verified by the likelihood ratio test (Resende, 2016). The heritability in the broad sense (h²) of the traits evaluated in the greenhouse and field were estimated, respectively, using the following formulas: h 2 = σ ^ g 2 σ ^ g 2   +   σ ^ e 2 and h 2 = σ ^ g 2 σ ^ g 2   +   σ ^ gs 2 s   +   σ ^ e 2 r s , where σ ^ g 2 is the genotypic variance, σ ^ gs 2 is the variance of the genotype × season interaction, σ ^ e 2 is the residual variance, r is the number of repetitions in each season, and s is the number of seasons. The selective accuracy (Ac) was obtained as follows: A c = 1 − PEV σ ^ g 2 , where PEV is the variance of the prediction error of the genotypic values, and σ ^ g 2 is the genotypic variance.

Correlations between traits and principal component analyses (PCA) were performed using BLUP means. The significance of the estimates was verified using the t-test at a 5% probability level (α = 0.05). The correlation estimates were visualized using the correlation network approach. These analyses were performed by the 'qgraph' and 'ggbiplot' packages of the R software version 3.6.0.

Genomic DNA was initially isolated from leaf tissues of inbred lines, as established by Coan et al. (2018). DNA samples were then sent to the Genomic Diversity Institute at Cornell University for genotyping-by-sequencing (GBS) according to the protocol described by Elshire et al. (2011). Sequencing was performed in the Illumina HiSeq 2000 sequencer (Illumina Inc., San Diego, USA). The sequences were deposited in the European Variation Archive (EVA) (GCF_902167145.1).

Sequencing data were analyzed using the Tassel 5.0 GBS v2 software (Glaubitz et al., 2014). The sequences obtained were aligned to the reference genome of Zea mays L. version AGPV3 (B73 RefGen v3) obtained from the MaizeGDB database (https://www.maizegdb.org). SNP markers quality control was performed using the VCFtools software version 0.1.15 (Danecek et al., 2011), removing the SNPs through the following criteria: i) non-biallelic, ii) minor allele frequency (MAF) less than 5%, iii) inbreeding coefficient less than 90%, and iv) missing data greater than 90%. Heterozygous SNPs were treated as missing data since the evaluated lines are highly inbred and, thus, these SNPs may come from sequencing errors. Missing data were imputed with the hidden Markov model (HMM) using the Beagle software version 5.0 (Browning et al., 2018). After filtering, a total of 273,775 high-quality SNPs were retained for subsequent analyses ( Supplementary Figure 1 ).

The population genetic structure was inferred using the Bayesian clustering model through the software Structure version 2.3.4 (Pritchard et al., 2000) based on the method described by Evanno et al. (2005). One hundred thousand MCMC (Monte Carlo Markov chain) iterations, a burn-in period of 10,000 iterations, an admixture model, and correlated allelic frequencies were used. Subgroup values (ΔK) between one and ten were tested, with ten independent interactions for each K value. The ideal number of K was determined using the Structure Harvester software version 0.6.92 (Earl, 2012).

PCA was performed using the Tassel software version 5.2.48 (Bradbury et al., 2007). Additionally, the Euclidean distance between the inbred lines was calculated, and the UPGMA (unweighted pair group method using arithmetic averages) hierarchical grouping was performed. These analyses were performed using the R software version 3.6.0 via the packages 'SNPRelate', 'factoextra', and 'ggplot2'. The kinship matrix was calculated based on the centralized identity-by-state (IBS) method (Endelman and Jannink, 2012) using the Tassel software version 5.2.48 (Bradbury et al., 2007).

Linkage disequilibrium (LD) analysis was performed using the 'LDcorSV' package of the R software version 3.6.0. LD was estimated on all chromosomes simultaneously and individually. In addition to the conventional measure of r², corrected r² was estimated by: i) population structure (r²Q), ii) kinship matrix (r² _K), and iii) population structure plus kinship matrix (r²QK) (Mangin et al., 2012). The results of the Bayesian clustering and IBS method were used to correct the population structure and kinship matrix, respectively, as presented above. LD was adjusted using the nonlinear regression method proposed by Hill and Weir (1988) using the nls function of the R software version 3.6.0. The LD decay was defined by the distance at which half of the maximum LD decayed (LD half decay). This estimate indicates the initial slope of the LD decay, and it was considered the most consistent in the comparative study conducted by Vos et al. (2017).

For multi-locus GWAS analyses, least square means (lsmeans) estimated from the lsmeans function of PROC GLM in SAS software version 9.0 (SAS Institute, Cary, USA) were used. Multi-locus GWAS analyses were performed using eight datasets: i) 2018-2019 summer season at low P (LP_18), ii) 2018-2019 summer season under normal P (NP_18), iii) 2019 fall/winter season in low P (LP_19), iv) 2019 fall/winter season in high P (NP_19), v) combination of 2018–2019 and 2019 seasons in low P (LP_C), vi) combination of 2018–2019 and 2019 seasons under normal P (NP_C), vii) greenhouse at low P (LP_G), and vii) greenhouse in normal P (NP_G).

Five multi-locus GWAS methods implemented in the ‘mrMLM.GUI’ package (Zhang et al., 2020) of the R software version 3.6.0 were used: i) mrMLM (Wang et al., 2016), ii) FASTmrMLM (Tamba and Zhang, 2018), iii) FASTmrEMMA (Wen et al., 2018), iv) ISIS EM-BLASSO (Tamba et al., 2017), and v) pLARmEB (Zhang et al., 2017). The population structure (Q) and kinship matrix (K) were included in the tested models to minimize the identification of false-positive associations and increase the statistical analysis power. Critical values for significant associations were LOD (logarithm of odds) ≥ 3 (or P = 0.0002) for all methods. To obtain more accurate results, only the QTNs detected by at least three different methods were considered truly significant and, later, used in the search for candidate genes.

The QTNs detected had their favorable alleles identified, that is, the favorable alleles that cause positive effects on the traits. Then, a heatmap based on Ward's method and Euclidean distance was performed to group the inbred maize lines into clusters with different amounts of favorable alleles. In addition, box plots were made to verify if the pyramiding of these favorable alleles would result in higher PUE and its components. These analyses were performed using R software version 3.6.0. through the 'pheatmap' and 'ggplot2' packages.

Candidate genes were selected based on the reference genome of Zea mays L. version AGPV3 (B73 RefGen v3) obtained from the MaizeGDB database (https://www.maizegdb.org). The search radius of candidate genes was established based on the results obtained from the LD half decay. The classical genes or those with known functions in maize were annotated using the MaizeGDB database. Additionally, the molecular functions of all candidate genes were annotated according to the Gene Ontology (GO) database (http://www.geneontology.org/).

The deviance analyses showed a significant effect (P ≤ 0.01) of genotype and genotype–P level interaction for all traits evaluated in the greenhouse trials ( Table 2 ). Low P reduced all characteristics, except for PUtE_g, PUpE_g, and PUE_g. The reductions in these traits ranged from –4.8 (RD) to –65.0% (PS_g), while a high increase was observed for PUtE_g (161.9%), PUpE_g (169.6%), and PUE_g (678.9%) ( Supplementary Figure 2 ). The estimates of h² ranged from 0.55 (RV) to 0.74 (RD and RSA) at low P, while under normal P conditions, the estimates varied from 0.63 (RV) to 0.86 (RD). Ac values ranged from 0.74 (RV) to 0.86 (RD and RSA) and 0.86 (RV and SB_g) to 0.92 (RD and RSA) under low and normal P conditions, respectively.

In the field trials, the deviance analyses presented significant triple interaction (P ≤ 0.05) between genotype–season–P level for EH, ED, EL, PS_f, HI, PUtE_f, and PUE_f. For genotype–P level interaction, significant effects (P ≤ 0.05) were observed for SB_f, EL, PS_f, HI, PUtE_f, PUpE_f, and PUE_f. All traits assessed showed significant effects of genotype (P ≤ 0.01) and genotype-season interaction (P ≤ 0.05). Low P led to reduced values for all evaluated traits, except PUtE_f, PUpE_f, and PUE_f. The reductions in trait values ranged from –4.8 (EL) to –37.5% (PHI), while PUtE_f (19.4%), PUpE_f (221.9%), and PUE_f (256.7%) increased greatly. The estimates of h² ranged from 0.44 (HI) to 0.72 (PH) under low P and from 0.52 (PS_f) to 0.82 (EH) under normal P conditions. Ac values ​​ranged from 0.66 (HI) to 0.84 (PH) and 0.76 (PS_f) to 0.91 (EH) under low and normal P conditions, respectively.

The correlations among traits in the greenhouse and field are shown in Figure 1 . In general, there was a greater correlation between the traits within than between the field and greenhouse trials. However, few differences were observed between low and normal P conditions, mainly due to the magnitudes of the correlation estimates rather than their direction. Strong positive correlations were observed between greenhouse traits under both P conditions, primarily root system-related traits, PUtE_g, and PUE_g. Under field conditions, high and positive correlation estimates were found for PH × EH, PUE_f × HI, PUE_f × PUpE_f, PUE_f × SB_f, PG × PHI, and PUpE_f × SB_f. Conversely, negative correlations were reported for PS_g × PUpE_g in the greenhouse under low and normal P conditions.

PCA of low versus normal P, low P, and normal P conditions are shown in Figure 2 . Together, the first two principal components (PC1 and PC2) explained 65.2% (low P), 70.5% (normal P), and 72.9% (low versus normal P) of the total existing variation. There was a clear distinction between low versus normal P, with PUE_g, PUE_f, PUtE_g, PUtE_f, PUpE_g, and PUpE_f vectors associated with low P, while the others were associated with normal P ( Figure 2A ). The field corn and popcorn lines showed distinct collective behavior under both P conditions. At low P, the field corn lines were associated with ED, SB_f, PUtE_f, PUpE_f, and PUE_f vectors, whereas the popcorn lines were mainly linked to the PH, EH, PG, and PHI vectors ( Figure 2B ). Under normal P, the field corn lines were associated with the ED, HI, SB_f, PUtE_f, PUpE_f, and PUE_f vectors, while the popcorn lines were primarily linked with the PH, EH, PG, PHI, and PS_f vectors ( Figure 2C ).

Based on the PCA, the first two components explained 19.27 and 2.54% of the total existing variation, respectively, and it is possible to observe a distinction between field corn and popcorn lines ( Figure 3A ). Similar results were found by the UPGMA and Bayesian groupings ( Figures 3B, C ). Bayesian analysis classified seven popcorn lines as admixture, with an ancestry coefficient lower than 0.6 for each subpopulation.

LD showed rapid decay and exhibited variations within each chromosome ( Supplementary Figure 3 ). Considering all chromosomes, the differences between uncorrected r² (r² = 0.16; LD half decay = 1.41 kb) and corrected population structure (r²K = 0.15; LD half decay = 1.27 kb) were small. On the other hand, r² was affected by the corrections of the kinship matrix (r²Q = 0.07; LD half decay = 0.89 kb) and kinship matrix plus population structure (r²KQ = 0.07; LD half decay = 0.89 kb), presenting similar results ( Supplementary Figure 4 ). Regarding individual chromosomes, Chr07 and Chr05 showed the fastest and slowest decay, respectively. From the results obtained from the LD half decay, the search distance for candidate genes was defined as ± 1 kb.

The histograms of the lsmeans used to identify the QTNs are shown in Supplementary Figures 5, 6 . A total of 5838 QTNs was discovered by the five different multi-locus GWAS methods evaluated ( Supplementary Figure 7 ). The pLARmEB method (n = 1710) identified the highest number of QTNs, followed by the FASTmrMLM (n = 1604), ISIS EM-BLASSO (n = 1224), mrMLM (n = 756), and FASTmrEMMA (n = 544) methods. Among these QTNs, 810 were co-detected by two methods, 255 by three, 44 by four, and seven by five. On the other hand, 3242 QTNs were found by only one of the evaluated methods. Thus, only 255 QTNs were considered highly reliable, as they were co-detected by at least three methods and, consequently, were maintained in the present study ( Supplementary Table 3 ).

The FASTmrMLM method was the most efficient (95.4%), as 292 of the 306 QTNs were detected by this method, followed by pLARmEB (94.1%), ISIS EM-BLASSO (85.0%), FASTmrEMMA (29.1%), and mrMLM (19.9%). The 306 QTNs discovered are distributed in the ten chromosomes (Chr) of maize, ranging from 15 (Chr08) to 60 QTNs (Chr01) ( Supplementary Figure 8A ). A greater number of QTNs was detected in the field trials (n = 227) compared to greenhouse trials (n = 79) ( Supplementary Figure 8B ). In the greenhouse, 35 and 44 QTNs were found under low and normal P conditions, respectively ( Supplementary Figure 8C ), whereas 118 (low P) and 107 (normal P) QTNs were detected in the field trials ( Supplementary Figure 8D ).

The 79 QTNs identified in the greenhouse were associated with RV (n = 3), PUpE (n = 4), PUE (n = 5), SB (n = 5), NR (n = 6), RD (n = 6), RL (n = 6), PS (n = 7), TB (n = 7), RSA (n = 8), PUtE (n = 9), and RB (n = 13) ( Supplementary Table 3 ). The phenotypic variation explained (PVE) by these QTNs ranged from 0.01 to 16.31%. In turn, the 227 QTNs identified in the field were associated with PUE (n = 13), EL (n = 14), ED (n = 16), PG (n = 18), PH (n = 18), EH (n = 19), HI (n = 19), SB (n = 19), PUtE (n = 20), PUpE (n = 22), PS (n = 24), and PHI (n = 25) ( Supplementary Table 3 ). The PVE by these QTNs varied from 0.01 to 18.62%.

A total of 16 pleiotropic QTNs were detected in the greenhouse (n = 10) and field (n = 6) trials ( Table 3 ). In the greenhouse experiment under low P, pleiotropic QTNs were associated with SB–TB (S2_170934541), PUtE–RL (S4_198262446), and RL–RSA (S5_97898187 and S9_11783675). Furthermore, a pleiotropic QTN (S1_46130565) was detected in three characteristics (RSA, NR, and RL). Under normal P conditions, pleiotropic QTNs associated with PUtE–PS (S1_1159069), RB–TB (S1_151079747), RSA–TB (S3_225847941), and RB–RSA (S5_212662145) were revealed. Moreover, a highly pleiotropic QTN (S10_119011145) was found in four traits simultaneously (PUE, RB, SB, and TB). In the field trials under low P, pleiotropic QTNs were associated with PUE–PUpE (S3_190210784) and HI–PUE (S8_148017105). Under normal P conditions, pleiotropic QTNs were identified in PG–PH (S1_259083391), HI–PUE (S5_10594899), and ED–PUE (S7_7624808).

The heatmap grouped the inbred maize lines into three distinct groups regarding favorable allele number ( Supplementary Figure 9 ). Group 1 was composed of 41 lines, characterized by having the smallest number of favorable alleles in the field experiments, mainly under low P conditions. Groups 2 and 3 were formed by 36 and 55 lines, respectively. In general, the lines in group 2 had the highest numbers of favorable alleles in the greenhouse and field trials in both P conditions. On the other hand, group 3 lines had the lowest numbers of favorable alleles, mainly in the greenhouse experiments. In general, there was no distinction between popcorn and field corn lines concerning the accumulation of favorable alleles.

Considering only the favorable alleles related to PUE and its components, a gradual increase in these traits was observed by pyramiding these favorable alleles ( Figure 4 ). Popcorn lines 17-P9-1-6 and 162-P1780 had the highest and lowest number of favorable alleles considering all traits, respectively ( Figure 5 ). The 17-P9-1-6 line, considered P-efficient, presented 58.33% (field in normal P), 57.14% (field in low P), 62.16% (greenhouse in normal P), and 58.27% (greenhouse under low P) of the favorable alleles. On the other hand, line 162P1780, considered P-inefficient, presented only 32.18% (field in normal P), 31.45% (field in low P), 27.02% (greenhouse in normal P), and 24.13% (greenhouse under low P) of the favorable alleles.

A total of 186 potential candidate genes distributed in all maize chromosomes were identified ( Supplementary File 1 ). Among them, 45 have already been annotated as classical or known genes in the maize genome, 21 of which have functions related to transcription factor activity, while the others are mainly associated with transport and transferase activity. According to the GO annotation, the primary molecular functions of the candidate genes included transcription factor, transferase, hydrolase, catalytic binding, ATP binding, DNA binding, and nucleotide-binding activity. In addition, several molecular functions directly related to P molecules were also identified, such as activities of phosphotransferase, phosphatase, pyrophosphatase, calcium-dependent phospholipid binding, all-trans-nonaprenyl-diphosphate synthase, glucose-6-phosphate isomerase, phosphorelay sensor kinase, phosphorylative mechanism, phospholipid binding, and phosphoenolpyruvate carboxylase.

Genotype–season–P level interaction was observed in most of the traits evaluated in the field, indicating that the differential behavior of genotypes depends on the combination of P levels and seasons studied. Additionally, there were more significant interactions between genotype–season compared to genotype–P level interactions, suggesting that season strongly influences the evaluated traits. In the greenhouse trials, there was a significant genotype-P level interaction for all characteristics, indicating that the differential behavior of genotypes also depends on the P level in early evaluations. The genotype–P levels interaction in maize has already been reported in several greenhouse and field studies, which confirms that studies related to PUE and its components should be carried out under specific P conditions (Xu et al., 2018; Sahito et al., 2020; Li et al., 2021).

The negative influence of low P on the characteristics evaluated in the greenhouse and field was already expected since P is an essential element for plant development, constituting key cell molecules such as ATP, nucleic acids, and phospholipids, in addition to acting in the central regulation of many metabolisms, including energy transfer, protein activation, and carbon and amino acids metabolic processes (Wang et al., 2021). In contrast, the lines showed higher PUE, PutE, and PUpE under conditions of low P availability. In general, plants present higher PUE under low P since, under high phosphate fertilization, part of this element is lost in erosive processes and/or made unavailable by adsorption and immobilization processes in the soil (Alewell et al., 2020).

Although PUE-related traits are complex and strongly influenced by the environment, the heritability and selective accuracy estimates obtained in the present study fell into moderate to high, indicating favorable conditions for improvement aimed at increasing PUE. Heritability is the central parameter of any breeding program, used to estimate selection response and explain the proportion of phenotypic variation due to genetics (Hallauer et al., 2010). In the present study, heritability estimates were higher under normal P than low P conditions. In general, stress conditions tend to reduce heritability estimates, mainly due to the increased environmental influence on plant phenotype (Vats, 2018). Xu et al. (2018), evaluating 11 PUE-associated traits in maize, also observed reductions in heritability estimates under low P compared to normal P conditions.

The existence of correlations among traits evaluated in the greenhouse and field trials would enable the early selection of genotypes with higher PUE. Unfortunately, however, no relevant correlations were observed between the traits evaluated in both trials, indicating that early selection in hydroponic systems may not be efficient in selecting plants with higher PUE under field conditions. Similar results were reported in other studies, suggesting that different genes and morphophysiological mechanisms act at different stages of plant development (Dissanayaka et al., 2018; Wang et al., 2019a). In fact, experiments in hydroponic systems or pots using nutrient solutions may present conditions that are very different from real field conditions since the environmental influence is much greater in the field than under greenhouse-controlled conditions (Araus et al., 2018).

P is relatively immobile in the soil, and its availability is greater at surface soil horizons, decreasing substantially at deeper soil horizons (Menezes-Blackburn et al., 2018). Thus, several studies have reported the importance of the maize root architecture system for higher PUE under field and greenhouse conditions (Jia et al., 2018; Klamer et al., 2019). The present study corroborated these results, showing a high correlation between the root-system-related traits (RB, RSA, RV, NR, and RL) with PUE in both P conditions under a hydroponic system. The indirect selection of PUE through root system traits is highly relevant in breeding programs, representing a faster and money-saving alternative, as it does not require laboratory analysis for P quantification. Further, high-yield phenotyping methods are already well established to characterize the root system of plants (Tracy et al., 2020).

Linkage disequilibrium (LD) is the non-random association of alleles at two or more loci and plays a central role in association mapping analyses (Romay et al., 2013). In the present study, we observed a rapid LD decay. In addition, we also report differences in LD decay in relation to chromosomes and correction methods assessed, with LD decay being mainly affected by population structure. Several studies have reported rapid LD decay in maize, ranging primarily from 0.1 to 10 kb (Yan et al., 2009; Romay et al., 2013; Bennetzen et al., 2018; Coan et al., 2018). In addition to the different approaches to estimating LD and LD decay, LD estimates may vary depending on the genomic region and germplasm evaluated (Vos et al., 2017). In maize, tropical and subtropical germplasms have faster LD decay when compared to temperate germplasms, as they have greater genetic diversity and the presence of rare alleles (Yan et al., 2009; Romay et al., 2013).

Multi-locus GWAS methods have been recently used to investigate the genetic basis of essential traits in several agricultural crops, such as maize (Lu et al., 2021), rice (Cui et al., 2018), bread wheat (Yang et al., 2020), soybean (Zhang et al., 2018), barley (Hu et al., 2018), and upland cotton (Su et al., 2018). In the present study, five different multi-locus GWAS methods were assessed. Among them, the FASTmrMLM and pLARmEB were the most efficient methods, as they detected more than 94% of the reliable QTNs in the present study, that is, QTNs co-detected by three or more methods. Although these methods present a combined two-step approach, the quantitative and qualitative differences in the identified QTNs are due to the different algorithms used in each method (Zhang et al., 2020). Thus, Zhang et al. (2019) suggested using combined results between the different methods of multi-locus GWAS to obtain more accurate results.

The present study identified a large number of small effects QTNs, confirming the complex and quantitative nature of PUE in maize. Using small-effect QTNs associated with traits of interest represents a smart strategy in genomic selection (GS) approaches since the use of these QTNs alone can replace the need for high-density genotyping through random SNPs, thus reducing genotyping costs (Lan et al., 2020). Furthermore, GS models using only SNPs known to be associated with traits of interest showed greater prediction accuracy, as they had lower background noise (Ali et al., 2020). In addition to their use in GS, these QTNs may be promising for application in breeding programs aimed at pyramiding favorable alleles through marker-assisted selection (MAS) (Boopathi, 2020).

In a QTLome study for low P tolerance in maize under limiting phosphorous conditions, Zhang et al. (2014) identified 23 meta-QTLs, that is, genomic regions more likely to harbor candidate genes responsible for traits related to tolerance to low P and/or PUE. Among the 306 QTNs discovered in the present study, 20 are located within the ranges of 11 meta-QTLs. In another analysis of meta-QTLs, Guo et al. (2018) identified 53 meta-QTLs associated with the maize root architecture system in the presence and absence of abiotic stresses (water deficit, high temperature, and P and N deficiency). Among the 54 QTNs found in the present study to be associated with root traits, 10 are located within the ranges of eight meta-QTLs. Thereby, the location of QTNs in genomic regions known to be associated with traits and characteristics of interest reinforces the significance and relevance of the QTNs pinpointed in the present study.

Among the 186 candidate genes detected in the present study, 45 have already been annotated as classical or known genes in maize. For instance, the GRMZM2G083841 gene (associated with PUpE in LP_G) codes for a phosphoenolpyruvate carboxylase, catalyzing the carboxylation of phosphoenolpyruvate to produce oxaloacetate and Pi, performing the primary fixation of atmospheric CO₂ in C4 plants (Nimmo, 2003). The GRMZM2G140614 gene (associated with PUE in NP_C) encodes a glucose-6-phosphate isomerase, an essential enzyme in the oxidative pathway of pentose phosphate responsible for producing energy-rich cofactors and generation of carbon skeletons required for biosynthetic reactions (Yu et al., 2000). Interestingly, GRMZM2G307119 (associated with RV in NP_G) was initially related to the formation of spikelet meristems (Chuck and Hake, 2005). However, Jansen et al. (2013) revealed that this gene plays a fundamental role in developing lateral roots in maize. In turn, the GRMZM2G177792 gene (associated with PG in LP_C) has a peroxidase activity and has been associated with resistance to biotic stresses (Santiago et al., 2016; Musungu et al., 2020). In addition to these genes, three main mechanisms and pathways were discovered: