Trunk circumference (TC), an indicator of tree vigor, was measured using a tape measure positioned 10 cm above the graft union. Due to differences in planting times, the trees were assessed at different ages, resulting in unbalanced data. We used TC data for trees assessed at the ages of 9 (TC9, 200 unique accessions) and 12 (TC12, 199 unique accessions) years (total of 207 unique accessions) due to the availability of a relatively large number of individuals assessed in these years.

Physiologically mature fruits were harvested from the outer tree canopy, where they were exposed to sunlight. The fruits were washed thoroughly with a detergent, treated with a fungicide dip (1.0 ml L-1 Fludioxonil (230g/L)) for five minutes at 52 °C to control anthracnose. They were then stored in a ripening room maintained at 22°C until they developed a soft texture. Fruit blush color (FBC) was assessed in 220 accessions over at least two seasons, using ten ripe fruits from each accession. FBC of the ripened fruit was rated on a categorical scale, in order from least to most desirable: no blush, orange, pink, pink-red, red, and burgundy. FBC categorical data was converted to a numerical scale as: no blush or yellow = 0, orange = 1, pink = 2, pink-red = 3, red = 4, and burgundy = 5.

The average fruit weight (AFW) in grams (g) was calculated across 222 accessions using the weight of ten fruits at the eating ripeness stage. Fruit firmness (FF) was measured in 221 mango accessions using an analogue firmness meter. Not all accessions were assessed for the three fruit quality traits in every season due to the irregular bearing of some cultivars and differences in planting seasons, resulting in unbalanced data.

Two approaches were used to validate genomic prediction models in this study. In the first cross-validation approach (parental validation), own phenotypes of the 41 gene-pool accessions that are being used as parents in the QMBP served as an independent dataset for model validation, while the remaining gene-pool accessions served as the training population. Predictive ability was estimated as the Pearson correlation between the phenotypes predicted by the linear mixed models (GEBVs) and the observed phenotypes of parental accessions, r ( y ,     ^ y ) .

To provide a more robust evaluation of model performance, a second validation approach involving random 5-fold cross-validation (5-fold CV) was also implemented. In this approach, the entire gene pool collection was randomly partitioned into five subsets in which each subset consisted of 20% of the accessions. For each fold, four subsets (80% of total individuals) were used for model training and the remaining fold (20% of the accessions) for model validation. Predictive ability was calculated as the Pearson correlation between the GEBVs and the observed phenotypes after each 5-fold CV run. To ensure stability and reliability of the predictive ability estimates, the 5-fold CV procedure was repeated five times. Thus, 25 correlation values were calculated for each model. For trunk circumference, only phenotypic data collected from trees aged 12 years were used for validation. The bias of predictions was calculated as the regression of phenotypes on GEBVs for individuals in the validation set.

Genomic predictive ability varied across traits, marker density, and validation strategy ( Table 2 , Supplementary Table 7 ), and generally aligned with the narrow-sense heritability estimates ( h ^ 2 ). When considering predictions based on WGS data and baseline GBLUP models (i.e., models without population structure correction or fixed-effect SNPs), higher predictive abilities (PA) were observed for highly heritable traits and lower predictive abilities for traits with lower h ^ 2 . Under the parental validation strategy, the highest predictive abilities were observed for FBC and AFW (PA = 0.67 for both traits), followed by TC (PA = 0.54), with FF showing the lowest predictive ability (PA = 0.41). A similar trend was observed in the 5-fold cross-validation (CV) strategy, where predictive abilities for AFW (0.65) and TC (0.57) were comparable to those from the parental validation ( Supplementary Table 7 ). However, the predictive ability for FBC increased substantially under the 5-fold CV strategy (0.80), while that for FF decreased markedly (0.28), relative to the parental validation results.

Results from the evaluation of marker density effects under the parental validation strategy using baseline GBLUP models revealed that predictive ability varied with density. Predictive ability ranged from 0.60 to 0.67 for FBC, 0.59 to 0.67 for AFW, 0.35 to 0.41 for FF, and 0.51 to 0.54 for TC ( Supplementary Table 5 ). Across all traits, predictive ability generally increased with marker density but plateaued beyond LD_20k (~20,000 SNPs), indicating little gains at higher SNP densities. Models incorporating a GRM estimated from the lowest-density marker set (LD_10k) exhibited substantially lower predictive ability compared to those using higher-density marker sets (LD_20k to WGS), which showed only marginal variation in predictive ability among themselves. For TC, differences in predictive ability were relatively stable across marker densities, with a maximum difference of just 0.03 between LD_10k and WGS. Under the 5-fold CV, differences in predictive ability across marker densities were minimal for all traits ( Supplementary Table 7 ).

