Mapping of dynamic quantitative trait loci for plant height in a RIL population of foxtail millet (Setaria italica L.)

Plant height (PH) is a crucial trait for strengthening lodging resistance and boosting yield in foxtail millet. To identify quantitative trait loci (QTL) and candidate genes associated with PH, we first developed a genetic map using a recombinant inbred line (RIL) population derived from a cross between Aininghuang and Jingu 21. Then, PH phenotyping data and four variations of best linear unbiased prediction (BLUP) were collected from nine environments and three development stages. Next, QTL mapping was conducted using both unconditional and conditional QTL methods. Subsequently, candidate genes were predicted via transcriptome analysis of parental samples at three developmental stages. The results revealed that the genetic map, based on re-sequencing, consisted of 4,360 bin markers spanning 1,016.06 cM with an average genetic distance of 0.23 cM. A total of 19 unconditional QTL, accounting for 5.23%–35.36% of the phenotypic variation explained (PVE), which included 7 major and 4 stable QTL, were identified. Meanwhile, 13 conditional QTL, explaining 5.88%–40.35% of PVE, including 5 major and 3 stable QTL, were discovered. Furthermore, four consistent and stable QTL were identified. Finally, eight candidate genes were predicted through RNA-seq and weighted gene co-expression network analysis (WGCNA). Those findings provide a crucial foundation for understanding the genetic mechanisms underlying PH development and facilitate molecular marker-assisted breeding of ideal plant types in foxtail millet.


Introduction
Foxtail millet (Setaria italica L., 2n=2x=18), with a cultivation history extending over 16,000 years in China (Doust et al., 2009), is considered an important source of food and fodder globally, particularly in the arid and semi-arid regions such as India and China (Singh and Prasad, 2020).Because of its higher resistance to drought, higher adaptability to infertile soils, and higher water use efficiency, foxtail millet not only is of immense significance for the advancement of ecological agriculture (Singh et al., 2024), but also plays a key role in adjusting cropping structures as a strategic reserve crop (Diao, 2019;Ramesh et al., 2023).However, increasing yield and environmental adaptability in foxtail millet, compared to other major cereal crops like maize, remains a significant potential, by developing short-statured and lodging-resistant varieties (Diao et al., 2014).
All previously mentioned QTL and genes were identified using a conventional mapping method to analyze phenotype data at maturity.However, the development of PH is a complex and dynamic process, influenced by complex genetic networks and environment factors.Numerous genes involved in PH development exhibit dynamic expression patterns across various development stages (Che et al., 2020;Fu et al., 2022;Wu et al., 2022).To fully understand the regulatory mechanisms of PH development, it is imperative to determine the temporal and spatial expression patterns of underlying genes.Consequently, it is of paramount importance to accurately identify QTL and how to regulate PH across entire development stages.
Multiple major and stage-specific QTL for PH in crops such as wheat, cotton, and rice have been identified by dynamic QTL mapping, indicating that PH is differentially expressed throughout the plant development process (Che et al., 2020;Fu et al., 2022;Wu et al., 2022).
In foxtail millet, only Mauro-Herrera and Doust (2016) carried out dynamic QTL mapping for PH in a RIL population and found that the QTL H9a, with a PVE ranging from 7.6% to 15.5%, had a significant and consistent impact during the vegetative growth, flowering, and harvest stages under varying environmental conditions.However, research on dynamic QTL mapping for PH throughout the entire developmental process is still scarce.
In the present study, both unconditional and conditional QTL mapping methods were employed to identify QTL associated with PH.This approach will enable us to uncover QTL that are specifically expressed only under certain conditions and those stably expressed across different environment conditions and developmental stages.By further integrating of RNA-seq and weighted gene co-expression network analysis (WGCNA), we aim to identify candidate genes potentially involved in regulating PH.These findings will not only deepen our understanding of the molecular mechanisms underlying the dynamic development of PH, but also lay a foundation for the fine mapping and cloning of QTL for PH in foxtail millet.

