MBL population was developed from crosses between four founder lines SC630 (PI533937), SC605 (PI534096), BTx642 (PI656029), and BTxARG-1 (PI561072). These founders were selected to capture genetic diversity across multiple botanical races (SC630 = kafir, SC605 = guinea, BTx642 = durra, and BTxARG-1 = caudatum) as well as broad phenotypic variation across various qualitative, agronomic and yield related traits (Kumar et al., 2023). In addition, the MBL was genetically characterized for some well-known, heritable traits such as seed color, plant color, and awns using QTL mapping and GWAS. Full details on the development and initial characterization of the MBL panel are described by Kumar et al. (2023).

As previously mentioned, a set of 708 F₆ RILs and the four founder lines were phenotyped for PT, DTA, PH, GY, TGW, TN, and YPP. The first and last plants in each row were not phenotyped since they served to eliminate confounding results caused by border effects. Panicle type (PT) was visually assessed (C = compact, SC = semi-compact, SO = semi-open, and O = open) at physiological maturity at a single field location (Simpson Research Farm, Pendleton, SC) and used in the analysis as numerical values in 1–4 scales, (1 = compact, 2 = semi-compact, 3 = semi-open, and 4 = open). DTA was measured as days after planting to when 50% of the plants in the plot were at mid-bloom. PH was measured at physiological maturity in centimeters from the ground to the apex of the primary panicle. The tiller number per plot (TN) was estimated by counting the total number of plants from a representative 1-m row. All panicles were harvested from this 1-m section and subsequently threshed to process grain yield (GY) and individual yield component traits. Harvested panicles were dried for 3–4 days in an electric dryer to a constant moisture content (∼12% moisture) and threshed individually with a BT-14E belt thresher (Almaco, Nevada, IA, United States). GY was estimated as a total grain weight of 1-m harvested plot, which was also used for estimating TGW by counting 1,000 grains of every individual using Model U electric seed counters (International Marketing and Design Co., San Antonio, TX, United States). To estimate the grain yield per panicle (YPP), GY was divided by number of panicles. All seven traits were divided into three major categories: 1) panicle morphology (PT), 2) agronomic (DTA and PH) and 3) yield related traits (GY, TGW, TN, and YPP).

MBL (708 F₆) RILs along with their four founders were grown at two field locations including Simpson Research Farm, Pendleton, SC near to Clemson University (designated as CU) and Pee Dee Research and Education Center, Florence, SC (designated as FL) during summer 2021. The first location was planted on 2 June 2021 at Simpson Research Farm (34.624954, −82.726496) in Pendleton, SC included 701 RILs, while the second location was planted on 7 June 2021 at the Pee Dee Research and Education Center (34.287834, −79.744063) in Florence, SC, which included 708 RILs. An alpha lattice field design was used across locations to evaluate the MBL RILs and founder lines, with two replications in each location and four incomplete blocks per replicate. Each incomplete block contained 182 RILs, four founders, and a F₁ hybrid check. Each genotype was grown as a single-row plot at 3 m in length and a row spacing of 0.76 m. A plant density of ∼130,000 plants ha⁻¹ was calculated based on a plant establishment of 75% using a seeding rate (El Naim et al., 2012).

Before planting the field trials, the seeds were treated with a blend of fluxofenim (Concep, herbicide antidote), clothianidin (Nipsit, insecticide), mefenoxam (Apron XL, fungicide), and fludioxonil (Maxim XL, fungicide). Pre-plant N-P-K was applied at a variable rate based on point soil samples and worked into the soil using conventional tillage. To prevent the germination of weeds, fields were sprayed just after planting with a pre-emergent herbicide containing atrazine and S-metolachlor. A post-emergent application of atrazine was administered approximately 40 days after planting. Sugarcane aphids (Melanaphis sacchari) were controlled with one or more applications of flupyradifurone (Sivanto Prime), and chlorantraniliprole (Prevathon) was administered in a single application to prevent corn earworm (Helicoverpa zea) and fall armyworm (Spodoptera frugiperda) infestation. The field trials were irrigated as required to prevent puzzling effects on genotypic performance due to maturity and varying degrees of drought tolerance.