Incorporating the top six principal components (PCs) as fixed effects to account for population structure resulted in substantial reductions in predictive ability for all traits ( Table 2 ). The decrease in predictive ability ranged from 0 to over 100% depending on marker set and validation approach employed, highlighting the dominant influence of population structure on genomic prediction within this gene-pool for these traits. Under parental validation, predictive ability based on WGS data decreased from 0.67 to 0.44 for FBC, from 0.67 to 0.30 for AFW, and from 0.41 to 0.30 for FF when population structure was accounted for ( Supplementary Table 6 ). The exception was TC, where a slight increase in predictive ability was observed, rising from 0.54 to 0.57. Notably, while population structure correction reduced predictive ability for FBC, this decline was substantially mitigated when the FBC-associated SNP on chromosome 15 was fitted as a fixed effect in GBLUP models. When only the top six PCs were included as fixed effects, predictive ability for FBC dropped by 34%. However, when both the first six PCs and the most significant GWAS-identified SNP were jointly fitted as fixed effects, the reduction in predictive ability was mitigated to just 7%.

Similarly, results from 5-fold cross validation revealed a marked decline in predictive ability after correcting for population structure ( Supplementary Table 8 ). However, unlike in the parental validation strategy where the predictive ability for TC remained stable despite population structure correction, the predictive ability in the 5-fold cross-validation declined sharply, dropping from 0.57 to 0.45 when using WGS data.

Reliable trait-associated SNPs (identified by at least two GWAS methods) showed clear effects on phenotypic variation, as revealed by GEBVs for the three genotypic classes: homozygous reference, heterozygous, and homozygous alternate allele ( Supplementary Figure 5 - Supplementary Figure 7 ). For FBC, the SNP on chromosome 15 (G/A) showed that cultivars with the GG genotype (e.g., ‘Ah Ha!’, ‘Tommy Atkins’, and ‘Irwin’) had significantly higher FBC ratings (p < 0.0005, mean GEBV = 2.0) compared to those carrying the A allele either in homozygous form (mean GEBV = 1.2; e.g., ‘Dashehari’, ‘Mallika’, and ‘Arumanis A’) or heterozygous form (mean GEBV = 1.0; e.g., ‘Maha Chanook’, ‘Alphonso’, and ‘Carabao Pep’). For AFW, the SNP on chromosome 17 (A/G) revealed that cultivars with the A allele in homozygous form had significantly lower fruit weight (p< 0.0005; mean GEBV = 322.0 g) than the heterozygous cultivars (mean GEBV = 415.2 g). For TC, the SNP on chromosome 7 (T/A) revealed that AA genotypes (e.g., ‘Manjeera’, ‘Lippens’) had significantly lower trunk circumference (p < 0.0005; mean GEBV = 45.5 cm) compared to cultivars carrying the T allele either in homozygous form (mean GEBV = 56.5 cm) or heterozygous form (mean GEBV = 52.4 cm). Notably, heterozygous (T/A) genotypes also had significantly smaller trunk circumference (p < 0.0005) than homozygous TT genotypes.

Models incorporating a GRM derived from variants preselected based on the highest ranked probability of effect as estimated using GWAS improved predictive ability across all traits, with improvements of up to 93% under parental validation ( Table 2 ). The magnitude of these improvements varied depending on the trait, density of GWAS-preselected variants, GWAS method applied, and whether population structure was accounted for. When using base models (i.e., models without population structure correction or fixed-effect SNPs) under the parental validation strategy, preselecting variants based on GWAS showed an advantage depending on the GWAS method used to identify variants, particularly for AFW and to a lesser extent for FBC, FF, and TC ( Figure 2A ). The predictive ability for AFW was markedly higher when using 20,000 TOP-BLINK GWAS-preselected variants, reaching 0.78, compared to 0.67 using the complete WGS dataset. In contrast, improvements in predictive ability for other traits were more modest, increasing from 0.67 to 0.71 for FBC using 100,000 SNPs from the TOP-FarmCPU set, from 0.54 to 0.59 for TC using either 20,000 or 50,000 SNPs from the TOP-BLINK or TOP-FarmCPU set, and from 0.41 to 0.45 for FF using 1,000 SNPs from the TOP-GLMM set. However, under 5-fold cross-validation using models that did not account for population structure, GWAS-based SNP preselection did not lead to improvements in predictive ability across any of the traits ( Figure 3A , Supplementary Table 7 ).