Plant materials
To carry out QTL mapping, we developed an F 2:7 recombinant inbred line (RIL) population consisting of 127 lines.This population was created using the single-seed descent method in a cross between Aininghuang (female parent, derived from natural variation of Ninghuang 1) and Jingu 21 [male parent, accession number: GPD Foxtail millet (2017)140009, obtained from Co60-irradiated Jinfen 52 dry seeds].Aininghuang is characterized by a dwarf phenotype of approximately 110 cm, whereas Jingu 21 typically attained a taller stature, approximately 180 cm.The RIL population, along with two parents, was cultivated in three distinct experiment sites in Shanxi Province (China), during three crop seasons (2020)(2021)(2022).The sites included Datong (DT,39.3°N,113.3°E),characterized by a cold-dry climate, with growing season temperatures ranging from 15°C to 25°C, and an annual precipitation of approximately 400 mm; Jinzhong (JZ, 37.6°N, 112.7°E), characterized by a semi-arid climate where summer temperatures often exceed 30°C and annual precipitation is approximately 450 mm; and Changzhi (CZ,36.2°N,113.1°E),which has a humid subtropical climate with temperatures ranging from 20°C to 28°C and an annual precipitation of approximately 550 mm.A randomized complete block design (RCBD) (Federer, 1955) with three replications was used.Each experimental block comprised two parents and all 127 lines of the RIL population.A wide-narrow row planting pattern, with a wide row of 0.48 m and a narrow row of 0.18 m, was used to cultivate all plant materials.Each plot was allocated to two rows with 2 m in length.Thinning was conducted to ensure a density of 25 plants per row once the third leaf emerged.

Phenotype measurement and analysis
PH was dynamically evaluated at three key development stages: jointing (T 1 ), heading (T 2 ), and harvest (T 3 ).At T 1 and T 2 , PH measurement was taken from the soil surface to the apex of the main stem; at T 3 , PH was measured from the soil surface to the tip of the main panicle.For phenotyping, five well-grown plants from the middle of one row, randomly selected in each line, were investigated.Net increases in PH were quantified as DT 1-2 and DT 2-3 , representing the increases from T 1 to T 2 and from T 2 to T 3 , respectively.Generalized heritability (H 2 ) was classified into three levels: low (less than 20%), medium (between 20% and 40%), and high (more than 40%) (Wang et al., 2018).

