QTL analysis in multiple sorghum mapping populations facilitates dissection of the genetic control of agronomic and yield-related traits in sorghum [Sorghum bicolor (Moench)]

The genetic architectures of agronomic and yield-related traits are expected to involve multiple loci that are unlikely all to segregate for alternative alleles in a single biparental population. Therefore, the identification of quantitative trait loci (QTL) that are expressed in diverse genetic backgrounds of multiple bi-parental populations provides evidence about both background-specific and common genetic variants. The purpose of this study was to map QTLs for agronomic and yield related traits using three connected mapping populations of different genetic backgrounds, to gain insight into the genomic landscape of these important traits in elite Ethiopian sorghum germplasm. The three bi-parental populations, each with 207 F 2:3 lines were evaluated using an alpha lattice design with two replications under two moisture stress environments. Data analysis was done separately for each population using composite interval mapping, finding a total of 105 QTLs. All the QTLs identified from individual populations were projected on a combined consensus map, comprising a total of 25 meta QTLs for seven traits. The consensus map allowed us to deduce locations of a larger number of markers than possible in any individual map, providing a reference for genetic studies in different genetic backgrounds. The meta QTLs identified in this study could be used for marker-assisted breeding programs in sorghum after validation. Only one trait, reduced leaf senescence, showed a striking bias of allele distribution, indicating substantial standing variation among the lines that might be employed in improving drought tolerance of sorghum.