The increases in predictive ability observed with GWAS-preselected variants relative to WGS data were much larger when population structure was accounted for ( Figure 2B ). Under the parental validation strategy, adjusting for population structure in GBLUP models led to a 93% improvement in predictive ability for AFW, increasing from 0.30 to up to 0.58 when 15,000 variants from the TOP-FarmCPU set were used instead of WGS data. Similar improvements in predictive ability were observed for FBC and FF, rising from 0.44 to 0.50 using either 50,000 or 100,000 TOP-FarmCPU SNPs for FBC, and from 0.30 to 0.35 using top 1,000 SNPs selected by GLMM for FF. In contrast, there was little variation in predictive ability for TC between models that included GWAS pre-selected variants with or without adjustment for population structure.

A comparable pattern was observed under the 5-fold cross-validation strategy in GBLUP models that included population structure correction ( Figure 3B ). Specifically, predictive ability increased by up to 29% for FBC (from 0.35 to 0.45), 50% for AFW (from 0.32 to 0.48), and 150% for FF (from 0.08 to 0.20), while TC showed a modest improvement of 11% (from 0.45 to 0.50). The highest predictive abilities under 5-fold cross validations were achieved using 10,000 SNPs from TOP-GLM for FBC, 1,000 SNPs from TOP-MLMM for AFW, 1,000 SNPs from all GWAS methods for FF, and 10,000 or more SNPs from either TOP-FarmCPU or TOP-GLMM for TC ( Supplementary Table 8 ). Notably, for FF and TC, the predictive abilities obtained using GWAS-preselected variants were comparable to those achieved using LD-pruned marker sets (LD_10k to LD_2mil).

Predictive abilities using GWAS-preselected variants showed substantial variation depending on marker density and validation strategy, with no consistent trend across traits ( Supplementary Tables 5 , 7 ). Under the parental validation strategy in models ignoring population structure, the highest predictive abilities were achieved using 20,000 SNPs from BLINK for AFW, 100,000 SNPs from FarmCPU for FBC, 1,000 SNPs from GLMM for FF, and either 20,000 or 50,000 SNPs from BLINK or FarmCPU for TC. In contrast, under 5-fold cross validation, predictive abilities remained relatively stable across different marker densities ( Supplementary Table 7 ). Differences in maximum predictive ability between GWAS models were generally small (< 0.03), except for FBC and AFW in models that accounted for population structure ( Supplementary Table 8 ). In these cases, the highest predictive abilities were achieved using 10,000 and 1,000 SNPs from TOP-GLM and TOP-MLMM, respectively. All subsequent results are based on the parental validation strategy using both the full WGS dataset and the optimal set of GWAS-preselected variants for each trait.

The impact of incorporating reliable markers as fixed effects on predictive ability varied depending on the trait, marker set, and whether population structure was accounted for ( Figure 4 , Supplementary Table 9 , Supplementary Table 10 ). Our findings indicate that incorporating a reliable SNP as a fixed effect in prediction models markedly improved predictive ability for FBC and TC, with gains of up to 0.26 and 0.07, respectively.

For FBC, incorporating the reliable trait-associated SNP on chromosome 15 as a fixed effect in GBLUP models without accounting for population structure resulted in an improvement in predictive ability ranging from 0.01 to 0.07 compared to models without the fixed-effect SNP. Notably, predictive ability increased from 0.67 to 0.74 with WGS data and from 0.70 to 0.77 using 100,000 TOP-BLINK markers when the FBC-reliable marker was included as a fixed effect. Strikingly, under population structure correction, the enhancement in predictive ability due to the inclusion of the FBC-reliable marker as a fixed effect was even more pronounced, with gains ranging from 0.12 to 0.26. In these population structure corrected models, predictive ability increased from 0.44 to 0.62 with WGS data, from 0.45 to 0.69 using 50,000 TOP-BLINK markers, from 0.50 to 0.62 using either 50,000 or 100,000 TOP-FarmCPU markers, from 0.37 to 0.51 using 15,000 TOP-MLMM markers, and from 0.33 to 0.48 using 10,000 TOP-GLMM markers. Further analysis using ~ 2 million SNPs showed that this FBC-associated SNP accounted for 36% of the genetic variance (results not shown).

