Genome-wide identification of quantitative trait loci for morpho-agronomic and yield-related traits in foxtail millet (Setaria italica) across multi-environments

Foxtail millet (Setaria italica) is an ideal model of genetic system for functional genomics of the Panicoideae crop. Identification of QTL responsible for morpho-agronomic and yield-related traits facilitates dissection of genetic control and breeding in cereal crops. Here, based on a Yugu1 × Longgu7 RIL population and genome-wide resequencing data, an updated linkage map harboring 2297 bin and 74 SSR markers was constructed, spanning 1315.1 cM with an average distance of 0.56 cM between adjacent markers. A total of 221 QTL for 17 morpho-agronomic and yield-related traits explaining 5.5 ~ 36% of phenotypic variation were identified across multi-environments. Of these, 109 QTL were detected in two to nine environments, including the most stable qLMS6.1 harboring a promising candidate gene Seita.6G250500, of which 70 were repeatedly identified in different trials in the same geographic location, suggesting that foxtail millet has more identical genetic modules under the similar ecological environment. One hundred-thirty QTL with overlapping intervals formed 22 QTL clusters. Furthermore, six superior recombinant inbred lines, RIL35, RIL48, RIL77, RIL80, RIL115 and RIL125 with transgressive inheritance and enrichment of favorable alleles in plant height, tiller, panicle morphology and yield related-traits were screened by hierarchical cluster. These identified QTL, QTL clusters and superior lines lay ground for further gene-trait association studies and breeding practice in foxtail millet. Supplementary Information The online version contains supplementary material available at 10.1007/s00438-022-01894-2.