Re-sequencing
Genomic DNA from leaves of the biparent and RILs was isolated via the CTAB method (Chen and Ronald, 1999).Resequencing was carried out according to the CASAVA 1.8 (Illumina, Inc., San Diego, CA, USA).Genomic DNA was sheared into ~350-bp fragments for library construction.Pairedend reads (150 bp) were sequenced using the HiSeq 2500 system (Illumina, Inc., San Diego, CA, USA), with the target sequencing depth for each sample set at 30×.Raw reads were filtered based on barcode sequences.After trimming, clean reads were aligned to the Yugu1 genome sequence (Setaria italica v2.2) using the Burrows-Wheeler Aligner (BWA) software (Li and Durbin, 2009).Duplicate marking was performed using Picard tools (Picard, http://sourceforge.net/projects/picard/),and Genome Analysis Toolkit (GATK) (McKenna et al., 2010) was used for InDel realignment and base recalibration as part of the preprocessing steps.Subsequently, GATK was utilized again for the detection and filtering of single-nucleotide polymorphisms (SNPs) to obtain the final set of SNP sites.Finally, SNPs identified between two parents were considered polymorphic, with those exhibiting an aa×bb pattern being selected for subsequent analysis.

SNP genotyping and genetic map construction
To enhance the quality of the genetic map, further filtration and selection of polymorphic SNP markers from the initial set of 712,243 SNPs were performed.Sliding scans across chromosomes for genotyping were determined by a window of 15 SNPs and a step size of 1 SNP.When the count of "aa" or "bb" genotypes within the window reached or surpassed 11, it was classified as "aa" or "bb"; otherwise, it was imputed and corrected to "ab" (Huang et al., 2009).In the RIL population, SNP loci were compared against parental genotypes for binning.Samples were organized based on their chromosomal positions, marking genotype transitions as recombination breakpoints and grouping corresponding SNPs into bins.Bins shorter than 10 kb and markers exhibiting significant segregation distortion were further excluded to minimize bias.Filtered bins were utilized for genotyping analysis and further segmented into various linkage groups with HighMap software (Sasaki and International Rice Genome Sequencing, P, 2005).

QTL mapping
PH data from the RIL population were collected from nine different environments (20CZ, 20JZ, 20DT, 21CZ, 21JZ, 21DT, 22CZ, 22JZ, and 22DT).In addition, four best linear unbiased predictions (BLUPs) were calculated to account for environmental variability.Four BLUPs were designated as follows: BLUP1 for CZ, BLUP2 for JZ, BLUP3 for DT, and BLUP4 representing a combined analysis across all nine environments.Each individual test was based on data from a single year, location, or stage, respectively.Unconditional QTL referred to the cumulative effects at T 1 , T 2 , and T 3 stages under various environmental conditions (Fu et al., 2022;Meyer et al., 2023), while conditional QTL referred to the net genetic effects during DT 1-2 and DT 2-3 (Zhu, 1995;Wu et al., 1997).In total, 65 tests were conducted, consisting of 39 tests for unconditional QTL mapping and 26 tests for conditional QTL mapping.Among the unconditional QTL mapping tests, 27 were conducted at T 1 , T 2 , and T 3 stages across nine different environments, and 12 were conducted at T 1 , T 2 , and T 3 stages using four BLUPs.Correspondingly, in the 26 conditional QTL mapping tests, 18 were carried out during DT 1-2 and DT 2-3 across nine environments, and 8 were performed using four BLUPs during DT 1-2 and DT 2-3 .
QTL mapping analysis was conducted using the Inclusive Composite Interval Mapping (ICIM) method of the IciMapping 4.2 software (Meng et al., 2015), with a stepwise distance of 1 cM and a PIN value set at 0.001.Candidate QTL were identified based on a threshold corresponding to a 0.995 confidence level through a permutation test conducted 1,000 times.QTL name was designated as "q+PH+chromosome number+'-'+number" (Mccouch et al., 1997).Unconditional QTL and conditional QTL were distinguished by prefixes "G" and "D", respectively.Major QTL referred to those identified at least in two tests with a LOD score greater than 3.0 and a PVE greater than 10%.Stable QTL were those identified in more than three tests (Sun et al., 2012;Fan et al., 2015).

RNA-seq analysis and prediction of candidate genes
RNA from two parents was prepared from the penultimate internode at T 1 and T 2 , and the internode below the panicle at T 3 , respectively, with three biological replicates.RNA extraction and cDNA library construction were conducted by Biomarker Technologies (Beijing, China) according to standard procedures.The cDNA libraries were sequenced on the Illumina HiSeq 2500 platform with paired-end 150-bp reads.After filtering raw reads with Trimmomatic (Bolger et al., 2014), 116.45 Gb of clean reads were obtained, with the percentage of Q30 bases in each sample being not less than 92.38% (Supplementary Table S1).The clean reads were mapped to the Yugu 1 reference genome using Hisat2 (Kim et al., 2015), and alignments were quantified with StringTie (Pertea et al., 2015).
Gene expression level was calculated using fragments per kilobase of transcript per million fragments mapped reads (FPKM) method.Differential gene expression analysis was conducted using DESeq2 (Love et al., 2014), with criteria set for a false discovery rate (FDR) < 0.01 and a fold change ≥ 2. Following the sequencing of the library preparations on an Illumina platform and the generation of paired-end reads, a comprehensive data analysis was undertaken.This processing included quality control, comparative analysis, functional annotation of genes, SNP calling, quantification of gene expression levels, and differential expression analysis.
To elucidate the biological significance of gene expression changes, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were performed using BMKCloud (www.biocloud.net).KEGG analysis was used to find the pathways of three stages, then DEGs across three stages were identified as candidate genes for dynamic PH.Furthermore, WGCNA was conducted based on PH phenotype and differential gene expression levels to identify significantly related modules, and then the functions of the genes within these modules were then analyzed using BMKCloud.

Expression analysis of candidate genes by qRT-PCR
Total RNA was extracted from the penultimate internode at T 1 and T 2 , and the internode below the panicle at T 3 , using RNAiso Plus (TaKaRa Bio Inc., Shiga, Japan), following the manufacturer's instructions.An equal amount of 2.0 mg of total RNA was reversetranscribed to cDNA using oligo dT primers in the PrimeScript II 1st Strand cDNA Synthesis Kit (TaKaRa Bio Inc., Shiga, Japan).Primer 3.0 (https://primer3.ut.ee/) was employed to design qRT-PCR primers based on sequences from the S. italica v2.2 genome (https://phytozome-next.jgi.doe.gov/info/Sitalica_v2_2).The qRT-PCR reaction was performed on an ABI7500 system in a 10.0-mL volume containing 1.0 mL of template cDNA, 5.0 mL of TB Green Premix Ex Tap II (TaKaRa Bio Inc., Shiga, Japan), 2.0 mL of primer (2 mmol L −1 ), and 2.0 mL of ddH 2 O.The reaction conditions were as follows: 95°C for 3 min, followed by 40 cycles of 95°C for 10 s, 58°C for 30 s, and 72°C for 30 s, followed by a melt curve at 65-95°C by increments of 0.5°C/s.The ACTIN of S. italica was used as the internal reference, and the relative expression levels of interested genes were calculated using the 2 −DDCT method (Kong et al., 2019a, b).All primer sequences used in the present study are listed in Supplementary Table S2.

Statistical analysis
Descriptive statistical analysis for PH of the biparent and RIL population was conducted using IBM SPSS Statistics 17 (SPSS, Chicago, USA).Analysis of variance (ANOVA) and H 2 were processed using the QTL IciMapping 4.2 software (Meng et al., 2015;Ma et al., 2019).BLUPs were estimated using the lme4 package in R software (de Los Campos et al., 2013).The 2 −DDCT of qRT-PCR data were calculated via Microsoft Excel 2013, and mean calculation, significance analysis, and bar graph creation were conducted using GraphPad Prism 9.5.Venn diagrams and heatmaps of differentially expressed genes (DEGs) were drawn using TBtools (Chen et al., 2023).

PH variation
The male parent Jingu 21 consistently exhibited taller PH than the female parent Aininghuang at three stages (T 1 , T 2 , and T 3 ) and during DT 1-2 and DT 2-3 across nine environments almost in all tests.The RIL population showed a continuous distribution across all nine environments, exhibiting relatively small skewness and kurtosis, and two-way transgressive segregation was observed almost at all stages (Supplementary Figure S1).Therefore, the PH of the population followed a normal distribution, making it suitable for QTL analysis (Table 1; Supplementary Figure S1).
Almost under all environments with the exception of Jingu 21 at 22DT, the growth speed of PH in two parents and RIL population rapidly increased during DT 1-2 , with an average growth of 61.99 cm, and then the speed of increase slowed down during DT 2-3 , with an average increase of 27.34 cm (Table 1), indicating that the rapid growth period of PH occurred before T 2 .ANOVA revealed that there were extremely significant differences on PH in the RIL population at T 1 , T 2 , and T 3 stages, attributing to genotype, environment, and genotype × environment interactions (p < 0.01).Furthermore, H 2 was 91.73%, 96.99%, and 97.81% at T 1 , T 2 , and T 3 stages, respectively, indicating that PH exhibited high heritability within the RIL population (Supplementary Table S1).

Genetic map construction
RAD-seq was conducted on two parents and 127 RIL lines, and the RAD-seq data were aligned to the Yugu 1 reference genome.Aininghuang yielded 7.76 Gbp of clean reads with an average coverage of 31×, while Jingu 21 obtained 7.66 Gbp of clean reads with an average coverage of 34×.For 127 RIL lines, 88.21 Gbp of clean reads were obtained, with an average coverage of 2.97×.Finally, a genetic map was constructed using 4,360 bin markers from 712,243 SNPs, spanning 1,016.06cM in length with an average interval of 0.23 cM between adjacent markers (Figure 1A).The longest linkage group, Chr. 9, spanned 144.37 cM and included 597 bin markers with an average genetic distance of 0.24 cM, whereas the shortest group, Chr. 6, spanned 92.48 cM and consisted of 300 bin markers with an average genetic distance of 0.31 cM.Collinearity analysis revealed an average Spearman coefficient up to 0.99 between the genetic and physical maps (Figure 1B; Supplementary Table S3).

Unconditional QTL mapping
A total of 19 unconditional QTL associated with PH, including 7 major QTL, were identified at three stages under nine environments and BLUP data (Table 2).Among those, eight, nine, and eight QTL were identified at T 1 , T 2 , and T 3 , respectively, of which five were repeatedly identified out at least at two stages.These QTL were distributed on chromosomes 1, 4, 5, 7 and 9, respectively, characterized by PVE ranging from 5.23% to 35.36%, and LOD scores spanning from 3.02 to 16.44.Among them, 13 QTL had positive additive alleles contributed by Aininghuang, and the remaining ones were attributed to Jingu 21.Notably, four QTL, namely, GqPH5-1, GqPH5-2, GqPH5-3, and GqPH1-1, were consistently identified at least under three different conditions and were considered stable QTL (Table 2, Figure 2).Specifically, GqPH5-2 was identified in 24 different tests, with PVE ranging from 12.67% to 34.62% and LOD scores ranging from 3.77 to 15.99.This QTL was consistently identified across two or three development stages under four BLUPs, exhibiting substantial PVE (24.91%-34.61%),which further underscores its reliability.In addition, GqPH5-1 also exhibited notable stability across seven tests, with PVE and LOD scores ranging from 20.88% to 35.36% and 7.85 to 16.44, respectively.These results showed that PH was controlled by different QTL at various development stages.

Conditional QTL mapping
During DT 1-2 and DT 2-3 , eight and five QTL were identified, respectively, but no overlapping QTL were identified between them, which explained 6.5%-40.35% of PVE with LOD values ranging from 3.01 to 25.11 (Table 3).The positive additive effects of all QTL were contributed by Aininghuang.In total, five major QTL and three stable QTL were further identified.Notably, DqPH9-1, one of the three stable QTL, was a major co-localized QTL, and expressed only during DT 2-3 , explaining 8.04%-40.35% of the PVE under various environments.This QTL was consistently identified during DT 2-3 under four BLUPs, with a substantial PVE ranging from 12.3% to 40.35%, which further underscores its reliability.Additionally, DqPH5-1 was consistently identified during DT 1-2 under four BLUPs, proving it is also a stable QTL.

Consistent QTL mapping
The comparative analysis between 19 unconditional and 13 conditional QTL showed that a total of 4 QTL, namely, qPH5-1, qPH5-2, qPH5-3, and qPH5-4, were consistent QTL due to their overlapping confidence intervals (Table 4, Figure 2).These QTL exhibited PVE ranging from 12.67% to 35.36% and LOD scores from 3.77 to 16.44.Notably, qPH5-1 and qPH5-2 were highlighted for stable expression across three experiment sites.Moreover, qPH5-2 was identified at multiple development stages, indicating that it played a crucial role throughout the entire developmental process.In addition, qPH5-1 exhibited stable expression across multiple environments and development stages, including that it may be a QTL with a continuous effect.However, qPH5-4 was identified only in DT, suggesting it may possess environmental specificity or adaptability.

RNA-seq analysis
To further identify the key candidate genes underlying PH, we conducted a full-length transcriptome analysis.A total of 131.29 Gb of clean reads, with a GC content of 50.67%-53.00%,were obtained after removing adapter and low-quality reads.The sequencing quality was high enough (Q20 ≥ 97.81%, Q30 ≥ 92.97%) for further analysis (Supplementary Table S4).Upon mapping the quality-controlled sequencing data from two parents to the Yugu 1 reference genome, it was found that over 91.5% of the clean reads corresponded to the reference genome (Supplementary Table S5).Most of the mapped reads were concentrated in the exon region of genes (Figure 3A; Supplementary Figure S2), indicating that the sequencing results were consistent with the characteristics of RNA-Seq.Principal component analysis (PCA) revealed that two distinct groups between two parents were clustered with minimal variation among three biological replicates (Figures 3B, C).
A total of 4,856, 5,715, and 5,759 DEGs were identified at T 1 , T 2 , and T 3 , respectively (Figure 3D).Among these, 239 upregulated genes (Figure 3E) and 475 downregulated genes (Figure 3F) were overlapped and identified across T 1 , T 2 , and T 3 , respectively, The chromosome-wise distribution of QTL for PH.suggesting that a shared group of genes played a continuous and crucial role throughout the entire process of PH development.

GO term and KEGG pathway enrichment analysis
GO analysis of DEGs in four consistent stable QTL regions showed that DEGs were predominantly involved in cellular processes and metabolic processes within the biological process (BP) category, cellular anatomical entity, intracellular and proteincontaining complex in cellular component (CC), and binding and catalytic activity in molecular function (MF) at three stages, respectively (Figures 4A, C, E).Additionally, the KEGG enrichment analysis demonstrated that significant enrichment was present in carbon fixation, photosynthetic organisms and glycosphingolipid biosynthesis-globo isoglobo series at T 1 (Figure 4B), glycerolipid metabolism and biosynthesis of unsaturated fatty acids at T 2 (Figure 4D), and carbon fixation in photosynthetic organisms and fructose and mannose metabolism at T 3 (Figure 4F).
To predict candidate genes, we analyzed pathways that were significantly enriched across three stages.A total of 24 pathways were significantly enriched at T 1 , T 2 , and T 3 stages, including fatty  5A).From these pathways, we identified 61 potential candidate genes for PH within four consistent stable QTL mapping intervals.Notably, five specific genes, namely, Seita.5G363400,Seita.5G372100,Seita.5G394300,Seita.5G402700, and Seita.5G404900,showed significant expression differences between two parents at three stages (Figures 5B, 6A).Compared to Aininghuang, Seita.5G372100,Seita.5G363400, and Seita.5G404900 in Jingu 21 were upregulated at T 1 , downregulated at T 2 , and upregulated again at T 3 ; Seita.5G394300 was upregulated at T 1 , then downregulated at T 2 and T 3 .Seita.5G402700 was downregulated at T 1 , followed by upregulation at both T 2 and T 3 .Those expression patterns suggested that the T 1 stage was a critical period for PH development, and these genes may be involved in key processes such as stem elongation, cell division, and hormone regulation on PH development (Figure 6A).