Incorporating the reliable TC-associated SNP on chromosome 7 as a fixed effect in GBLUP models also improved predictive ability, with gains of up to 0.07 in models without population structure control, and up to 0.06 when population structure was accounted for. For example, predictive ability increased from 0.54 to 0.58 with WGS data, from 0.59 to 0.64 using 20,000 or 50,000 SNPs from either BLINK or FarmCPU, and from 0.54 to 0.61 with 20,000 TOP-MLMM markers when the reliable TC-associated SNP was included as a fixed effect in models without population structure correction. A similar pattern was observed when population structure was accounted for, with the predictive ability for WGS data increasing from 0.57 to 0.61, from 0.61 to 0.66 using 50,000 SNPs from either BLINK or FarmCPU, and from 0.55 to 0.62 using 15,000 TOP-MLMM markers. In contrast, for AFW, adding fixed-effect markers to the prediction models did not improve predictive ability.

Combining GWAS-preselected variants with fixed-effect SNPs substantially improved predictive ability for FBC and TC compared to models using WGS data alone or GWAS-preselected variants alone. The highest predictive abilities for these traits were achieved using this integrated approach both with and without population structure correction ( Table 2 ; Figure 4 ). For example, substituting WGS data with TOP-BLINK markers improved predictive ability for FBC from 0.67 to 0.70 ( Figure 2A ). Incorporating the FBC-associated reliable SNP on chromosome 15 as a fixed effect increased predictive ability with WGS data from 0.67 to 0.74 (a 0.07 increase). Notably, combining 100,000 GWAS-preselected variants from BLINK with the fixed-effect SNP yielded a substantial improvement, boosting predictive ability by 0.10 (from 0.67 with WGS data to 0.77 using a combination of GWAS-preselected variants and the fixed-effect SNP). A similar trend was observed when population structure was accounted for, with the highest predictive ability (0.69) achieved by incorporating the FBC-reliable SNP as a fixed effect in a GBLUP model based on a GRM derived from 50,000 TOP-BLINK markers. This predictive ability represents a substantial improvement, exceeding that of WGS data alone by 0.25 and that of TOP-BLINK markers alone by 0.24, and surpassing WGS data with a fixed-effect SNP by 0.07.

A similar trend was observed for TC, where the highest predictive abilities were achieved by integrating GWAS-preselected variants with fixed-effect SNPs in a single GBLUP model. Specifically, including 20,000 or 50,000 GWAS-preselected variants from TOP-FarmCPU or TOP-BLINK alongside fixed-effect SNPs improved predictive ability by 0.1, increasing from 0.54 with WGS data to 0.64. This enhancement in predictive ability surpasses the gains of 0.06 and 0.05 obtained when using either WGS data plus fixed-effect SNPs or GWAS-preselected variants alone. Notably, a comparable pattern emerged in GBLUP models that accounted for population structure, with predictive abilities remaining nearly identical to those observed in models without population structure correction.

When considering the optimal marker density for GWAS-preselected variants under parental validation, defined as the density yielding the highest predictive ability, TOP-BLINK and TOP-FarmCPU both achieved the highest predictive abilities in eight of the fourteen trait-by-scenario combinations (four traits and four scenarios [GBLUP, GBLUP + fixed SNP, GBLUP + fixed PCs, and GBLUP + fixed PCs + fixed SNP]). TOP-GLMM produced the highest predictive ability in three combinations, while MLMM did not result in the highest predictive ability in any of the scenarios. In contrast, the performance of GWAS models under 5-fold cross-validations was comparable across all traits when population structure was ignored. However, under models that accounted for population structure, GWAS-preselected variants identified using the MLMM and GLMM yielded the highest predictive ability for AFW and FBC, respectively.

This study identified five distinct and statistically significant associations for FBC ( Table 3 , Supplementary Figure 2 ). Notably, a MYB114 transcription factor was located just 0.5 kb from a key FBC-associated marker on chromosome 15, consistently identified by three different GWAS methods. MYB transcription factors are widely reported to regulate fruit skin color in multiple fruit tree species, including mango (Kanzaki et al., 2020; Wilkinson et al., 2025), apple (Plunkett et al., 2019; Sun et al., 2021), pear (Cong et al., 2021; Zhang et al., 2021) and kiwifruit (Ampomah-Dwamena et al., 2019). These genes are central regulators of the anthocyanin biosynthesis pathway, which plays a critical role in pigmentation of fruit peels (Gao et al., 2021). The well-established role of anthocyanin accumulation in contributing to red skin coloration in fruits is consistent with previous findings in mango. Wang et al. (2020) demonstrated that anthocyanin biosynthesis genes were significantly upregulated in the peel of red-skinned mango cultivars compared to yellow- or green-skinned types. Additionally, Kanzaki et al. (2020) reported that exposure to light stimulus increased the expression of MiMYB1 and MiMYB4 transcription factors in reddened mango fruit, further highlighting the involvement of MYB transcription factors and light exposure in regulating peel coloration.