Introduction
Foxtail millet (Setaria italica), a member of genus Setaria, subtribe Cenchrinae, tribe Paniceae, subfamily Panicoideae, family Poaceae (Kellogg. 2017), has been domesticated and cultivated extensively as a staple food crop for over 11,500 years in China (Wang et al. 2012;Hunt et al. 2021). Foxtail millet has also been developed into a model organism for studying architectural evolution of monocotyledon plant (Doust. 2007;Mauro-Herrera et al. 2016), C 4 photosynthesis (Brutnell et al. 2010;Mamidi et al. 2020;Chatterjee et al. 2021), nutritional properties (Muthamilarasan et al. 2016), abiotic tolerance and bioenergy in cereal grasses (Kole et al. 2015). In the process of domestication from its ancestor green foxtail (Setaria viridis), foxtail millet has undergone a strong artificial selection to generate a wide range of phenotypic changes in branching, plant height, shattering, flowering time and seed size . One of the most critical steps Communicated by Bing Yang.
1 3 in its domestication was the retention of panicle integrity (i.e. non-shattering) (Schroder et al. 2017;Odonkor et al. 2018). The panicle traits, determining its inflorescence architecture, include primary branch number and density, primary branch length, numbers of branching orders, spikelet number and bristle (sterile branchlet) number, which are also the key morphological characters distinguishing S. italica and S. viridis (Doust et al. 2002(Doust et al. , 2005Hussin et al. 2020). Identification of QTL and genes underlying these phenotypes is essential to elucidate genetic mechanism of development of Setaria and further breeding application with genome editing tools. Then interspecific segregation population helps to mine QTL/ genes of phenotypes that have been selected significantly during domestication. From F 2 , F 2:3 and RIL populations of an interspecific cross between S. italica acc. B100 and S. viridis acc. A10, QTL for basal branching (tillering), axillary branching, inflorescence architecture related traits, shattering, flowering time, culm height, number of nodes and plant biomass were meticulously detected, and a few genes, such as tb1, SD-1, Sh1 and qSH1, were identified and evaluated by comparative mapping with other cereal crops (Doust et al. 2004(Doust et al. , 2005(Doust et al. , 2006Mauro-Herrera et al. 2013;Odonkor et al. 2018). These QTL and putative genes identified in S. viridis served as the basis for foxtail millet domestication.
The existing landraces and contemporary bred cultivars display different panicle types including cylindricalshaped type, conical-shaped type, spindle-shaped type, and tip-branched type, which could be a result of the variation in primary branch length and density in different parts of the inflorescence. Gene NEKODE1 responsible for the tip-branched panicle is mapped to the physical position around 13.6-14.4 Mb on chromosome 9 (Masumoto et al. 2016). Furthermore, Hussin et al. (2020) found that a novel member of MADS-box transcription factors, SiMADS34, involves in the regulation of panicle width, primary branch length, number of primary branches, panicle length and grain weight. Apart from inflorescence branching, QTL of agronomic and yield-related traits, namely heading data (HD), panicle weight (PW), panicle length (PL), panicle diameter (PD), flag-leaf length (FLL), plant height (PH), stem diameter (SD), stem node number (SNN), code number (CN), code grain number (CGN), thousand-grain weight (TGW), neck length (NL), leaf color (LC), bristle color (BC) and anther color (AC) have been positioned on a high density map derived from different bi-parents intraspecific population Fang et al. 2016;Zhang et al. 2017;Wang et al. 2017Wang et al. , 2019b. Among them, QTL for plant height and heading time were elucidated in detail through linkage and bulked segment analyses (BSA) and the relevant genes, such as Heading date 1 (Hd1), KELCH REPEAT,, Roc4 and Seita.1G242300, were predicted by homologous comparisons with close species (Mauro-Herrera et al. 2013;Yoshitsu et al. 2017;Jaiswal et al. 2019;He et al. 2021).
In China, landraces or bred cultivars of foxtail millet has the striking ecotypes or eco-regions that can be classified into the Northeast Plain, North China Plain, Inner Mongolia Plateau and the Northwest Plateau according to the natural climatic conditions of the foxtail millet-producing areas and the heading dates of various groups of foxtail millet varieties (Wang et al. 2012;Jia et al. 2015;Diao et al. 2017). Thus, local adaptation is an important factor in foxtail millet evolution that should be considered in the species breeding programs. Siprr37 with a transposon insertion was identified as a gene responsible for the adaptation of foxtail millet to the environmental conditions of the early spring sowing region (the Northeast Plain) . Notably, using varieties from different eco-regions to construct a segregation population and dissect quantitative trait loci/gene is necessary and imperative for foxtail millet breeding.
In the present study, we used two cultivars, one from the Northwest Plateau eco-region and the other from the North China Plain eco-region, to cross and construct the RIL population which then was used for genome-wide resequencing and constructing an updated high-density bin map. Seventeen main morpho-agronomic and yield-related traits in breeding practice were phenotyped under four to ten environments. Combined with phenotypes and genotypes, single environment QTL, multi-environment QTL, QTL clusters and six superior lines were identified, revealing the important genomic regions of 17 traits in foxtail millet. These results laid a foundation for fine mapping, identification of candidate genes, elaboration of molecular mechanism of development and breeding application in foxtail millet.

Plant materials
The progenies of the RIL population derived from a cross of Yugu1 and Longgu7, first reported in Liu et al. (2020), were further trialed and phenotyped in 2018, 2019, 2020 and 2021. The F 2:10 individuals and the parents were grown at Sanya (SY) in the winter of 2018 to the spring of 2019 and two different irrigated fields (DHa: irrigation prior to sowing and at the seedling stage; DHi: irrigation prior to sowing and at the seeding, jointing and filling stages) at DH in 2019, respectively. In 2020, the same sites at DH were used for testing with F 2:11 lines. For investigating tiller of the RIL population, F 2:12 were planted in two environments similar to DHa and DHi at Dunhuang and two field (TGh: soil environment with a high phosphorus content; TGl: soil environment with a low phosphorus content) at Taigu (TG) in 2021. The year, irrigation and soil phosphorus conditions in different geographic locations were combined to form 14 test environments (Fig. S1). The geographic location, elevation and other relevant information of the test sites were attached in Fig. S1. All data before 2019 were combined for subsequent analyses.

Phenotyping
Sequential arrangement method in single factor experimental design was adopted in all test environments and the arrangement order was Yugu1, Longgu7, RIL1, RIL2, …, RIL10, Yugu1, Longgu7, RIL11, RIL12, …, RIL20, Yugu1, Longgu7, … …, Yugu1, Longgu7, RIL141, RIL142, …, RIL164, Yugu1, Longu7 (Fig. S2). In total, 164 RIL lines were used. ~ 80 plants of each seed line were plants. The planting density was 450,000 individuals per hectare. Morpho-agronomic traits were classified into four patterns: 1) stem traits (length of the main stem (LMS), diameter of the main stem (DMS), node number of the main stem (NMS), peduncle length (PL), tiller number (TN)), 2) leaf traits (flag leaf length (FLL) and flag leaf width (FLW)), 3) paniclerelated traits (main panicle length (MPL), main panicle diameter (MPD), spikelet density (SD), grain number per spikelet (GNS) and bristle length (BL)), 4) growth period (from emerging to ripen (GP)) and 5) yield-related traits including straw weight per plant (SWP), panicle weight per plant (PWP), grain weight per the main panicle (GWP) and 1000-grain weight (TGW). At the end of filling stage, flag leaf length, flag leaf width and bristle length of 10 plants were measured with a ruler with 1 mm accuracy. After the plant ripening, 15 plants per line were selected randomly and pulled out manually with roots to measure other traits. The diameter of the main stem and main panicle diameter were measured by electronic vernier calipers with accuracy of 0.01 mm. Other length traits were measured with a ruler with accuracy of 0.1 cm. The tiller number, node number of the main stem and spikelet number per panicle were counted visually. Spikelet density was calculated by dividing spikelet number per panicle by main panicle length. After plant drying out, ten spikelets per line were randomly selected to measure grain number per spikelet. The 1000 seeds counted using automatic granule counting machine were weighed to get individual seed weight. Then 15 labeled individual plants were weighed together for straw weight per plant, panicle weight per plant, and grain weight per main stem.