WGCNA
To explore the gene expression network associated with PH, WGCNA was conducted using a total of 38,494 expressed genes that were detected.These genes were divided into 14 modules, of which 2,   B).Scatter plot analysis revealed that there was a positive correlation between module membership (MM) and gene significance (GS) in the two modules (Figures 7C, E).KEGG enrichment analysis revealed that the lightsteelblue module was significantly enriched in photosynthesis, photosynthesis-antenna proteins, carbon metabolism, biosynthesis of amino acids, and thiamine metabolism pathways, and the salmon2 module was significantly enriched in starch and sucrose metabolism, DNA replication, glycosaminoglycan degradation, protein processing in endoplasmic reticulum, and glycosphingolipid biosynthesis-ganglio series, among others (Figures 7D, F).Based on the KEGG enriched pathways in both modules, four genes, namely, Seita.5G350500,Seita.5G356600,Seita.5G399300, and Seita.5G394300, were identified as candidates for PH.Through the integration of KEGG enrichment analysis of DEGs in RNA-seq with WGCNA, three genes, namely, Seita.5G350500,Seita.5G356600, and Seita.5G399300, were consistently identified, underscoring their potential significance in modulating PH.

Expression analysis of candidate genes by qRT-PCR
To clarify the precise expression level of candidate genes, we randomly selected four candidates, namely, Seita.5G350500,Seita.5G356600,Seita.5G363400, and Seita.5G372100, to carry out expression analysis using qRT-PCR in the two parents (Figure 6B).The results demonstrated that almost all genes were differentially expressed between two parents at three development stages, which were consistent with the results of RNA-Seq (Figure 6A), indicating that these genes were key candidates.