This study identified 11 novel SNPs significantly associated with AFW ( Table 3 , Supplementary Figure 3 ), a key trait for influencing consumer appeal and market value, and therefore a major target for improvement in mango breeding programs (Bally et al., 2021). Our results support the role of hormone-mediated cell division in determining fruit weight, consistent with findings in other horticultural fruit tree species (Karim et al., 2022; Li et al., 2024; Zhang et al., 2006). Notably, two auxin response factors were identified within 12 kb of an AFW-associated SNP on chromosome 13, suggesting a likely regulatory role of auxin signaling in fruit weight variation. Auxin response factors have previously been implicated in apple fruit weight variation through modulation of cell division and expansion (Devoghalaere et al., 2012). Additionally, an AFW-associated SNP on chromosome 7 was located ~110 kb from two cell division control proteins, reinforcing the mechanistic link between cell division during mango fruit development and fruit size. Similar associations have been reported in sweet cherry (Prunus avium L.), where a fruit size QTL was closely linked to a gene governing cell number, underscoring the conserved nature of these genetic mechanisms across species (De Franceschi et al., 2013).

We identified eight unique marker-trait associations for TC across seven chromosomes ( Table 3 , Supplementary Figure 4 ). Our analyses identified a GATA transcription factor located within 70 kb of a TC-associated SNP on chromosome 7. GATA transcription factors have been reported to regulate tree growth in Populus (An et al., 2020). BLAST analysis revealed that the GATA transcription factor identified in our study shares 86% sequence similarity with the one reported in Populus (An et al., 2020), suggesting similar regulatory mechanisms in mango tree growth. GATA transcription factors are known to modulate the expression of auxin efflux carrier genes, facilitating the basipetal movement of auxins to the roots (An et al., 2020, An et al., 2014). In our study, two auxin efflux carrier genes were located just 6 kb and 16 kb from TC-associated SNPs on chromosome 2, further supporting the potential regulatory role of auxin transport in mango tree growth.

Prior studies strongly support a model in which plant dwarfism results from reduced expression of PIN genes (auxin efflux carriers) in stem bark tissues, leading to impaired auxin transport to the roots. This disruption limits root growth and cytokinin biosynthesis, ultimately constraining shoot development (An et al., 2017; Li et al., 2018). These mechanisms align with previous studies in apple, where use of dwarfing inter-stock (M9) led to decreased expression of auxin efflux carrier genes in stem bark tissues, suppressing the basipetal movement of auxins and leading to reduced root and shoot development (Zhang et al., 2015). Similar findings in pear demonstrated significantly higher expression levels of the PcPIN-L auxin efflux carrier gene in standard-size trees compared to dwarf types (Qi et al., 2020), further reinforcing the role of auxin transport in tree growth regulation.

In addition to the GATA transcription factor and auxin efflux carrier proteins, we identified a growth-regulating factor gene located approximately 21 kb from a TC-associated SNP on chromosome 2. This, together with the proximity of auxin efflux carrier genes and the GATA transcription factor to TC-associated SNPs, suggests the involvement of a coordinated regulatory network governing tree growth in mango. The trait-associated markers identified for FBC, AFW and TC in this study represent a valuable resource for marker-assisted breeding in mango, pending validation in independent populations.

Our findings underscore the superior statistical power of multi-locus GWAS methods compared to single-locus approaches. Consistent with previous studies (Cebeci et al., 2023; Huang et al., 2019; Minamikawa et al., 2018), multi-locus GWAS methods, particularly BLINK and FarmCPU, identified more significant marker-trait associations than the single-locus GWAS approach (GLMM). The increased power of multi-locus GWAS stems from their ability to account for LD between SNPs (as in BLINK) while simultaneously testing multiple markers, enhancing the detection of small-effect loci associated with a trait (Segura et al., 2012; Wang et al., 2016). BLINK and FarmCPU, which use multi-locus strategies and iterative inclusion of pseudo-QTNs, tend to capture both large- and small-effect loci more robustly, especially in polygenic traits. In contrast, MLMM’s stepwise regression approach appears to underperform, likely due to over-adjustment for population structure and the inherent sensitivity of its sequential covariate inclusion, which may mask genuine signals.