The MBL population was previously sequenced as described in Kumar et al. (2023). Briefly, DNA was extracted by Intertek (Alnarp, Sweden) from desiccated leaf punches collected from individual RILs and four parents at Florence, SC field site, with most plants at the grain filling stage. DNA samples were sent to Diversity Arrays Technology Pty Ltd. (Canberra, Australia (https://www.diversityarrays.com/); for Diversity Arrays Technology sequencing at low density (DArTseqLD). For DArTseqLD analysis, DNA wasdouble digested using PstI and MseI (Kilian et al., 2012), and amplified fragments were bulked and sequenced by the Hiseq2500 (Illumina^® Inc., San Diego, CA, United States). Resulting FASTQ data were processed using fastp (Chen et al., 2018) to remove barcodes and low-quality sequences before aligning reads with BWA (Li and Durbin, 2009). Variants were called using Genome Analysis Toolkit (GATK) (McKenna et al., 2010) best practices (DePristo et al., 2011; Van der Auwera et al., 2013). In brief, due to the nature of restriction digest, duplicates were not marked, but instead aligned reads went straight to base-quality recalibration. For recalibration, whole-genome sequencing data (Boatwright et al., 2022) with high quality (30x coverage) were used from the Sorghum Association Panel (Casa et al., 2008). The recalibrated BAM were then subjected to individual-sample variant calling to generate gVCFs before consolidating all gVCFs into a database for joint variant calling (Poplin et al., 2018). SNPs were hard filtered for quality (QD < 2.0, InbreedingCoeff < 0.0, QUAL < 30.0, SOR > 3.0, FS > 60.0, MQ < 40.0, MQRankSum < −12.5, and ReadPosRankSum < −8.0), missing data (50%), and minor allele frequency (>0.05) using both GATK and BCFtools, prior to GWAS. Beagle was used to impute missing genotype data in the VCF file assembled from GATK.

Pearson’s correlation coefficient matrix was generated using metan package and the corr_plot and plot.corr_coef functions were used for the visualization of correlation matrices for each trait in R software (R Foundation, 2020). Best linear unbiased predictions (BLUPs) values for each trait were calculated using random effect of the genotypes following lme4 package in R (Bates et al., 2015). BLUPs values of each trait were used as response variables in GWAS analyses. The variance components for genotypes (i.e., RILs) were estimated using the lme4 package in R (Bates et al., 2015). All effects were treated as random. The lme4() function within this package optimized the linear mixed model using restricted maximum likelihood and was implemented to determine variance components for each random effect. Because there were two locations and 1 year, replicates and locations were used in the broad-sense heritability calculation in place of year along with interaction between genotype and location to estimate the variance caused by genotype x environment interaction as shown below. H 2 = σ G 2 σ G 2 + σ G × R 2 R + σ G × L 2 L + σ E 2 R L where G is genotype, R is replicate, L is location, and E is error.

GWAS were performed using a Memory-efficient, Visualization-enhanced, and Parallel-accelerated (rMVP) GWAS program (Yin et al., 2021) installed in the R programming language (R Foundation, 2020). The rMVP package was designed to process more efficiently the large GWAS datasets, quickly evaluate population structure, and implement parallel-accelerated association tests to dramatically improve computation time. Further, rMVP provides access to several of the most popular models including the mixed linear model (MLM; Zhang et al., 2010), and fixed and random model circulating probability unification (FarmCPU; Liu et al., 2016) model, which were both used for this study. The use of MLM permits a single-locus analysis, where individuals are included as random effects and the degree of correlation among individuals is determined using a kinship (K) matrix. The use of the MLM further provides shrinkage to the model such that potential false positives due to shared ancestry are no longer significant. An MLM can be described using Henderson’s matrix notation as follows: Y = X β + Zu + e , (1) where Y is the vector of observed phenotypes; β is an unknown vector containing fixed effects, including the genetic marker, population structure (Q), and the intercept; u is an unknown vector of random additive genetic effects for individuals/lines; X and Z are the known design matrices for fixed and random effects, respectively; and e is the unobserved vector of residuals. The u and e vectors are assumed to be normally distributed with zero mean and unit variance.

FarmCPU represents a multi-locus model that iteratively uses fixed and random effect models to generate sets of pseudo-quantitative trait nucleotides (QTNs) to use as covariates and control for false positives during analysis (Liu et al., 2016). FarmCPU provides benefits over traditional MLM as it performs a multi-locus analysis, may be efficiently computed, and removes confounding between kinship and the testing marker. By iterating a fixed effect model to identify significant pseudo-QTNs to use as covariates in a random effect model using a restricted kinship matrix like the SUPER algorithm (Wang et al., 2014a) to further refine the set of included covariates by maximizing the likelihood of the random effects model. Iterations cease when no change occurs in the estimated set of pseudo-QTNs. The significant marker trait associations, corresponding to putative SNPs for each trait were determined using Bonferroni-corrected p-value threshold 1.3e⁻⁵. This threshold was calculated using 0.05/m, with m being the number of markers at 3,751.

The linkage disequilibrium (LD) decay was estimated using PopLDdecay (Zhang et al., 2019) program within a 10 Mb window. LD decay was plotted for individual chromosomes as well as genome-wide using the custom R scripts (Boatwright et al., 2021), where the coefficient of determination (r ² ) between markers located on each chromosome was measured to estimate the LD relationship between loci. The r ² was plotted on the y-axis and physical distances (Mb) on the x-axis.

GWAS were also performed for yield and yield component traits (GY, TGW, TN, and YPP) within a subset of RILs with genetic backgrounds fixed for the three dwarfing genes (Dw1, Dw2, and Dw3). For performing this analysis, we used two haplotypes (alternative alleles of each gene) based on the closest associated SNP with each gene (Dw1, Dw2, and Dw3) and performed GWAS analysis.

Functional annotation was performed using reference genome BTx623 V3.1.1 to identify candidate genes associated with significant SNPs, which were identified through GWAS. Similarly, all the significant SNPs identified were used for validating the locations of MTAs based on the sorghum QTL atlas (Mace et al., 2019; aussorgm.org.au).

Pleiotropic effects were assessed for all seven traits (PT, DTA, PH, GY, TGW, TN, and YPP) using the R package mashr, which uses a multivariate adaptive shrinkage approach to identify significant pleiotropic effects across the traits (Urbut et al., 2019). The estimated effect sizes and standard errors for every significant SNP marker in the MLM and FarmCPU for the above traits were filtered using a local false sign rate (LFSR) < 0.1 based on a condition-by-condition analysis using mashr in R (Stephens, 2017). The LFSR represents the probability of incorrectly assigning the direction of an effect. The LFSR provides a superior measure of significance over traditional multiple-testing corrections such as Bonferroni or False Discovery Rate (Benjamini and Hochberg, 1995) due to its robust estimation process (Stephens, 2017). A control set of estimated effects and standard errors were also randomly selected from the 3,751 markers to estimate the covariance between SNPs for each phenotype. Using this control set, a correlation matrix was estimated using mashr (Urbut et al., 2019) to control for any confounding effects arising from correlated traits. For testing pleiotropy across traits, canonical and data-driven covariance matrices were used. Posterior probabilities were calculated for each SNP by fitting a mash model on all tests. Bayes factors were extracted and plotted from mash results using the CDBN genomics R package (MacQueen et al., 2020). Variants exhibiting Bayes Factors greater than 10 were considered as demonstrating significant pleiotropic effects.

The contrasting features among founders foreshadowed the wide range of phenotypic diversity of the MBL as summarized in Table 1. SC605 was an early flowering parent and reached anthesis at 57 days after planting, (DAP) while BTxARG-1 flowered significantly (16 days) later at 73 DAP. The RILs of MBL population showed a wide range of variation for flowering time. Early flowering line reached anthesis at 49 DAP at CU location compared to 50 DAP at FL location. Similarly, the late flowering line reached anthesis at 91 DAP at CU location compared to 89 DAP at FL location. The mean anthesis of the MBL RILs was 62 DAP at CU instead of 65 DAP at FL location (Supplementary S1). The range of PH among the founders was relatively narrow from 101 cm (BTx642) to 119 cm (SC630). BTxARG-1 had the highest grain yield followed by SC630 (kafir). Conversely, SC605 (guinea) was a poor yielder and displayed a high tillering capacity.

In addition, the MBL population showed transgressive segregation for the majority of the traits (Table 1). The range of PH was wide in the RILs (69–201 cm) and an average of the RILs was significantly higher at 130.7 cm compared to parental lines. The results of variance component analysis demonstrated that variances due to genotype (line), genotype x location, and blocks within replication and location had significant contributions to total phenotypic variance for each trait (DTA, PH, GY, TGW, TN, and YPP) of the genotype (Supplementary Table S1). Environmental effects on phenotypic trait values were largely from blocks within replication and variance due to location for all the traits except GY, however, it showed a highly significant effect due to genotype x location. The MBL population showed wide and continuous distribution for the six quantitative traits (DTA, PH, GY, TGW, TN, and YPP) as expected (Figure 1). Estimates of broad sense heritability were lowest for TN (0.34) and highest for PH (0.84) (Table 1). Heritability was on the higher side for TGW (0.77) while moderate for GY (0.63) and YPP (0.60). Since, PT was an ordinal trait thus excluded from this analysis.

Panicle type (PT) did not show any relationship with DTA, TGW, and TN, while it showed poor positive correlations with PH (r = 23), GY (r = 0.12), and YPP (r = 0.11) (Table 2). DTA was negatively correlated with TN (r = −0.40) and TGW (r = −0.17), but no significant relationship was found with PH or GY. PH showed a significant and positive correlation with GY (r = 0.50) and yield component traits (TGW, TN, and YPP). GY had a strong positive correlation with TGW (r = 0.44) and YPP (r = 0.78) while moderate with TN (r = 0.31). TGW was positively correlated with TN and YPP. TN displayed a negative relationship with YPP. Overall, strong and positive correlations were observed among PH, GY, TGW, and YPP across the locations (Supplementary Figure S1).

GWAS were performed using MLM as a single-locus and the FarmCPU as a multi-locus model to predict genotype-by-phenotype associations. BLUPs of phenotypes were used to minimize errors across the multi-environment data more efficiently to identify QTL or MTAs in GWAS or QTL mapping studies. The MTAs were defined based on the rate of average LD decay observed in our population. In MBL, genome-wide LD fell around 2.5 Mb, therefore any two or more loci detected apart from 2.5 Mb distances were considered different MTA. Altogether, GWAS identified 70 significant associations (52 MTAs) for seven targeted traits. These 52 MTAs were distributed on 32 genomic regions (i.e., linkage blocks) across the 10 chromosomes of sorghum, except chromosome 8 (Figure 2, Supplementary Table S2, and Supplementary Figure S2). Of these 70 associations, 46 significant SNPs were identified by FarmCPU alone, only 10 were identified by MLM and 14 associations were commonly identified by both models. Several significant MTAs across the traits were corroborated with previous studies conducted for the corresponding traits while some of them were novel associations.

The FarmCPU model identified three significant MTAs for PT, which were located on Chr1 (∼24 Mb), Chr2 (∼58 Mb), and Chr6 (∼43 Mb). No significant association was identified by MLM. Of these three associations, a sole SNP (Chr6:43,244,873) overlapped with PH and GY. This genomic region at ∼43 Mb was within 0.50 Mb of the dwarfing gene Dw2 (Sobic.006G067700).

In total, 12 MTAs were significantly associated with DTA, which were located on seven different chromosomes (Chr1, Chr2, Chr3, Chr4, Chr5, Chr6, and Chr10). Among these, nine MTAs were uniquely detected by FarmCPU, whereas two MTAs were detected using MLM, and a single SNP was commonly detected by both models. Two significant MTAs between DTA and PH overlapped, which were located on Chr3 (Chr3:62,831,536) and on Chr5 (Chr5:2,588,266) (Figure 2; Supplementary Table S2, and Supplementary Figure S2).

In total, 10 significant MTAs (21 SNPs) were identified for PH, which were located on seven different chromosomes (Chr2, Chr3, Chr4, Chr6, Chr7, Chr9, and Chr10). A single MTA carrying 10 significant SNPs was associated with Chr9, which spanned a ∼1.3 Mb region (∼57 Mb) within a LD block. Out of 10 significant SNPs, one SNP was positioned on Chr9 (Chr9:57,030,394) near Dw1 (Sobic.009G229800). Another MTA was detected on Chr6 (Chr6:43,244,873; Chr6:44,567,620) in the same LD block being ∼0.43 Mb to Dw2 (Sobic.006G067700). For PH, a significant MTA was identified on Chr7 (Chr7:59,606,838), in proximity (0.21 Mb) to Dw3 (Sobic.007G163800). In addition to these PH loci, seven additional MTAs that passed the significant threshold were identified including three on Chr3, two on Chr2, and one each on Chr4 and Chr10. Interestingly, a common MTA was detected on Chr3 (Chr3:62,831,536) around ∼62 Mb for both the agronomic (DTA and PH) traits (Supplementary Table S2, and Supplementary Figure S2).

Altogether, GWAS identified 27 significant MTAs for grain yield and its components (GY, TGW, TN, and YPP) including several shared associations among the traits (Figure 2 and Supplementary Table S2, and Supplementary Figure S2). Four MTAs (five SNPs) were associated with GY, which were located on Chr6 (∼43 Mb), Chr7 (∼59 Mb), and Chr9 (∼57 Mb). All three significant loci overlapped with three major dwarfing genes (Dw1, Dw2, and Dw3). Of these five SNPs, three were independently detected by FarmCPU and a sole SNP (Chr7:59,591,981) was detected by MLM, whereas only one SNP (Chr7:59,606,838) detected by both models (FarmCPU and MLM) within the same LD block. All four MTAs showed shared associations with all traits (PT, PH, GY, TGW, TN, and YPP) except DTA.

Five MTAs (eight SNPs) were identified for TGW in total, which were located on five different chromosomes (Chr1, Chr2, Chr3, Chr9, and Chr10). Of these five associations, one MTA (four SNPs) was significantly associated with a locus on Chr9 at ∼57 Mb that spanned a genomic region (0.75 Mb). This genomic region showed significant association with multiple traits including TGW (PH, GY, and TN). In total, GWAS identified 12 SNPs (10 MTAs) significantly associated with TN, which were located on seven different chromosomes (Chr1, Chr3, Chr4, Chr5, Chr7, Chr9, and Chr10). All the significant associations were detected by FarmCPU. Three genomic regions, each located on Chr5 (∼25 Mb), Chr7 (∼59 Mb), and Chr9 (∼57 Mb) showed common associations with multiple traits (TN, DTA, PH, GY, and TGW).

Application of FarmCPU to YPP resulted in the identification of eight MTAs, which were located on six chromosomes (Chr2, Chr3, Chr4, Chr5, Chr7, and Chr9). Of these eight MTAs, a locus on Chr2 (∼59 Mb) showed a significant association with PH. Similarly, another significant MTA on Chr7 (∼57 Mb) showed association with multiple traits including YPP, PH, and GY (Figure 2 and Supplementary Table S2, and Supplementary Figure S2).

All the significant SNPs identified through GWAS following two models (FarmCPU and MLM) were used to assess pleiotropic effects for all the seven traits using the mashr program. A total of 15 significant MTAs (29 SNPs) were detected with pleiotropic effects on multiple traits (Table 3; Supplementary Table S3, and Supplementary Figure S3). These MTAs were distributed on all the chromosomes of sorghum except Chr2 and Chr8. In this analysis, PH showed nine shared associations with six additional traits (PT, DTA, GY, TGW, TN, and YPP) involving five genomic regions located on five different chromosomes (Chr2, Chr3, Chr6, Chr7, and Chr9). GY showed four shared associations with five yield components (PT, PH, TGW, TN, and YPP) that involved three genomic regions located on Chr6, Chr7, and Chr9. Similarly, TGW also showed four shared associations with three traits (PH, GY, and TN) involving a repeatedly detected locus on Chr9 (∼57 Mb). TN showed three shared associations with four traits (DTA, PH, GY, and TGW) involving three chromosomes (Chr5, Chr7, and Chr9). YPP showed two shared associations with two traits (PH, and GY) involving two genomic regions located on Chr2 (∼59 Mb) and Chr7 (∼59 Mb) (Table 3; Supplementary Table S3). Overall, the most common genomic regions located on Chr7 (∼59 Mb) and Chr9 (∼57 Mb) were significantly associated with multiple traits (PH, GY, TGW, and YPP).

GWAS were also performed for yield related traits within a subset of RILs with genetic background fixed for the three dwarfing genes (Dw1, Dw2, and Dw3) that were segregating in the MBL population. Using this approach, two additional MTAs were identified on two different chromosomes (Chr1 and Chr7) that were all independent from Dw1, Dw2, and Dw3 (Supplementary Table S4). Both these MTAs were novel genomic loci, each was located on Chr1 (∼11 Mb) and Chr7 (∼7.5 Mb). A total of eight haplotypes were formed based on the three major dwarfing genes (Dw1, Dw2, and Dw3) identified in the four founders of this population (Table 4; Supplementary Table S4).

Interestingly, we observed several colocalized or overlapped significant MTAs between four yield influencing traits (PH, GY, TGW, and YPP), which was reflected by the relationship among these traits (Table 2; Supplementary Figure S3A,B). The overall phenotypic pair-wise correlations among the four traits (PH, GY, TGW, and YPP) were significantly positive, except TN and YPP, which showed a poor and negative correlation. We identified 15 MTAs commonly associated (shared associations) with multiple traits through pleiotropic analysis particularly between positively correlated traits (PH, GY, TGW, and YPP). These correlated traits can be simultaneously improved by selecting a single trait via indirect phenotypic selection. Conversely, TN did not show any correlation with YPP and resulted in no colocalization of QTL. Similarly, PT, DTA, and PH exhibited no significant correlation between them or with the above traits. Therefore, colocalization was rarely observed between these traits (PT, DTA, and PH), which suggests that they are largely under independent genetic control in this population. QTL mapping studies made similar observations in context of the colocalization of QTL between multiple traits in sorghum (Mackay et al., 2009; Boyles et al., 2017b; Olatoye et al., 2020).

Panicle type (i.e., morphology) is an important trait that facilitated racial classification in sorghum. Panicle morphology also influences traits related to crop adaptation such as grain maturity, grain yield, and grain size (Brown et al., 2006; Hmon et al., 2013). However, the genetic architecture of panicle morphology is not completely understood in sorghum, only a few genes have been characterized so far (Morris et al., 2013; Hmon et al., 2014; Wang et al., 2021). The four founders of the MBL belong to four different races and each founder has different panicle morphology. SC630 (kafir) has a compact panicle in contrast to the open panicle type founder (SC605; guinea), whereas BTx642 (durra) and BTxARG-1 (caudatum) are semi-open and semi-compact types, respectively. We identified three significant MTAs for PT, located on three different chromosomes (Chr1, Chr2, and Chr6). A candidate gene (Sobic.006G067700) was detected within the 0.50 Mb region of the dwarfing gene (Dw2), which has been mapped earlier in sorghum (Klein et al., 2008; Hilley et al., 2017). The same locus overlapped with panicle length (Zhang et al., 2015). Other genetic loci associated with panicle length have been previously reported in sorghum on Chr2 (Shehzad and Okuno, 2015), and Chr6 (Zhou et al., 2019) in the overlapping regions. The association found on Chr1 is novel.

The MBL population demonstrated wide continuous distributions for DTA and PH, which indicates the quantitative nature of these traits (Figure 1). GWAS identified the most significant associations (12 MTAs) for DTA, which were located on seven different chromosomes. This is not surprising, because the genetic loci associated with DTA have been reported earlier on all the ten chromosomes of sorghum (Boatwright et al., 2022). Of the 12 significant associations of DTA, a locus (QDTFL1.53) on Chr1 (∼57 Mb) overlapped in previous studies (Mace et al., 2013; Burks et al., 2015). Similarly, two additional loci, one on Chr4 (QDTFL4.18) and another on Chr6 (QDTFL6.56) overlapped in earlier studies by Mace et al. (2013), and Sangma (2013), respectively. The rest of the genetic loci identified for DTA were novel.

As we know, three major dwarfing genes have been previously reported in sorghum on different chromosomes, such as Dw1 (Sobic.009G229800) on Chr9 (Brown et al., 2008; Klein et al., 2008), Dw2 (Sobic.006G067700) on Chr6 (Wang et al., 2014b; Higgins et al., 2014; Burrell et al., 2015), and Dw3 (Sobic.007G163800) on Chr7 (Morris et al., 2013; Girma et al., 2019; Boatwright et al., 2022). For PH, 10 MTAs were identified, which were located on seven different chromosomes (Chr2, Chr3, Chr4, Chr6, Chr7, Chr9, and Chr10). In this study, the most significant SNPs were located on three chromosomes (Chr6, Chr7, and Chr9), which overlapped with major dwarfing genes/QTLs known in sorghum such as Chr6 (Dw2), Chr7 (Dw3), and Chr9 (Dw1). In addition, three common genetic loci associated with PH on the three chromosomes (Chr1, Chr7, and Chr9), have been mapped earlier by Boyles et al. (2017b) using two biparental populations sharing founders (BTx642 and BTxARG-1) of the MBL population. Additional genetic loci identified in this study were also reported earlier studies in the overlapping regions, each located on the Chr3 (QHGHT3.3) (Hart et al., 2001; Feltus et al., 2006; Wang et al., 2014b), Chr4 (Wang et al., 2014b), Chr9 (QHGHT9.30) (Felderhoff et al., 2012; Takai et al., 2012; Wang et al., 2014b), and Chr10 (Parra-Londono et al., 2018) using biparental population and diverse panel of sorghum.

GY is a complex trait determined by many yield components like grain weight, tiller number per unit area, and grain yield per panicle (Boyles et al., 2016; Boyles et al., 2017b; reviewed in Baye et al., 2022). Since GY and its component traits are polygenic in nature, they exhibit continuous phenotypic distributions (Figure 1). For GY, three of the four MTAs detected each on Chr6 (∼43 Mb), Chr7 (∼59 Mb), and Chr9 (∼57 Mb) were overlapped with three major dwarfing genes namely Dw2, Dw3, and Dw1, respectively. These observations indicated the potential influence of dwarfing genes on other traits. Genetic loci associated with GY have been overlapped on Chr6 (Leiser et al., 2014), Chr7 (Guindo et al., 2019), and Chr9 in sorghum (Sabadin et al., 2012; Reddy et al., 2013; Reddy et al., 2014; Boyles et al., 2017b). A novel significant MTA was detected on Chr9 (∼52 Mb). Interestingly, these three genomic regions exhibited significant shared associations with all the phenotypic traits attempted in this study, except DTA.

Grain weight is one of the major yield components and highly heritable traits (Boyles et al., 2016; 2017b; reviewed in Baye et al., 2022). In total, five significant MTAs were identified for TGW, which were located on five different chromosomes (Chr1, Chr2, Chr3, Chr9, and Chr10). Of these five, three MTAs overlapped on Chr1 and Chr2 (Boyles et al., 2017b; Patil et al., 2019), and Chr3 (Sabadin et al., 2012). Another locus for TGW has been mapped in the same region on Chr9 (∼57 Mb) in previous studies (Reddy et al., 2013; Boyles et al., 2016). The above locus also showed a shared association with multiple traits (PH, GY, and TN). However, a locus associated with TGW on Chr 10 (∼35 Mb) was identified for the first time in our study. In continuation of grain yield components, TN is also an important trait, that determines the grain yield by increasing the number of tillers or panicles per unit area in sorghum. For TN, 10 significant MTAs were identified, which were located on seven different chromosomes (Supplementary Table S2). Of these 10 associations, three MTAs were previously mapped using multiple populations in the overlapping regions on Chr1, and Chr3 (Alam et al., 2014; Kong et al., 2014), and Chr4 (Alam et al., 2014). A genomic region associated with multiple traits (TN, GY, and YPP) was identified on Chr9 (∼57 Mb), which was reported earlier for tiller numbers (Feltus et al., 2006; Zhang et al., 2015). However, another locus associated with the same set of multiple traits (TN, GY, and YPP) on Chr7 (∼59 Mb) was identified in this study, which is a novel association for tiller numbers. YPP directly influences the overall grain yield in sorghum by several factors such as high grain weight and grain number per panicle. Altogether, eight significant MTAs were identified for YPP, which were located on six chromosomes (Chr2, Chr3, Chr4, Chr5, Chr7, and Chr9). A locus was identified on Chr2 (∼59 Mb) that was previously mapped in the overlapping region (Shehzad and Okuno, 2015). Two genetic loci each on Chr7 (∼59 Mb) and Chr9 (∼57 Mb) were identified for YPP. Both the common genomic regions associated with grain yield components located on Chr7 and Chr9 have been previously mapped in the same region (Boyles et al., 2017b) using two biparental RIL populations sharing founders (BTx642 and BTxARG-1) of the MBL population.

High resolution mapping is a critical step for the identification of novel genes and narrows down the genetic distance between the candidate genes associated with complex traits in crop plants (Huang et al., 2015; Scott et al., 2020; Kumar et al., 2023). As we reported earlier (Kumar et al., 2023), the genome-wide LD decayed much more quickly in the MBL compared to a biparental population (Supplementary Figure S4), which is consistent with previous comparisons between biparental and MBL populations (Boyles et al., 2017a). Based on the previous mapping studies, the MAGIC populations are more efficient to narrow down the genomic regions compared to biparental populations (Huang et al., 2012; Kumar et al., 2023). Here, we compared our results with an earlier QTL mapping study performed by Boyles et al. (2017b), where two different biparental RIL mapping populations were used for QTL discoveries sharing two of the four founders (BTxARG1 and BTx642) of the MBL. Comparing the common genomic regions identified on three chromosomes (Chr2, Chr7, and Chr9) in both the studies for the same set of quantitative traits (PH, TGW, and YPP), MBL placed all the common QTL in narrow genomic regions as compared to biparental populations (Boyles et al., 2017b). The potential of the MAGIC populations in facilitating high resolution mapping have been reported earlier in crops with even more complex genomes like wheat (Huang et al., 2012; Stadlmeier et al., 2018) and cotton (Islam et al., 2016).

Despite MBL parents all having short stature (Table 4), three of the four founders displayed dominance at a unique locus SC630 (Dw1), BTxARG1 (Dw2), and SC605 (Dw3). As a result, all eight gene combinations were present across MBL progeny, which led to significant height variation. Because plant height confounds grain yield and related traits, these dwarfing genes tended to dominate explained phenotypic variance and thus masked the effects of other genetic variants segregating in the population. To identify novel alleles associated with yield related traits, GWAS were performed using a subset of RILs that contained matching plant height gene combinations. Using this approach, seven significant genomic loci were identified on five different chromosomes (Chr1, Chr3, Chr4, Chr6, and Chr7) for yield influencing traits. These genomic regions overlapped with multiple yield component traits (GY, TGW, TN, and YPP). Of these seven, the two significant loci (Chr1 and Chr7) identified appear to be novel allelic combinations. This MBL population was designed to identify and recombine useful genetic variation present in diverse maintainer (B-line) germplasm. Desirable allelic combinations were indeed elucidated based on phenotypic comparisons of all the MBL lines for yield related traits. Considering plant height allelic combinations, five best-performing MBL progeny lines were selected (Table 4). These favorable recombinants should be hybridized with elite pollinator parents to test for general combining abilities to help determine their value in sorghum hybrid development. To reinforce favorable recombination, various MBL progeny exhibited transgressive segregation for all the phenotypic traits (DTA, PH, GY, TGW, TN, and YPP), which indicates the distribution of positive and negative alleles in the founders. Such segregation and allele shuffling provide opportunities for trait improvement through pyramiding desirable alleles from selected progeny lines.