B A
(A) Expression heatmap of eight candidate genes related to PH at three stages.(B) Expression analysis of four candidate genes in biparent via qRT-PCR at three stages.*statistically significant at p < 0.05; **statistically significant at p < 0.01; ns, not significant.

BLUP of quantitative trait
BLUP significantly improved predictive accuracy in both animal and plant populations (Macedo et al., 2020;VanRaden and Cole, 2020).Considering the significant influence of interactions between genes and environment on quantitative traits, phenotype data were typically collected from diverse environments to mitigate environmental impacts and enhance predictive accuracy through phenotype regression analysis of numerous genetic variants (de Los Campos et al., 2013;Ankamah-Yeboah et al., 2020).In this study, BLUP values for three experiment sites (CZ, JZ, and DT) as well as for overall were calculated.The results demonstrated a high degree of agreement between the data and the model, thereby establishing a solid foundation for QTL mapping.

Combination analysis of unconditional and conditional QTL mapping
Traditional QTL mapping provided information on cumulative effects at specific growth stages (Fu et al., 2022;Meyer et al., 2023).However, this approach, which focused on the final values of quantitative traits, neglected the net or incremental genetic effects of QTLs across different development processes (Zhu, 1995;Wu et al., 1997).By integrating unconditional and conditional QTL mapping, it was possible to acquire the genetic loci that affected quantitative traits at various development stages and, further, to clearly elucidate the expression patterns and effects of underlying loci throughout the developmental process (Yang et al., 2006;Cui et al., 2011;Wu et al., 2022).
In the present study, we identified eight, nine, and eight unconditional QTL at T 1 , T 2 , and T 3 stages, and eight and five conditional QTL during DT 1-2 and DT 2-3 , respectively, indicating selective expression at different development stages.Among those, GqPH5-3, GqPH5-5, GqPH6-1, GqPH1-1, and DqPH9-1 overlapped with prior studies (Zhang et al., 2017;He et al., 2021), suggesting that those QTL may exhibit consistent functions or effects across different studies.Remarkably, QTL occurrence was more pronounced in the early stages of PH development, indicating that QTL expression was particularly vigorous during the initial stages, which was consistent with previous studies (Wang et al., 2019;Che et al., 2020;Wu et al., 2022).