Our results, particularly for AFW where BLINK and FarmCPU identified more significant marker-trait associations, highlight the value of this multi-method GWAS strategy. These findings are consistent with reports by Minamikawa et al. (2018) and Kumar et al. (2019), who identified a higher number of trait-associated SNPs in pears (Pyrus pyrifolia) by employing multiple GWAS methods rather than relying on a single approach. Integrating results across multiple GWAS methods is a powerful strategy to identify additional marker-trait associations as no single method is optimal for all traits. Moreover, loci detected by the different methods do not completely overlap (Zhou et al., 2023). The use of a combination of complementary GWAS methods not only strengthens statistical robustness but also strengthens confidence in associations consistently detected across analyses, making these associations strong candidates for marker development and functional validation.

In our study, we observed that increasing marker density beyond a certain threshold, even up to WGS level, did not yield further improvements in predictive ability ( Table 2 ). These findings are consistent with previous studies (Bedhane et al., 2021; Moghaddar et al., 2019; Raymond et al., 2018b, Raymond et al., 2018c; van Binsbergen et al., 2015; VanRaden et al., 2017) that also found little or no improvement in prediction accuracy when using WGS variants compared to lower density or high-density SNP chips. A plausible explanation is that WGS data include many variants that are not in strong LD with the causative loci (van Binsbergen et al., 2015). These non-informative markers may not capture the QTL effects or accurately reflect genetic relationships at causal loci, potentially undermining the performance of GP models through over-shrinkage of QTL effects. This phenomenon likely reflects the balance between capturing true causal variation and overfitting to random, non-informative variation. Our analyses using different LD-pruned subsets (e.g. LD_2mil, LD_800k, etc.) indicated that predictive ability tended to plateau or even decline when the number of markers exceeded an optimal threshold. This threshold is inherently linked to the underlying LD structure and genetic architecture of the trait in question.

This study showed that GRMs constructed using GWAS-preselected variants resulted in higher predictive abilities across the four studied traits compared to GRMs built using all WGS variants ( Table 2 , Figures 2 , 3 ). These findings highlight that preselecting WGS markers likely to be in LD with causal mutations, while excluding those that do not capture genetic relationships at causal loci, can improve genomic predictive ability. Thus, it appears that including markers not in LD with causative mutations in GRM construction may cause the realized genetic relationships to diverge from true relationships at causal loci, thereby reducing the performance of GBLUP models. However, when markers preselected for their potential causal effects are used, the GRM is dominated by SNPs in high LD with QTL for the target trait. Thus, the trait-specific GRM may better capture the genetic relationships among individuals at unobserved causal loci, potentially enhancing the accuracy of genomic predictions. Our results are consistent with those of Tan and Ingvarsson (2022) who showed that when the top 1% of markers from GWAS are selected, the accuracy of genomic predictions can be increased significantly. Chen et al. (2023) also showed that performing GP using a GRM built using 100 preselected markers resulted in improved prediction accuracies compared to models based on all markers.

While our results clearly demonstrate that the integration of GWAS-preselected variants improves predictive ability, we acknowledge that validation confined to a single, relatively small dataset may limit the external applicability and generalizability of our findings. Such internal validation alone does not adequately account for potential biases introduced by population-specific genetic structure or unique environmental factors. Although we employed a 5-fold cross-validation strategy to strengthen robustness of our model assessment, external validation in large, independent datasets such as a full-sib population remains essential. Such validation would verify whether the observed improvement in predictive performance genuinely reflects enhanced capture of causal genetic variation.

While the use of GWAS-preselected variants increased genomic predictive ability in our analyses, this approach still suffers from the assumption of the GBLUP model that all markers contribute an equal and individually small proportion of the total genetic variance (Meuwissen and Goddard, 2010). However, increasing evidence supports the hypothesis that SNPs in high LD with causal mutations explain more genetic variance than those in low LD (Meuwissen et al., 2024). Incorporating fixed-effect SNPs into GBLUP models appeared to improve predictive ability for both FBC and TC, likely by capturing variation associated with major QTLs ( Figure 4 ). This strategy enabled us to account explicitly for the effects of markers with large estimated effects, potentially helping to separate their contribution from those assumed under the infinitesimal model. While these results suggest benefits from including such markers, it remains important to recognize that the identified SNPs may not represent true causal variants, and further validation in an independent population such as a full-sib family would be needed to confirm their functional significance. The differentiation between large- and small-effect QTLs appears to model better the true genetic architecture of traits, leading to more accurate prediction models. This is especially true when markers in LD with major genes are treated as fixed effects (Li et al., 2019). Our findings are consistent with prior studies. For example, Kostick et al. (2023) demonstrated a substantial improvement in the predictive ability of ‘percent red overcolor’ in apple, which increased from 0.33 to 0.80 upon inclusion of a fixed-effect SNP at a fruit color locus. Similarly, Nsibi et al. (2020) reported a 25.8% increase in prediction accuracy for apricot (Prunus armeniaca) fruit color (hue angle) after incorporating two major QTLs as fixed effects.