High-throughput sequencing, sequence alignment and variant calling
The same batch of high-throughput sequencing of the parental lines and 164 F 2:8 RILs from NCBI under an SRA accession number PRJNA562988 were used as we reported earlier (Liu et al. 2020). More strict analyses were used in this study than previous described (Liu et al. 2020). Briefly, the raw data were filtered by cutadapt (version 1.13) and trimmomatic (version 0.36) software to remove the residual adapter and low-quality sequence (Q score < 30, reads length < 50 bp). Then, the high-quality reads were aligned to the foxtail millet reference genome S. italica v2.0 (http:// plants. ensem bl. org/ Setar ia_ itali ca/ Info/ Index) using BWA (version 0.7.15-r1140) with MEM algorithm. SamTools software (version 1.3.1) was used to convert the alignment results into BAM format. And SortSAM in Picard tool (version 1.91) was applied to sort the reads in BAM files. Then, using RMDUP in SamTools removed PCR duplication, resulting in a BAM file for further coverage and coverage depth statistics and variant calling. GATK (version 3.7) software package was employed to conduct variant calling for all samples, including SNP and InDel. The resulting variation was further screened according to the following conditions: 1) The proportion of missing genotype in individuals is less than 25%; 2) The frequency of minor alleles is not less than 20%; 3) Observed heterozygosity of individuals is less than 25%.

Genotyping and construction of linkage maps
Hidden Markov Model (HMM) described by Xie et al. (2010) was adopted to perform genotyping of RIL. Briefly, the parental genotypes were inferred by the linkage relationship between SNP/InDels in the RIL population, and genotypes of each RIL were converted to A and B via inferred genotypes of Yugu1 and Longgu7. Simultaneously, the actual genotypes of Yugu1 and Longgu7 were used to judge the authenticity of the inferring parental genotypes. Based on HMM, the missing genotypes were filled, and the wrong genotypes were corrected. Then the identical SNP and InDel genotypes in the interval were merged to generate bin markers. Furthermore, 74 SSR unevenly distributed on 9 chromosomes of foxtail millet reported previously by Fang et al. (2016) were selected to assess the accuracy of the bin genotype. Then the genetic map distance between markers including bin and SSR was calculated using the Kosambi mapping function in R/onemap (version 2.1.3) and the linkage map was constructed through applying R/Link-ageMapview (version 2.1.2).

QTL identification
MapQTL 6.0 was applied to detect QTL with multiple-QTL models (MQM) mapping. A threshold of log of odds (LOD) ≥ 2.0 indicated the existence of QTL. 1-LOD supporting interval was designated for the confidence intervals (Wang et al. 2019a). The sign " + " represented that allele was from Yugu1, while "-" indicated that allele derived from Longgu7. QTL with partially or fully overlapped supporting intervals were regarded as the identical QTL, and QTL detected under two or more test environments was termed as the stable QTL. The letter "q" combining with the trait abbreviation, the chromosome number and the QTL serial number was used to represent the QTL identity.

Statistical analysis
Descriptive statistics and analysis of variance (ANOVA) were performed for each trait using SPSS Statistics 17.0. Pearson correlation coefficients between traits in different trials were calculated and visualized using R/corrplot package (version 0.90). Package ggplot2 (version 3.3.5) and NbClust (version 3.0) in R language were used to perform the distribution and cluster analysis for all traits, respectively. In the cluster analysis, the function scale was used to standardize the phenotypic values of 17 traits, and then the Nbclust package performed the hierarchical cluster analysis, in which distance and method arguments were euclidean and average, respectively (Charrad et al. 2013). Lem4 package (version 1.1-23) in R language was applied to count heritability for each trait, and the formula is as follows: where v g , v ge , v gy and v e are genetic variance, variance of interaction between genotype and environment, variance of interaction between genotype and year and environmental variance, respectively. L and Y are the number of environments and years, respectively.

Phenotypic variation
All traits were classified into five categories based on plant different organs and growth habit, namely stem, leaf, panicle, growth period and yield-related traits. Among stemrelated traits, LMS, DMS, NMS, PL and TN were measured separately in 10, 9, 9, 8 and 8 testing environments, respectively. The values of these traits in general increased with the increase of irrigation and growth period. Almost all the traits conformed to the normal distribution and were typical quantitative traits (Table S1, Fig. S3). However, growth period and tiller number appeared more different under all tested environments. Growth period was longer than other trails in HN, and displayed two peaks under DHa and DHi in 2020, suggesting that it may be regulated by multiple majorminor effect loci. Mean of tiller number was higher at WW, SY and TG than other sites. Four yield-related traits were greater at HN and DH than other tested trials. Panicle weight per plant and grain weight per plant were lowest at SY than other places. On stem related traits, length of the main stem and nodes of the main stem were greater at WW than other environments. All traits are significantly affected by environment factors and a significant (P < 0.01) interaction between genotype and environment (Table S2). But the relationships among all traits were significantly different. There were significantly positive relationships between yield-related traits and panicle-related traits, such as panicle length and the main panicle diameter, and a common negative relationship between grain number of spikelet and spikelet density on all tested environments (Fig. 1, Fig. S4). Only flag leaf length and flag leaf width were measured for leaf-related traits in nine of the fourteen tested environments. Both traits were significant positively and higher at DHi in 2019 and 2020. Among all traits, heritability of MPD, LMS and SWP was higher than 0.9. GP, FLL, PL, MPL and SD were between 0.8 and 0.9. The range of FLW, DMS, NMS, GNS, PWP and TGW was 0.7 to 0.8. BL and GWP were lower than 0.7. The heritability of TN was the lowest (Table S3).

Sequencing, SNP/InDel identification, and bin map construction
Resequencing for 164 F 2:8 lines and two parental lines was completed previously (Liu et al. 2020). After filtering, there was a total of 1,929,632,243 clean reads with average of 11,624,291 reads per line (Table S4). The 87.8% of the reads were uniquely mapped to the reference genome, with an average 7.45 × depth and 92.8% genome coverage for each line (Table S4). In total, 1,325,599 SNPs and 167,818 InDels were identified from all the samples (Table S5). The number of SNPs and InDels on chromosomes was decreasing in the following order: chromosome 8, 7, 2, 3, 9, 5, 6, 1 and 4, indicating the corresponding order of chromosome genetic differences between two parental lines (Fig. 2a, Table S5). Over 70% of variations were in intergenic regions, and only 11.1% in genic regions. Remaining 15.3% and 1.9% of variations were in upstream/downstream and UTR, respectively (Fig. 2b). The application of the HMM algorithm generated 2,297 bin markers. A bin map from these markers was constructed using R/onemap (Table 1, Fig. 2c). Then, R/LinkageMapveiw was used to construct a linkage map with these bin markers and 74 SSR markers, covering a total of 1315.1 cM on the nine linkage groups of Setaria.
The interval between adjacent markers ranged from 0.3 to 13.5 cM, with an average interval of 0.56 cM (Table 1, Fig. S5). Analysis of pairwise recombination fractions and LOD scores showed that most of markers had no aberrant phenomenon on the bin map except for few co-segregation markers (Fig. 2d, Fig. S5).

QTL mapping
A total of 447 loci formed 221 QTL for seventeen traits were detected across multi-environments (

QTL for flag leaf length
For flag leaf length (FLL), a total of 15 QTL were identified in nine environments and mapped on all chromosomes (Fig. 3), explaining 5.5 to 17.3% of the phenotypic variance. Twelve QTL regions were repeatedly identified in at least two environments, including the only one QTL qFLL2.2 detected in six environments. The Yugu1 alleles had a positive effect for 10 QTL identified except for qFLL1.1, qFLL1.2, qFLL5.3, qFLL7.2 and qFLL7.3 (Table S6).

QTL for tiller number
There were 9 tiller QTL identified at 8 environments, positioning on chromosome 1, 4, 5, 7, 9, and explaining 6-17.6% of the phenotypic variance. Of these, qTN5.2 and qTN7.3 were the two most stable QTL for tiller detected in six environments. qTN5.1 and qTN7.1 were identified in four and three environments, respectively. qTN1.1 and qTN7.2 were the two dual-environments QTL. And the remaining 3 QTL were identified in a single environment. Alleles increasing tiller number were from Yugu1 except for qTN7.1 and qTN9.1 (Fig. 3, Table S6).

QTL for peduncle length
Fourteen QTL for peduncle length (PL) were detected in eight environments and distributed on chromosome 1, 2, 5, 6, 7, and 9 (Fig. 3, Table S6). The phenotypic variation explained by these QTL ranged from 5.8 to 18%. qPL5.2, qPL1.1 and qPL6.2 were identified in five, four and three environments, respectively. qPL5.3 and qPL6.1 were detected in two environments. And the remaining 9 QTL were detected under only one environment. Except for qPL1.1, the additive effects for the other QTL were contributed by the Yugu1 (Table S6).

QTL for length of the main stem
Seventeen QTL for the length of the main stem (LMS) were detected in ten environments and mapped on all chromosomes, accounting for 5%-36% of phenotypic variance. Among these, ten QTL including the most stable QTL (qLMS6.1) in this study were detected in at least two environments, while the others were environmentspecific QTL. The Longgu7 alleles had a positive effect for qLMS1.1, qLMS2.1, qLMS4.1, qLMS4.2, qLMS5.1, qLMS5.2, qLMS6.1 and qLMS9.2, whereas the additive effects for the other QTL were contributed by the Yugu1. (Fig. 3, Table S6).

QTL for diameter of the main stem
Fourteen QTL associated with diameter of the main stem (DMS) were detected in six environments and mapped on chromosomes 1, 5, 6, 7 and 9, and explained 5.6-15.1% of the phenotypic variation. Of these QTL, qDMS9.4 was identified across four environments. qDMS5.4 and qDMS6.2, qDMS5.1 and qDMS5.2 were detected three, two environments, respectively. The remaining QTL were identified in only one environment. The additive effect of all the QTL originated from Yugu1, with the exception of qDMS5.1 (Fig. 3, Table S6).

QTL for node number of the main stem
Thirteen QTL for node number of the main stem (NMS) were detected in seven environments and located on chromosome 1, 2, 6, 7, 8 and 9 (Fig. 3). The explained phenotypic variation ranged from 5.6 to 22.8%. Among them, qNMS1.3 and qNMS6.3 were detected in four and three environments, respectively. Seven QTL were identified in two environments, and the others were environmentspecific QTL. The additive effects for these QTL were contributed by the Yugu1 alleles except qNMS1.1, qNMS1.2, qNMS2.1 and qNMS6.3 (Table S6).

QTL for main panicle diameter
Fourteen QTL for main panicle diameter (MPD) were identified across nine environments and located on all chromosomes except chromosomes 2 and 3, explaining 5.9%-23.7% of phenotypic variance. Of these, qMPD5.1 and qMPD5.2 were detected across six environments, and qMPD9.3, qMPD5.3, qMPD9.1 were detected in five, four and three environments, respectively, which all of the additive effects were from Yugu1. qMPD6.1 and qMPD7.1 were detected in two environments. The additive effect of qMPD6.1 was derived from Yugu1, while qMPD7.1 was from Longgu7. The remaining 7 QTL were identified in a single environment, and all of the additive effects were derived from Yugu1 except qMPD7.2 (Fig. 3, Table S6).

QTL for spikelet density
Fourteen QTL for spikelet density (SD) were detected across four environments, and explained 5.5-11.9% of the phenotypic variation of these QTL, two (qSD3.1 and qSD3.2), two (qSD1.1 and qSD6.1), three (qSD8.2, qSD9.1 and qSD9.2) were identified in five, three and two environments, respectively. The favorable allele increasing SD was from Yugu1. The remaining 7 QTL were only identified in a single environment. The additive effects for these QTL were contributed by the Yugu1 alleles, with the exception of qSD7.1 (Fig. 3, Table S6).

QTL for grain number per spikelet
Eight QTL for grain number per spikelet (GNS) were detected across three environments and mapped on chromosome 1, 3, 5, 8 and 9, which accounted for 5.8-13.2% of the phenotypic variation. Among them, qGNS8.2 and qGNS9.1 were identified in two environments, while the other six QTL were detected under a single environment. The additive effects of qGNS5.3 and qGNS9.1 were derived from Yugu1 and the other QTL were from Longgu7 (Fig. 3, Table S6).

QTL for bristle length
Five QTL for bristle length (BL) were detected in five environments and distributed on chromosome 1, 2, 7, 8 and 9 (Fig. 3). The phenotypic variation explained by the five QTL ranged from 6.3 to 35.2%. qBL1.1 and qBL7.1 were identified in five and two environments, respectively. And the remaining QTL were detected in a single environment. The additive effects of qBL1.1 and qBL7.1 were from Longgu7, while the other QTL were derived from Yugu1. qBL1.1 had the highest percentage of phenotypic variation and was detected in all test environments, which was an important locus to dissect the genetic mechanism of bristle length.

QTL for panicle weight per plant
Fifteen QTL for grain weight per plant (PWP) were detected in nine environments and mapped on the other chromosomes except for chr.8, explaining 5.6-19.5% of the phenotypic variance (Fig. 3, Table S6). Ten QTL were identified in the single environment. qPWP2.2, qPWP9.2 and qPWP5.3, qPWP6.1 were detected in two and three environments, respectively. qPWP5.2 was identified across five environments. The favorable alleles of qPWP2.1 and qPWP2.2 for increasing the trait value came from Longgu7, while favorable alleles of the other QTL were derived from Yugu1. Compared with the published (Liu et al. 2020), qPWP5.2 were identified for a new and more stable QTL controlling panicle weight per plant (Table S6).

QTL for grain weight per the main panicle
Fourteen QTL for grain weight per the main panicle (GWP) were identified on eight environments and unevenly distributed on eight chromosomes except for Chr.8, which explained 6.1-15.1% of the phenotypic variance (Fig. 3, Table S6). Of these QTL, eight were detected in the single environment, four (qGWP2.1, qGWP5.3, qGWP6.1 and qGWP9.1) and two (qGWP5.1 and qGWP5.2) were detected in two and five environments, respectively. The favorable alleles of four QTL (qGWP2.1, qGWP2.2, qGWP7.1 and qGWP9.2) came from Longgu7, while the favorable alleles of the remaining QTL were derived from Yugu1. Compared with the published (Liu et al. 2020), after updating the linkage map and adding the test environments, two new stable QTL (qGWP5.1 and qGWP5.2) were mapped on chromosome 5 (Table S6).

QTL clusters for multiple traits
The QTL cluster was defined as a chromosome region harboring multiple QTL for various traits within ~ 20 cM (Tan et al. 2018). And phenotypically correlated traits are often mapped to similar genomic regions .
Of 221 QTL for seventeen traits, 130 with overlapping intervals formed 22 QTL clusters (Table 2). Chromosome 9 and 5 carried separately five and four QTL clusters and chromosomes 2, 3, 6 and 7 had two QTL clusters, whereas there were three and one clusters on chromosomes 1 and 4, 8, respectively (  (Table 2). These paired traits in the same cluster had significant correlation in the most of test environments (Fig. S4), suggesting there was the 1 3 existence of pleiotropy or tight linkage in genomic regions of these QTL clusters.

Superior lines
To select out dwarfing lines and lines with excellent comprehensive characters, we performed cluster analysis with plant height and peduncle length as one kind of variables and the other 15 traits as another kind of variables. Six Fig. 4 The stable QTL identified in two to six combined environments recombinant inbred lines, including RIL35, RIL48, RIL77, RIL80, RIL115 and RIL125 with transgressive inheritance in plant height and peduncle length, panicle morphology and yield related-traits were identified as superior lines (Fig.  S6). Of these, RIL35 was classified into a distinctive branch according to the length of the main stem and peduncle length, both of which had significantly positive correlation in most tested environments (Fig. 1, S4, S6a). Furthermore RIL48, RIL77, RIL80, RIL115 and RIL125 with better comprehensive characters were clustered into another distinctive branch via the remaining 15 phenotypic values (Fig. S6b, Table S7). T-test showed that length of the main stem and peduncle length of RIL35 were significantly shorter than those of parents (Fig. 5). RIL35 carried dwarf genotypes on all 17 QTL regions for length of the main stem and 11 out of 14 QTL regions for peduncle length (Fig. S7a). Similarly, RIL115 having obvious over-parent characteristics in straw weight per plant, panicle weight per plant and grain weight per plant, combined 38 favorable alleles of 43 QTL regions for the three traits (Fig. 5, S7f). RIL125, an over-parent line in spikelet density and tiller, carried 13 of the 14 QTL alleles decreasing spikelet density, 6 of 9 QTL alleles increasing tiller number, respectively (Fig. 5, S7b). RIL77 combined 11 favorable alleles in 14 QTL regions to produce over-parent phenotype in length of main panicle (Fig. 5, S7c). The bristle length of RIL48 and RIL 80 were significantly different from those of parents (Fig. 5). The two lines carried opposite genotypes in the most stable QTL (qBL1.1) and a single environment QTL (qBL8.1) regions, indicating that the two QTL were the most major loci affecting bristle length ( Fig.  S6 d and e). These results suggested that over-parent characteristics of these traits in superior lines were resulted from an accumulation of favorable alleles.

An updated linkage map
The previous linkage map was constructed based on F 2 and RIL population from a cross between Yugu1 and Longgu7 (Fang et al. 2016;Liu et al. 2020). In the present study, we reanalyzed and used a more stringent selection of higher quality SNPs than previously published

Correlations among traits
Seventeen morpho-agronomic and yield-related traits are the most important in foxtail millet breeding and also the pivotal indicators distinguishing different foxtail millet ecotypes. Interrelationships between them affect foxtail millet morphogenesis and application in practice. QTL and genes for the traits in the same plant organ or those which are intercorrelated are more likely to be located in the same genomic region. In foxtail millet, Zhang et al. (2017) positioned the QTL for PL, FLL, and PH that all clustered on the chromosomes 5 and 9. Zhi et al. (2021) identified 34 co-located QTL clusters involving nine traits related to panicle architecture and grain yield. Wang et al. (2019b) suggested that intercorrelated traits were probably associated with the growth module and genetic regulation pathway. In this present study, we identified a similar trend and correlation, for example, with cluster1.

Favorable QTL alleles and superior lines
The favorable alleles for a trait do not necessarily come from the more favorable parent (Wang et al. 2019a). In foxtail millet QTL mapping, Wang et al. (2017) identified that the additive effects of 7 out of the 11 major QTL for plant height, main panicle length, main panicle diameter, first main internode diameter, second main internode diameter, and third main internode diameter and main panicle weight per plant were from inferior parent. Among the 26 QTL for plant height identified by He et al. (2021), the dwarfing favorable alleles of 8 QTL came from the parent with higher plant height. Zhi et al. (2021) identified 159 QTL for panicle architecture and grain yield-related traits, of which the additive effects of 60 QTL were from Ai88 and the additive effects of 89 QTL were from Liaogu1. In the present study, the 74 and 147 QTL alleles increasing phenotypic values originated from Longgu7 and Yugu1, respectively. These results indicated that both the superior and inferior parent could contribute QTL alleles that increase the trait values. The accumulation of parental favorable alleles leads to the generation of over-parent lines in progeny populations (Wang et al. 2019a). Kulkarni et al. (2020) identified that six RILs possessing the major yield related QTL and fertility restorer loci Rf3 and Rf4 alleles were complete restorers, which can be useful in hybrid rice breeding. Wang et al. (2020) identified that phenotypically superior RIL47 with the major QTL genotypes linked to rice photosynthesisrelated traits (PRT) could be considered for genetic improvement of PRT under cold water stress. In this study, six superior recombinant inbred lines, RIL35, RIL48, RIL77, RIL80, RIL115 and RIL125 combining multiple favorable alleles from bi-parents and generating transgressive inheritance in plant height, tiller, panicle morphology and yield related-traits were identified for superior lines, which could be applied into foxtail millet breeding programmes.

The necessity of fine mapping for the stable QTL
One of the main aims of QTL mapping is to identify and clone causal genes affecting target traits (Liu et al. 2012). QTL mapping is affected by the resolution of genetic map, the size of the population, and the accuracy of the phenotyping. The construction of high density or ultra-high density improved the precision of QTL mapping in foxtail millet Wang et al. 2019b;He et al. 2021;Tian et al. 2021). However, the amount of variations between parental lines and the degree of uniform distribution on chromosomes further affects the resolution of the map. Therefore, chromosome segments with dense markers and small haplotype blocks are preferred for mapping QTL into relatively small intervals. Previously, based on the 333 RIL population and 3744 bin markers, He et al. (2021) mapped and predicted a plant-height QTL gene (QTG), Seita.1G242300 in the interval with dense variations through homologous comparison with relative species. In this study, utilization of high-density markers improved the resolution of mapping. However, qGP9.1 with the smallest interval on the terminal of chromosome 9 still harbored over 140 genes, so it is difficult to map and clone quantitative trait genes on the primary population with a small number of individuals or intervals with fewer variations. Fang et al. (2017) established an F 2 secondary population with 2484 individual plants from a cross between recombinant line RIL014 and CCRI35 and fine mapped qFS07.1 into a 62.6 kb genome region containing four annotated genes, which greatly improved the precision of candidate gene identification. Then, fine mapping of qGP1.1 or all stable QTL will be necessary through constructing secondary population with the larger number of individuals. Moreover, qLMS6.1, a most stable QTL within 33,367,330-34,859,035 on chromosome 6, was mapped across nine environments (Table S6, Fig. 3). The QTL for plant height was also identified in fewer environments by Jia et al. (2013), Fang et al. (2016), Zhang et al. (2017) and He et al. (2021). Seita.6G250500 underlying qLMS6.1, an ortholog to Os08g44590 in rice, is a putative candidate gene, involving in GA biosynthesis or signaling pathways . Similarly, the qLMS6.1 region contains 188 annotated genes (https:// phyto zome. jgi. doe. gov/ pz/ portal. html# !info? alias= Org_ Sital ica). whether Seita.6G250500 is the major effect gene is still uncertain, so fine mapping of the QTL will be necessary for the precise identification of causal gene and functional verification.

Conclusion
An updated high-density genetic map harboring 2297 bin and 74 SSR markers was constructed, which covered 1315.1 cM with an average distance of 0.56 cM between consecutive markers. Two hundred twenty-one QTL including 103 stable QTL and 22 QTL clusters for 17 morphoagronomic and yield-related traits were identified across 14 environments. Six recombinant inbred lines were identified as superior materials in plant height, tiller, panicle morphology and yield related-traits. This study provided insights into the genetic dissection of 17 morpho-agronomic and yieldrelated traits in foxtail millet.