Prediction of candidate genes via integration of transcriptome and WGCNA
Transcriptome analysis was effective in identifying genes that responded significantly to specific conditions, and WGCNA focused on identifying gene sets that may not show significant changes in expression individually but work together in biological processes (Yao et al., 2023;Xie et al., 2024).In the present study, we identified 61 DEGs associated with PH from RNA-seq, among which five genes, namely, Seita.5G404900,Seita.5G363400,Seita.5G394300,Seita.5G372100, and Seita.5G402700, were continuously differentially expressed throughout the development stages and considered as candidates.Furthermore, three genes, Seita.5G350500,Seita.5G356600, and Seita.5G399300, were identified to be differentially expressed through WGCNA.In summary, through the integration of transcriptome analysis and WGCNA, we identified eight candidate genes potentially influencing PH.
Notably, the "Green Revolution" gene SD1 in rice, the homology of Seita.5G404900, has been proven to play an essential role in regulating PH (Phillips et al., 1995;Huang et al., 1998;Hedden and Phillips, 2000;Sakamoto et al., 2001;Yamaguchi, 2008;Zhu et al., 2023).Also in rice, OFP2, the homology of Seita.5G363400,decreased PH by interacting with KNOX and BELL genes to suppress gibberellin biosynthesis (Schmitz et al., 2015).RBOHH, the homology of Seita.5G372100,influenced PH by regulating reactive oxygen species (ROS) levels mediated by DELLA proteins in response to both biotic and abiotic stresses in Arabidopsis and rice (Achard et al., 2008;Zhu et al., 2024).Therefore, we randomly selected four candidates, namely, Seita.5G350500,Seita.5G356600,Seita.5G363400, and Seita.5G372100, to carry out expression analysis using qRT-PCR.The results showed that all of them exhibited differential expression between the two parents at three development stages (Figure 6B), which was consistent with the RNA-Seq, indicating that those candidates played a significant role in the development of PH in foxtail millet.