Critically, the effectiveness of using fixed-effect SNPs relies on their LD with a QTL, as reported by Li et al. (2019). In this study, the fixed-effect SNPs that enhanced predictive abilities were consistently identified by three GWAS methods (reliable SNPs), strengthening the evidence that these SNPs are likely in LD with underlying QTLs.

In our study, we demonstrated that while the utilization of GWAS-preselected variants or fixed-effect SNPs can enhance predictive ability, further improvements can be achieved through the integration of preselected variants with fixed-effect SNPs ( Table 2 , Figure 4 ). Traditional GBLUP models employing a single GRM constructed from GWAS-preselected variants do not fully capitalize on the predictive potential of large-effect SNPs due to the inherent assumptions of the infinitesimal model, which overly constrains their contribution to the total genetic variance. By contrast, our approach, combining preselected variants and fixed-effect SNPs, benefits from more accurate estimation of genomic relationships at causative loci. If all markers explain the same proportion of the total genetic variance, as is the assumption of the infinitesimal model, there would be no notable reduction in heritability when significant SNPs from GWAS are fitted as fixed effects in GBLUP models. However, our analyses demonstrated a notable reduction in additive genetic variance due to the anonymous markers and heritability when the fixed-effect SNP for FBC was included in GBLUP models, suggesting that a substantial portion of the additive genetic variance was explained by this SNP potentially due to its LD with the causative mutation. For mango breeding, fixed SNPs associated with FBC and TC provide particularly strong gains in predictive ability and should be prioritized for marker-assisted prediction pipelines.

While our findings demonstrate that integrating GWAS-preselected variants with fixed-effect SNPs can enhance genomic predictive ability, several limitations warrant discussion. First, the relatively modest training population size used in our study may limit statistical power to detect small-effect loci and increase the risk of overfitting, raising concerns about the external validity of this approach. Additionally, the specific population structure of our study may not fully represent the broader genetic landscape of mango germplasm, potentially affecting the transferability of our findings to more diverse populations. If between-subpopulation genetic variance differs across populations, the benefits of marker preselection and fixed-effect SNP integration may not be universally applicable. Future studies should validate these results in larger, independent datasets and assess the approach’s robustness across different genetic backgrounds to ensure broader applicability.

Several inconsistencies in predictive ability across varying densities of GWAS-preselected SNPs and different GWAS models highlight the practical challenges of selecting an appropriate GWAS method for variant preselection and determining the optimal number of SNPs to include. Such inconsistencies have important downstream implications, as the choice of GWAS method and preselected variants directly influences the construction of the GRM and the inclusion of fixed-effect SNPs in prediction models, ultimately affecting prediction accuracy. To address these inconsistencies and leverage the complementary strengths of individual GWAS methods, an ensemble-based approach that aggregates summary statistics from multiple GWAS models may offer a more robust solution. Such an approach could combine p-values, effect sizes, or marker rankings to prioritize SNPs that are consistently identified across methods, thereby balancing both sensitivity and specificity. Although ensemble GWAS has primarily been applied to the identification of causative variants (Zhou et al., 2023), its potential for SNP preselection in genomic prediction remains untapped. Meanwhile, ensemble genomic prediction models which aggregate predictions from multiple methods, have demonstrated improved accuracy in maize (Tomura et al., 2025), common bean (Chiaravallotti et al., 2025), and across cattle, wheat, and human datasets (Gu et al., 2024), underscoring the potential of model integration at various stages of the genomic prediction pipeline. While ensemble GWAS remains underexplored, a practical strategy for breeders is to prioritize markers consistently identified across multiple GWAS methods and benchmark the resulting models through cross-validation. This ensures that selected SNPs are both reproducible and practically useful in applied breeding programs.