most important cereal crops used for food and feed in different parts of the world (Biswas et al. 2001). The crop is well adapted to hot dry environments and regarded as a model for studying drought resistance among the grasses. It has broad ecological adaptation where it is commonly grown in diverse climates, from the dry lowland areas of the semi-arid tropics receiving very little annual precipitation, to high rainfall and humid regions, as well as in the more temperate regions of the world.
It is well-known that genetically controlled variation in agronomic and other characters influence grain yield, the most important agronomic trait of sorghum (Sanchez et al. 2002). Understanding the genetic control of those traits and applying the knowledge in sorghum breeding programs might be instrumental to develop improved germplasm. Important agronomic traits in sorghum, as in most crop species, are quantitatively inherited, which complicates genetic analysis. Classical genetic studies have contributed information regarding the manner of genetic control of some of these characters (Haussmann et al. 2002). The discovery of new quantitative trait loci (QTLs) for those traits is an effective method to support sorghum improvement.
Partitioning of grain yield into its component traits could help to increase the efficiency of genetic gain through selection and understand the genomic basis of local adaptation in sorghum . Moreover, secondary traits including phenology and agronomic related traits have been recognized as important constituents of grain yield in sorghum (Kapanigowda 2011;Blum 2005;Borrell and Hammer 2000) and the use of such secondary traits can improve selection efficiency (Araus et al. 2002). Plant height, flowering time, head exsertion and many other traits were reported to have association with grain yield in sorghum Feltus et al. 2006). The present study also makes use of these traits. Genetic manipulation of these traits can positively affect grain yield and could be very helpful to breeders (Yuan et al. 2008).
Understanding the genetic bases of agronomic, phenological and yield related traits may accelerate sustained improvement of sorghum. Mapping QTLs that control variation in traits of agronomic importance is a key part of marker-assisted breeding. Quantitative trait locus (QTL) mapping has been useful in identifying and localizing important genomic regions controlling quantitative traits in a wide range of species (Tanksley et al. 1989). It has become a routine approach for genetic studies of complex traits in plants (Li et al. 2003). Efforts to develop large numbers of molecular markers and high density genetic maps using bi-parental populations have become routine for many crop species, permitting one to simultaneously define gene action and breeding value at hundreds and often thousands of loci distributed across genomes. The results from such mapping studies provide greatly improved estimates of the number of loci, allelic effects and gene action controlling different traits (Yu et al. 2008;Holland 2007), valuable information for investigation of the genetic control of complex traits (Li et al. 2011).
Genetic mapping studies using single bi-parental populations are effective in detecting QTLs but suffer from low resolution and limited allelic diversity. Because only two parents are involved in population development, QTL from single bi-parental population can only be assumed to be relevant to the cross and environment(s) in which they were mapped. Agronomic and yield related traits are very complex traits and any single bi-parental population is unlikely to segregate for all genetic loci influencing it. The identification of QTLs in different genetic backgrounds provides evidence whether a common genetic variant(s) determine those traits in different sorghum genotypes. Such common genetic variants (Meta QTLs) have broad based expression and are potential candidates for marker-assisted transfer of traits into locally adapted materials. Populationspecific QTLs can also be exploited through markerassisted pyramiding to accumulate alleles conferring wider adaptation.
Several linkage mapping studies have been conducted in sorghum to dissect the genetic mechanisms controlling different traits using individual bi-parental mapping populations (Shehzad et al. 2013;Bibi et al. 2012;Gebisa et al. 2007;Haussmann et al. 2002;Sanchez et al. 2002;Kebede et al. 2001;Xu et al. 2000). As a result, genomic regions (QTLs) associated with agronomic, phenologicaland grain yield related traits were identified and reported (Mace et al. 2019). For example, variation in maturity has been suggested to be controlled principally by six genes (Rooney and Aydin 1999;Quinby 1966) and dwarfism by four genes (Quinby and Karper 1945) although QTL mapping studies have found many Page 3 of 22 24 Vol.: (0123456789) more genomic regions to influence these traits (Zhang et al. , 2013. Traits such as plant height, flowering time, leaf senescence and panicle weight have also been characterized using different segregating populations (Nagaraja et al. 2013;Srinivas et al. 2009;Brown et al. 2006;Hart et al. 2001). However, identification of QTL that are expressed in more than one mapping population of diverse genetic background is lacking. Keeping the above points in view, the objective of this study was to dissect the genetic architecture of grain yield and agronomic traits using linkage analysis and identify potential mechanisms of complex traits in sorghum using Ethiopian elite genetic materials.

Descriptions of experimental sites
The present study was conducted during the 2016/2017 cropping season at Kobo (Northern Ethiopia) and Meiso (Eastern Ethiopia) research stations. Meiso is located at 14° 23′ N, 37° 46′ E, at an altitude of 1394 masl and Kobo is located at 12° 09′ N, 39° 38′ E, at an altitude of 1468 masl. The study sites represent moisture stress and major sorghum producing areas in Ethiopia (MoA 1998). The monthly weather conditions of the two sites during the experimental season are presented in Fig. 1. Although there was a slight difference between the rainfall received at Kobo and Mieso during the crop growth period, in both locations over 35% of the total rainfall was received within two months of the cropping season ( Fig. 1) suggesting that the crop experienced considerable drought stress in most of its growth stages. According to the FAO soil classification system, the major soils of Meiso are vertic cambisol and haplic luvisol. The Kobo site is mostly occupied by Eutricvertisol/ Eutricfluvisol soils, and the dominant soil texture is clay (FAO 2007).

Genetic materials
Three mapping populations developed from crosses of 76T1-23 × Baji, Meco-1 × Birmash and 76T1-23 × Birmash, each with 207 F 2:3 lines and two parents were used for this study. The parents are contrasting for yield and other agronomic traits. The maternal parents Meco-1 and 76T1-23 are low-yielding, early maturing and relatively drought tolerant varieties whereas the paternal parents Baji and Birmash are high-yielding, late maturing and relatively susceptible to moisture stress. These four elite sorghum varieties belong to the caudatum race and are widely cultivated in different parts of Ethiopia. For each population, F 1 seeds were obtained by crossing two parents and the F 1 plants were selfed to produce a large number of F 2 , then F 2:3 progenies were developed using single seed descent.

Experimental design and field management
The experiment was laid out in a 10 × 21 alpha lattice design with two replications. Twenty one genotypes were assigned in each of the 10 incomplete blocks for all populations. The total size of a single block was 63 m 2 (15.75 m × 4 m). The experiment was carried out during 2016/17 cropping season. Seeds of the F 3 and their parents were sown by dropping in the furrows opened by the plough in two rows, with 75 cm spacing between rows. Seedlings were thinned three weeks after sowing with 20 cm spacing between individual plants which resulted in a plant density of 40 plants per plot. The net plot size was 0.75 m × 4 m. Inorganic fertilizers, DAP (Di-ammonium phosphate) and Urea were added at rates of 100 and 50 kg/ha as side dressing during sowing and three weeks after sowing, respectively. Chemicals were sprayed to protect the experiment from insect and disease damage. In both sites, all recommended agronomic practices were followed throughout the cropping seasons.

Phenotyping
Data were collected for seven agronomic and yield related traits in sorghum (Tables S1, S2). Five plants from the two rows were randomly chosen for recording agronomic and yield related traits. Days to flowering (DF) and days to maturity (DM) were the number of days from emergence to flowering and days to maturity respectively. Plant height (PH) was measured as the distance from the soil level to the tip of the main stem panicle. Panicle weight (PW) was measured as the average weight of dry from five randomly selected plants in grams (g). 1000-seed weight (TSW) was measured as the average weight of 1000 dried seeds in grams sampled randomly from each plot after harvesting, weighed and adjusted to a standard moisture content (12%). Grain yield per panicle (GYP): Grains harvested from five randomly selected plants of each genotype were dried, weighed and adjusted to a standard moisture content (12%) and their average was expressed in grams as grain yield per panicle. Leaf senescence (LS) was recorded by visual ratings on a scale of 1 to 5 based on the degree of leaf death at maturity on a plot basis following sorghum descriptors (IBPGR 1993).

DNA extraction
Leaf samples were collected from two weeks old five randomly selected seedlings and pooled together. Genomic DNA of the pooled tissues were extracted following the CTAB (hexadecyltrimethylammonium bromide) protocol (Maroof et al. 1984). The quality and purity of the isolated DNA was checked using agarose (0.7%) gel electrophoresis and Nano drop reading. Only high-quality DNA (of high molecular weight and non-degraded) was used for subsequent genotyping. After checking the quality, DNA was diluted with 50 ul TE and shipped for genotyping to the Plant Genome Mapping Laboratory at University of Georgia, USA.

Genotyping-by-sequencing (GBS)
Marker development GBS libraries were prepared using a PstI and MspI enzyme system (Poland et al. 2012) with modifications (Dong et al. 2018). Each DNA sample was digested with the rare cutting enzyme PstI and the frequent cutter enzyme MspI and ligated to a unique barcoded adapter and a common adapter. Samples were pooled and the 200-500 base pair (bp) size fraction was extracted from a 2% agarose gel after electrophoresis and purified using a Qiagen Gel Extraction Kit (Qiagen, Hilden, Germany). The purified DNA was PCR amplified using NEB 5X Master Mix (New England Biolabs Inc., Ipswich, MA, USA) and the PCR products were extracted as above to eliminate primer dimers. Parental DNA samples were replicated six times for sequencing so as to improve read depth and reduce missing data, which is especially important for correctly calling heterozygous loci in the progeny. All libraries were either sequenced twice on a Miseq or once on a Nextseq to increase the sequencing depth.

SNP calling
The reads obtained were first de-multiplexed according to the sample barcodes and adapter sequences were removed using a custom perl script. SNP discovery was performed using the GBS pipeline in TAS-SEL version 5.0 (Bradbury et al. 2007). Minimum quality score of 10, minimum call rate of 0.8, minor allele frequency of 0.05 and kmer length of 100 bp were considered during SNP discovery. SNPs were initially called across all three populations, obtaining a total of 16,064 SNPs. Then the raw SNP dataset was separated into three datasets for population 76T1-23 × Baji, Meco-1 × Birmash and 76T1-23 × Birmash, respectively (Tables S4, S5 and S6). After calling SNPs in the Tassel GBS pipeline, SNPs were further filtered in R (version 3.2.4). The distributions of parental genomes were checked (Fig. 2) and low coverage SNPs (reads ≤ 4) were converted into missing data to mitigate heterozygote undercalling errors. SNPs with greater than 30% missing data were discarded. Rare individuals with greater than 50% missing data, likely due to poor DNA quality, were also discarded. Goodness of fit tests were conducted after filtering out SNPs with > 30% missing data. In the 1:2:1 ratio test, missing data were simply left as missing without counting for genotype classes. Markers with p < 0.05 in chi-square test were declared as segregation distorted markers. We used p = 0.001 as the cutoff to retain many statistically segregation distorted markers to complement the biological meaning of segregation distortion. Finally a total of 2643 SNPs were retained for the three populations. Physical positions of generated SNPs were obtained based on alignment against the Sorghum bicolor genome v1.4 assembly (Paterson et al. 2009).

Genetic map construction
Genotypic data generated in this study was used for constructing genetic maps for each population using R/qtl (Broman 2003). For co-segregating markers (with identical genotypes), i.e. mapping to the same location without providing additional information, one was removed. Recombination fractions and LOD scores were estimated for all pairs of markers within each population. Then, the SNP markers were assigned to linkage groups (LGs) (Table S3) at a minimum logarithm of odds (LOD) score of 6 and maximum recombination fraction of 0.35. Genetic distances between SNPs within each LG were estimated using the Kosambi mapping function (Kosambi 1944). SNP orders within each LG were compared to corresponding physical positions on the S. bicolor genome assembly (Paterson et al. 2009).

Phenotypic analysis
A multi-environment analysis of variance (ANOVA) was first conducted across the two environments for each population using R software (version 5.6) to see the effects of environment, environment by genotype interaction and genotype. As incomplete blocks contain fewer genotypes than the total number of genotypes to be compared, block adjustment was done. Given the distinct conditions across the two environments, trait best linear unbiased predictions (BLUPs) ( Table S7) were used to estimate the values for each line within each environment and used for QTL analysis using a mixed linear model (Bates et al. 2015). Analysis of phenotypic data was based on the following model: where G is genotype, R is replication, B(R) is block nested within replication, error is random error.
The random effects are estimated by BLUP (Wright 1968). Broad-sense heritability was calculated from an ANOVA fitting effect of genotype (G) and environment (E), as where Hb 2 is broad sense heritability, σ 2 G is genotypic variance, σ 2 E is error variance, r is the number of replications (Nyquist 1991).

QTL analysis
Genotypic and phenotypic data obtained in the present study were used for QTL analysis using QTL Cartographer v.2.5 (Wang et al. 2012). Given the distinct conditions across the two environments, trait BLUPs were estimated for each line within each environment and used for QTL analysis using a mixed linear model implemented in lme4 package in R software (Bates et al. 2015). QTL analysis was performed for individual environments and across environments for the three populations separately. Composite interval mapping (CIM) was performed by selecting Model 6 with the default window size 10 cM and control marker number 5. To obtain more precise results the default walk speed was reduced to 1 cM. Finally, the location of a QTL was declared as a position where the LOD score value ≥ 2.5 based on 1000 permutation test (Wang et al. 2005). When two adjacent LOD peaks fell in a common support interval, only one QTL with the highest peak was considered present, and its position was taken as the location of the QTL (Dufey et al. 2012). The proportion of phenotypic variance explained by a single QTL was obtained by the square of the partial correlation coefficient (R 2 ). Estimates of the additive effects of the QTLs were obtained by fitting a model including all putative QTLs for the respective trait (Marinoni et al. 2017). The identified QTL were designated with an italicized symbol composed of a q, a trait abbreviation, the chromosome number in which the QTL is located, a hyphen, and in cases where more than one QTL controlling a trait were detected in the same LG, they were numbered serially (Sivakumar et al. 2016). QTLs were classified as major if the phenotypic variance explained was larger than 10%, and minor when it was less than 10% (Marinoni et al. 2017).
Meta QTL analysis was performed using Biomercator software v2.1 (Danan et al. 2011). All QTLs identified in the individual populations were projected on the consensus linkage map. SNP markers overlapping on at least two genetic maps were selected as anchor markers and used to integrate corresponding linkage groups on individual linkage maps. Then all QTLs identified from the three populations were projected onto the integrated map based on the chromosomal position, LOD score and proportion of phenotypic variance (R 2 ) explained by each QTL. The lowest Akaike Information Criterion (AIC) value was used to select the best QTL model for each chromosome (Hirotugu 1974). The selected model was used to fit the distribution of the "projected QTL" onto a chromosome, and to cluster them to determine the Page 7 of 22 24 Vol.: (0123456789) number of QTLs underlying the distribution of the observed QTL using the QTLclust command (Blanc et al. 2006). Finally, all the meta QTLs (mQTL) identified were compared to the Sorghum QTL Atlas database described in Mace et al. (2019). Known sorghum genes within mQTL region were also identified using adjacent markers and searched against the sorghum genome.

Data availability
In addition to data within this manuscript, we present supplemental files via FigShare. Table S3, contains individual linkage maps of the three populations. Table S4, Table S5 and Table S6 contain the genotypic data of the populations: 76T1-23 × Baji, Meco-1 × Birmash and 76T1-23 × Birmash, respectively. The raw phenotypic data of the three populations: 76T1-23 × Baji, Meco-1 × Birmash and 76T1-23 × Birmash used in the analysis are in Table S1 and  Table S2. Table S7, contains BLUPs of the three populations. All supplemental data are available centrally at FigShare (https:// figsh are. com/ search? q= techa le). Raw sequencing data are available in the NCBI Sequence Read Archive under project accession BioP roject ID PRJNA687679. Data analysis scripts have been deposited to GitHub (https:// github. com/ hxdong-genet ics/ Ethio pian-Sorgh um-F3).

Phenotypic variation and heritability
In this study, the analysis of variance of each population showed that, all seven traits were significantly affected by genotype and genotype-by-environment interaction. In the two environments, average flowering of the F 2:3 populations occurred earliest from population 76T1-23 × Baji at Kobo (67 days) after emergency while late flowering was from population Meco-1 × Birmash at Kobo (83 days) (Table 1). Overall, population 76T1-23 × Baji flowered approximately 13 days earlier than population Meco-1 × Birmash. Similarly, early maturity was registered from population 76T1-23 × Baji, followed by 76T1-23 × Birmash and Meco-1 × Birmash. Average grain yield per panicle ranged from 43.92 g for population 76T1-23 × Birmash to 57.61 g for population The estimated heritability of traits was higher than 0.39 in all three populations (Table 1). Heritability estimates ranged from 0.39 to 0.88, with GYP presenting the lowest heritability (0.39-0.62). The heritability of days to maturity (0.55-0.75) was consistently higher than that of GYP, but lower than those of days to flowering (0.67-0.88) and 1000-seed weight (0.69-0.82) across the three populations. High broad sense heritability was also recorded for plant height (0.60-0.74).
A total of 105 putative QTLs (27, 42 and 36 from the three populations, 76T1-23 × Baji, Meco-1 × Birmash and 76T1-23 × Birmash, respectively) were identified for seven agronomic traits over the two environments. The Meco-1 × Birmash population provided higher numbers of QTLs (42) than the other two populations. The QTLs related with traits and their relative positions on each chromosome were detected in all chromosomes. However, most QTLs were concentrated on Chrs. 1, 3, 4, 5, 6, 8 and 9. Quantitative trait loci (QTLs) detected for drought tolerance related traits and their QTL information are summarized and presented in Tables 2, 3 and 4. For days to flowering (DF), a total of eleven QTLs were detected from the three populations over the two environments with two each on Chrs. 3, 7, 8, and 10 and one each on Chr. 2, 3 and 5. Two significant QTLs (qDF4-1and qDF7-1) were identified at M eiso on Chr. 4 and 7 which accounted for 16% of the total phenotypic variance. Similarly, one major QTL (qDF2-1) explaining 13% of the total phenotypic variance and four other minor QTLs were detected on Chrs. 3, 4, 8 and 9 for days to flowering from Meco-1 × Birmash population. These QTLs had additive effects that ranged from -1.09 to 1.78 and dominance effects that ranged from -2.60 to 2.65 respectively, in total explaining 65.16% of phenotypic variance. Favorable alleles for 7 QTLs with negative additive effects came from 76T1-23 (6 QTLs) and Meco-1 (1 QTL), with those for 4 QTLs with positive additive effects from Birmash. From population 76T1-23 × Birmash, five significant QTLs were identified over the two locations for days to flowering on Chrs. 5, 7, 8 and 10. These QTLs explained a cumulative phenotypic variance of 35.37%. The largest effect QTL on Chr. 7, qDF7-1 explained 13% of the total phenotypic variance.
For days to maturity (DM), eight QTL were identified on Chrs. 3, 4, 5, 6 and 9 over the two environments from the Meco-1 × Birmash population. The phenotypic variation explained by each QTL ranged between 4 and 18%, collectively accounting for 67.68% of the total phenotypic variation. These QTLs had additive effects ranging from -1.88 to 2.23 and dominance effects ranging from -2.60 to 2.58. Favorable alleles for 3 QTLs with negative additive effects came from Meco-1, with those for 5 QTLs with positive additive effects were from Birmash. The largest effect QTLs are qDM1-1 and qDM6-1, explaining 18% and 12% of the total phenotypic variance. A days to flowering (qDF9-1) and days to maturity (qDM9-2) QTL were co-located on Chr. 9 from Kobo. Relatively fewer QTLs were detected from populations 76T1-23 × Birmash and 76T1-23 × Baji than from Meco-1 × Birmash. Four QTLs were identified at Meiso and Kobo on Chrs 3, 4, 6 and 9 with LOD scores of 3.26, 2.83, 4.20 and 4.30 respectively; explaining 28.14% of the cumulative phenotypic variance from population Meco-1 × Birmash. Meanwhile, two significant QTLs on Chr. 4 and 10 at Kobo with varying magnitude of effects were detected from population 76T1-23 × Baji. These QTLs had an additive effect that ranged from 0.86 to 1.21 and dominance effects that ranged from -2.61 to -1.38. For plant height (PH), two major and three minor QTLs were detected on Chrs. 2, 3, 5, 6 and 7 from the individual environment QTL analysis in population 76T1-23 × Baji. These QTLs had additive effects that ranged from -5.38 to 15.34 and dominance effects that ranged from -7.25 to 10.13. The phenotypic variance accounted by these QTL was 55.65%. Favorable alleles for 3 QTLs with negative additive effects derived from 76T1-23, with those for 2 QTLs with positive additive effects from Baji. In addition, one region on Chr. 1 showed one major QTL (qPH1-1) strongly associated with plant height from the pooled  analysis. The favorable allele for qPH1-1 came from 76T1-23. Six QTLs were detected for plant height on Chrs. 2, 3, 5, 6 and 10 from the individual environments in population Meco-1 × Birmash, collectively explaining 49.04% of phenotypic variance. Regions on Chrs. 5 and 6 showed relatively larger effect on plant height and were consistently identified under both environments. From population 76T1-23 × Birmash, two major and three minor QTLs were detected for plant height on Chrs. 2, 4, 5, 6, and 9 at Meiso and Kobo. Two of the largest effect QTLs qPH6-1and qPH2-1were identified on Chr. 6 and 2, respectively explaining 37% the cumulative phenotypic variance. Ten QTLs were strongly associated with grain yield per panicle (GYP) on Chrs. 2, 3, 6, 7 and 8 in the Meco-1 × Birmash population, more than that of the other two populations. Two QTLs (qGYP6-1 and qGYP7-1) were consistently detected over the two environments, and seemed to be particularly important for determining grain yield per plant under moisture stress conditions. A significant QTL qGYP6-1 was detected on Chr. 6 from the two populations explaining a cumulative phenotypic variance of 17%. Five QTLs on Chrs. 2, 4 and 6 were identified for grain yield across the two environments from population 76T1-23 × Baji. A QTL (qGYP6-1) from Meiso on Chr. 6 was co-located with QTL qPW6-2 for panicle weigh from Kobo. From population 76T1-23 × Birmash, four QTLs were detected over the two environments, one (qGYP3-1) being consistently identified on Chr. 3 under the two environments explaining a cumulative phenotypic variance of 31.96%. The QTL (qGYP3-1) with the largest effect on Chr. 3, explained 24% of the total phenotypic variance. Favorable alleles for 8 QTLs with negative additive effects came from 76T1-23 (4 QTLs) and Meco-1 (4 QTLs), with those for 9 QTLs with positive additive effects from Baji (2 QTLs) and Birmash (8 QTLs).
For leaf senescence, a total of 14 QTLs were identified from the three populations over the two environments. The QTL qLS1-1 on Chr. 1 was consistently identified in the two environments on exactly the same position from population 76T1-23 × Birmash,  accounting for 19.03% of the total phenotypic variation. A significant QTL qLS2-1on Chr. 2 was also stable in the two environments and co-located with QTL qTSW2-1 for TSW. Almost all the QTLs had negative additive effects except one; indicating that alleles for lower leaf senescence for this trait was contributed by the parent 76T1-23, as expected. Four QTLs (one major and three minor) on Chrs. 6, 1 and 5 were identified from population 76T1-23 × Baji, accounted for 17.03% of the total phenotypic variation. From population Meco-1 × Birmash, five significant QTLs were identified at Meiso for LS on Chrs. 2, 3, 5 and 10 with LOD scores ranging from 2.82 to 5.2, collectively accounting for 28.28% of the total phenotypic variance. Only one QTL (qLS3-1) was identified on Chr. 3 from the combined data, which mapped close to qLS3-1 from Meiso. For panicle weight per plant (PW), four QTLs were identified on Chrs. 2, 3, 5 and 6 from individual environments in population 76T1-23 × Baji. The QTLs for panicle weight at Kobo and Mieso were different. One major QTL (qPW2-1) was identified on Chr. 2 from Kobo explaining 25% of the phenotypic variation. Three QTLs were identified on Chrs. 2, 3 and 8 for panicle weight per plant from population Meco-1 × Birmash. In addition, one major QTL was identified on Chr. 3 explaining 12% of the phenotypic variations with LOD score of 3.35 from Mieso. More QTLs were detected for panicle weight from popula-tion76T1-23 × Birmash than from the two other populations, with nine QTLs identified on Chrs. 3, 4, 5, 6, 7, 8 and 10 in the two environments. The QTLs had additive effects that ranged from -6.33 to 7.91, dominance effects that ranged from -7.56 to 7.86, with phenotypic variance of 61.83%. Favorable alleles for 5 QTLs with negative additive effects came from 76T1-23, with those for 4 QTLs with positive additive effects from Birmash. QTLs qPW3-2 on Chr. 3 at Kobo coincided with grain yield QTL qGYP3-1at Meiso.
For 1000-seed weight (TSW), a total of 12 QTLs were identified from the three populations over the two environments. The QTLs were mainly concentrated on Chrs. 2, 5, 6, 8 and 10. The QTL (qTSW8-1) on chromosome 8 was consistently detected in the two environments in population 76T1-23 × Baji. Four QTLs were detected from population Meco-1 × Birmash, of which two of the QTLs (qTSW7-1 and qTSW8-1) were consistently detected over the two environments in population 76T1-23 × Birmash in an over lapping positions. Five QTLs were detected from population 76T1-23 × Baji. These QTLs explained a cumulative phenotypic variance of 40.10%. The largest effect QTL on Chr. 8, qTSW8-1 explained 11.12% of the total phenotypic variance. The QTL had additive and dominance effect of 3.32 and -1.55, respectively. Favorable alleles for 7 QTLs with negative additive effects came from 76T1-23, with those for 5 QTLs with positive additive effects from Birmash.

Meta QTL Mapping (mQTL)
All QTLs identified from individual populations were projected on the consensus map for meta-QTL analysis (Fig. 4). The meta-analysis reduced the total number of QTLs from 105 to 25 (Table 5). Five of these mQTLs were specific to grain yield per panicle, five to days to maturity, four to days to flowering, three to plant height, four to panicle weight and the remaining four to leaf senescence. These QTLs were mainly concentrated on Chrs. 2, 3, 4, 5, 6, 7, 9 and 10. The number of mQTL identified on each chromosome varied from 2 on Chr. 7 to 5 on Chr. 3, with an average of 3.12. The mean phenotypic variance explained by each mQTL varied from 6.94 to 16.00% and the overall average was 10.90%.
For grain yield per panicle, five mQTLs were identified on Chrs. 2, 3, 5 and 6, explaining a cumulative 52.33% of phenotypic variance. The regions on these chromosomes harbored QTL shared by two populations. Five mQTLs for days to maturity and four mQTLs for days to flowering were detected on Chrs. 2, 3, 4, 7, 9 and 10. For leaf senescence, four mQTLs were identified on Chrs. 3, 5, 7 and 10, explaining a cumulative 41.57% of phenotypic variance. A total of 3 mQTLs were detected for plant height distributed on Chrs. 2, 5 and 6, explaining a cumulative 28.46% of phenotypic variance. Four mQTLs were detected for panicle weight on Chrs. 3, 5 and 6 with the lowest and highest variance of 6.94 and 13.13%, respectively. Meta QTLs were detected on almost all chromosomes for various traits, however; Chr. 5 harbored mQTL for leaf senescence, plant height, grain yield per panicle and panicle weight.

Discussion
In order to uncover genomic regions associated with agronomic traits including grain yield, we developed three interconnected bi-parental (F2:3) mapping populations by crossing four parents and evaluated them using the same phenotyping protocol. Meco-1-1 and 76T1-23 are early maturing and relatively drought tolerant varieties whereas the paternal parents Baji and Birmash are high-yielding, late maturing and relatively susceptible to moisture stress. The studied traits showed allele distribution in both directions, indicating that both parents contributed favorable alleles for the traits to some degree for example, as found by Paterson et al. (1988). Only one trait, reduced leaf senescence, showed a striking bias of allele distribution, with 13 of the 14 favorable alleles coming from the early maturing 76T1-23.
In our study, the estimated heritability of all seven traits was high (> 0.39) in all three populations. High broad sense heritability were recorded for plant height, days to flowering and 1000-seed weight across all populations, indicating that most phenotypic variation appeared to be genetically determined. However, the detected QTLs together explained small portions of the phenotypic variation for some traits, while the heritability values were higher, suggesting that all the genetic variation is not explained by those QTL. Several assumptions may support these results. The small size of the population and sample size used in this study may lead to an underestimation of QTL numbers, an overestimation of QTL effects, and biases in the estimated proportions of the variance. Thus, some QTL with low individual effects may remain undetected (Vales 2005;Utz et al. 2000).
The diversity of QTLs among the three populations illustrates the importance of using multiple populations in QTL studies Relatively more QTLs were detected in the Meco-1 × Birmash population than that of  76T1-23 × Baji, and 76T1-23 × Birmash, with 40.00% of the QTL for yield and drought tolerance related traits identified from this population, perhaps due to the relatively higher genetic distance between the two parents. In comparison, the 76T1-23 × Baji and 76T1-23 × Birmash populations contributed fewer QTLs (25.71% and 34.28% of the total, respectively). The observation that each population contributed QTLs for drought tolerance related traits further supports the importance of using multiple populations for QTL studies that involve complex traits. Agronomic traits are a complex traits and its genetic architecture is expected to involve multiple loci and potentially multiple alleles at each locus (Yu et al. 2008). Hence, any single bi-parental population is unlikely to segregate for all genetic loci influencing the trait.
In the present study, fewer associations were detected above the significance threshold for DF than the other yield related traits, which suggests that variation in DF is controlled by either fewer larger-effect loci or true associations failed to reach the significance level. A highly significant peak (qDF7-1) on Chr.7 for DF was found in an overlapping position with previously reported QTLs such as QDTFL7.14 (Higgins et al. 2014) and QDTFL7.18 (Bouchet et al. 2017).Two of the flowering QTLs; qDF4-1and qDF8-1 were also detected in the same position with previously reported QTLs QDTFL4.21 (Felderhoff et al. 2012) and QDTFL8.19 (Wang et al. 2014), respectively. The identification of these consistent QTLs in different genetic backgrounds provides evidence that common genetic loci exist in these regions which determine flowering time in sorghum. Overall, there is a possibility that the primary QTLs responsible for flowering time in sorghum under moisture stress might be located on Chrs. 2, 4, 7, 9 and 10; and the QTLs may also refer to the same genomic regions.
The largest effect QTL (qDM6-1) for days to maturity was detected in different genetic backgrounds with stable expression, explaining 12% of the total phenotypic variance. It was detected at 9.4 kb downstream of the putative Ma6 gene on Chr. 6, a repressor gene of flowering in long days (Murphy et al. 2014). Several studies have also reported QTLs controlling maturity in sorghum near this position (Shazia et al. 2013;Nagaraja et al. 2013;Yousra et al. 2011). Three significant QTLs with varying magnitudes of effect were detected consistently on Chrs. 3, 4 and 6 in two populations. The consistency of these QTLs across the two environments and higher phenotypic variance explained provides added confidence in their real value for future breeding work. However, DF, days to flowering; DM, days to maturity; PH, plant height; LS, leaf senescence; PW, panicle weight; GYP, grain yield per plant; TSW, thousand seed weight; R(%), percentage of phenotypic variance explained by the QTL a few QTLs were located at positions not previously associated with maturity, perhaps representing novel QTLs that were not previously known to play a role in drought tolerance. A days to maturity QTL (qDM9-2) and days to flowering QTL (qDF9-1) were mapped almost to the same region/co-located on Chr. 9, indicating that these two traits could be controlled by the same gene or genes located close to each other (Zhang et al. 2013). The co-localization of QTLs for these two important traits suggests that improvement of DF may also result in enhanced maturity under drought stress. A highly significant GYP QTL (qGYP3-1), on Chr.3 was consistently identified in the two environments and overlapped with previously reported QTLs such as QGYLD3.14 (Hufnagel et al. 2014) and near to QGYLD3.16 (Nagaraja et al. 2014). Another region with strong evidence of effect on GYP is on Chr. 6 (qGYP6-1), overlapping with QGYLD6.27 (Gelli et al. 2016) and QGYLD6.21 (Sukumaran et al. 2016) previously reported. Despite the complex genetic architecture of grain yield, this comparison provided evidence for the presence of stable QTLs across different environments and genetic backgrounds. Five QTLs for yield and its component traits are clustered with positive alleles contributed from parents Baji and Birmash, possibly indicating roles of pleiotropic genes in the expression of all these traits. Fine mapping of this locus is necessary to identify the causative gene(s). Several population-specific QTLs were also detected from the three populations which could be exploited through marker-assisted QTL pyramiding to accumulate alleles conferring wider adaptation to similar environments.
A QTL (qLS2-1) for leaf senescence on Chr. 2 was consistently identified in two environments exactly at the same position in population 76T1-23 × Birmash. Similarly, Almeida et al. (2014) detected two large-effect QTLs on  Mb) and Chr.10 (120.54-146.55 Mb) for LS, explaining sizeable phenotypic variance (13 and 21%, respectively). These QTL were detected near to the present QTLs, indicating the importance of these regions in controlling leaf senescence. Moreover, some QTLs identified for LS were co-located with DF and DM QTLs that were consistently identified in more than one location from previous studies (Johnson et al. 2015;Tao et al. 2000), suggesting that these QTL can also be considered consistent since these were highly related.
Recently, two QTLs on Chrs. 10 and 4 explaining 10.5% and 12% of the phenotypic variance, respectively were reported (Nagaraja et al. 2014;Tao et al. 2000), which supports the present findings. Different researchers have also reported QTLs for LS on the same linkage group which could be the same gene although it is difficult to match both regions due to lack of common markers. Therefore, there is an indication that QTLs responsible for leaf senescence could be located on Chrs. 2 and 10 based on the observations from individual population QTL analysis.
A reliable mQTLs were detected for gene cloning and MAS in sorghum The three independent genetic mapping studies generated 105 QTLs associated with seven traits, distributed over ten chromosomes. Since all three mapping populations were segregating for the traits, it was worthwhile placing all QTL identified in the individual maps onto the consensus map. The map positions of those QTLs detected from the three populations were variable due to many factors including genetic background and population size, exemplifying that a single study can only be taken as suggestive, unless it is based on a large set of experiments. Meta QTL analysis helps to identify the most precise and concise QTLs based on the result from individual populations that can be further pursued for MAS or to predict candidate genes. The meta-analysis reduced the total number of QTLs from 105 to 25 mQTLs. The reduction in QTL number could be because of the reduction in the length of the mean confidence interval from the individual QTLs to the meta-QTLs.
Several mQTLs were detected near known sorghum maturity genes. The most significant mQTL (mqDM6-1) for days to maturity was detected 1.2 kb upstream of the putative Ma6 gene (CONSTANS-like 4; Sobic.006G004400) on Chr. 6. Ma6 was tentatively identified as SbGHD7, a repressor of flowering in long days (Murphy et al. 2014). Two significant peaks, one for days to flowering (mqDF7-1) and other for leaf senescence (mqTSW7-1) were detected near the DREB1A gene (Sobic.007G181500), which encodes a dehydration-responsive element-binding transcription factor on Chr. 7 ( (Liu et al. 1998).The mQTL peaks were located 1117.20 kb and 11.45 kb upstream of the gene, respectively. One association for GYP (mqGYP3-2) was only 215.78 kb from P5CS2, a gene that encodes an enzyme responsible for proline biosynthesis (Kishor et al. 1995), and was found to be highly expressed in a stay green sorghum line compared with a senescent line (Johnson et al. 2015). Another peak association for leaf senescence (mqLS3-1), was only 151.13 kb away from SbGI (Harmon et al. 2018), encoding a Gigantea protein (GI) involved in flowering time control and a wide range of other physiological activities.
From the present study associations near to Ma6 were detected, a major flowering gene in sorghum, implied the usefulness of these populations and associated data. In addition, several associations were detected near to known sorghum maturity genes DREB1A gene, which encodes a dehydrationresponsive element-binding transcription factor (Liu et al.1998), P5CS2, a gene that encodes an enzyme responsible for proline biosynthesis (Kishor et al. 1995) and SbGI (Harmon et al. 2018), a gene that encode a Gigantea protein (GI) involved in a wide range of physiological activities. Hence, this is a confirmation for the reliability of the QTLs detected in our study when used for gene cloning and MAS.
We also identified several associations that did not appear to be associated with previously reported QTLs or known candidate genes, and may represent novel genes not previously known to play a role in determining those traits. Alternatively, some might be explained by structural variation in the sorghum genome, resulting in annotated genes from the reference genome being in different genomic locations within the sorghum F2:3 populations, thus putting the association peaks in new locations. Those consistent QTLs, which can be regarded as hotspots with agronomical importance, are attractive regions for further study to identify closely-linked DNA markers for the causative genes involved in the genetic control of drought tolerance for candidate gene analysis and marker-assisted breeding.

Conclusion
We identified 105 QTLs responsible for sorghum yield and drought tolerance related traits using three connected populations grown under moisture stress conditions. By constructing a consensus linkage map, 25 mQTLs were identified that are common across the three genetic backgrounds. We confirmed the reliability of these associations through multi-faceted analyses, including comparing to previous QTLs and known genes. Several mQTLs were detected near known sorghum maturity genes and QTLs previously reported. In the future, it is of great interest to clone genes underlying mQTL regions and uncover molecular mechanisms of sorghum drought tolerance.