Conclusion
PH is a critical trait influencing lodging, stress resistance, and yield in foxtail millet.In the present study, a total of seven major unconditional QTL and five major conditional QTL for PH were identified using a high-density genetic map with 4,360 bin markers based on a RIL population, of which four QTL were simultaneously identified via unconditional and conditional QTL mapping.Within the four consensus QTL intervals, eight candidate genes were predicted through RNA-seq and WGCNA.This study laid the foundation for fine mapping and cloning of QTL for PH in foxtail millet.
Science Research Project (20210302124233), the Minor Crop Molecular Breeding Platform Special Project of the Shanxi Academy of Agricultural Sciences (YGC2019FZ3), the Central Government Guides Local Science and Technology Development Funds Project (YDZJSX2022A043), the earmarked fund for Modern Agro-industry Technology Research System (2024CYJSTX04-01), and Shanxi Province Science and Technology Innovation Team Project (2015013001-09).
(A) Distribution of polymorphic markers in the genetic map constructed from the RIL population.(B) Collinearity analysis between the genetic map and the physical map.
FIGURE 4 GO term and KEGG pathway enrichment analyses of DEGs obtained by RNA sequencing in four consistent stable QTL regions.(A, C, E) GO terms enriched for three aspects of all DEGs at T 1 (A), T 2 (C), and T 3 (E) stages, respectively, in which MF represents molecular function, BP represents biological processes, and CC represents cellular components.(B, D, F) The top 20 terms for KEGG pathway enrichment analyses of DEGs at T 1 (B), T 2 (D), and T 3 (F) stages, respectively.
FIGURE 5 (A) KEGG pathway enrichment analyses of DEGs at T 1 , T 2 , and T 3 stages in four consistent QTL regions.(B) Expression heatmap of 61 DEGs.
FIGURE 7 Identification of candidate genes via WGCNA.(A) Cluster dendrogram.(B) Module-trait relationships in WGCNA.(C, E) Correlation between MM (Module Membership) and GS (Gene Significance) in the lightsteelblue (C) and salmon2 (E) modules, respectively.(D, F) KEGG enrichment of genes within consistent QTL in the lightsteelblue (D) and salmon2 (F) modules, respectively.GS represents the correlation of each gene within the module, and MM represents the correlation between a single gene and its module.

TABLE 1
PH performance of RIL and their parents.

TABLE 2
Unconditional QTL of PH identified in different environments.

TABLE 3
Conditional QTL of PH identified in different environments.