Our findings demonstrate that the predictive ability of models based on GWAS-preselected variants varies depending on the GWAS methodology employed. The superior performance of BLINK and FarmCPU compared to MLMM and the GLMM indicates their greater power in ranking markers based on LD with QTLs, thereby enabling the selection of more informative SNPs for genomic prediction. Beyond detecting a higher number of trait-associated SNPs than the MLMM and GLMM, these methods likely provide a more refined prioritization of markers with strong trait relevance. This superior performance can be attributed to their ability to effectively eliminate confounding effects between testing markers and both population structure (Q) and kinship (K) by dividing the multi-locus linear mixed model (MLMM) into components using either a fixed-effects model (FEM) and a random effects model (REM, pseudo-QTNs) in FarmCPU, or a fixed-effects model (FEM, for selecting pseudo-QTNs) and Bayesian Information Criterion (BIC) in BLINK (Huang et al., 2019; Liu et al., 2016). The use of pseudo-QTNs selected using REM in FarmCPU and FEM in BLINK as covariates effectively control false positives while retaining power to detect true associations. These features likely increase the probability of detecting SNPs that surpass the Bonferroni threshold as well as prioritizing biologically informative variants for use in genomic prediction.

The observation that, in some cases, differences in predictive ability across GWAS methods and varying densities of preselected SNPs were minimal suggests possible redundancy among SNP sets, shared association signals across GWAS methods, or the inherently polygenic architecture of the traits. One possible explanation is that methods such as BLINK and FarmCPU initially fit a general linear model (GLM), and when no significant associations are detected, they may default to reporting GLM results (Zhiwu Zhang, personal communication). This can result in overlapping sets of preselected SNPs across methods, which may explain the similar or comparable predictive abilities observed among BLINK, FarmCPU, and the GLMM for fruit firmness and trunk circumference under parental validation. A second contributing factor to the minor differences in predictive ability may be the presence of shared association signals across GWAS methods, where overlapping SNPs are selected due to consistently low p-values, suggesting potential relevance to the trait despite not reaching strict statistical significance. A third contributing factor is marker redundancy, which may occur even when the sets of GWAS preselected variants differ, if the SNPs are in LD and tag the same underlying QTLs. As a result, different sets of preselected SNPs may contribute similar genetic information to the prediction model, resulting in minimal variation in predictive ability. These modest differences are also consistent with the polygenic architecture of fruit quality traits and tree growth, where predictive ability is distributed across many loci rather than being driven by a few large-effect variants (Dong et al., 2024; Srivastav et al., 2023).

Our analysis revealed a marked decline in predictive ability when population structure was accounted for in prediction models ( Figures 2 , 3 ), a pattern consistent with that reported by Guo et al. (2014) for wheat and rice. Our findings indicate that, for these traits in the gene-pool population, a considerable portion of predictive ability is derived from across sub-population genetic variance (i.e. the model’s ability to classify individuals into their respective sub-populations), rather than solely from within sub-population genetic variance (i.e. predictive ability attributable to LD between markers and QTLs). This result is consistent with the observations of Daetwyler et al. (2012), who reported a decline in GEBV accuracy when population structure was accounted for and argued that the reduced accuracy reflects the predictive power attributable to LD between markers and QTLs.

The relatively larger gains in predictive ability with GWAS-preselected variants when population structure was accounted for, compared to models without control for population structure, likely reflect the greater contribution of LD information once the confounding effects of population structure are minimized. A previous study in the Australian mango breeding population found that TC, FBC and fruit blush intensity are strongly associated to population structure (Wilkinson et al., 2022). To avoid spurious associations, separating trait-associated loci from loci associated to ancestry is particularly important in this population. Because population structure was already accounted for during GWAS (through inclusion of PCs as fixed effects), the preselected variants are more likely to tag causative QTLs or be in meaningful LD with them, rather than merely reflecting population stratification. In contrast, WGS data contain many markers that may not be in LD with causative loci but can still contribute to predictive ability by capturing population structure. When population structure is explicitly controlled for in the prediction model, these markers provide little useful genetic signal and may introduce noise, leading to a sharper decline in predictive ability compared to models using trait-informative GWAS-preselected variants.

While our findings demonstrate a marked decline in predictive ability after accounting for population structure using fixed PCs, this sharp reduction may reflect over-correction for population structure arising from double-counting population structure effects (Hong et al., 2025). As argued by Janss et al. (2012), incorporating fixed PCs derived from the same GRM used in the random component of the model can redundantly adjust for population structure, thereby diminishing predictive ability by removing genuine genetic signals alongside confounding effects. Future studies should evaluate methods that address this issue, such as the reparameterized GBLUP model of Janss et al. (2012), which enables natural partitioning of across-subpopulation genetic variance due to population structure and within-subpopulation genetic variance that is of primary interest to breeders. Hong et al. (2025) advocated for accounting for population structure using PCs as random effects to avoid the over-correction that may occur when PCs are fitted as fixed effects in GBLUP models. However, in our study, fitting PCs as fixed effects provided conservative estimates of predictive ability, which are likely more transferable to homogeneous breeding populations or across-population predictions.