Longitudinal genome-wide association study reveals early QTL that predict biomass accumulation under cold stress in sorghum

Introduction Sorghum bicolor is a promising cellulosic feedstock crop for bioenergy due to its high biomass yields. However, early growth phases of sorghum are sensitive to cold stress, limiting its planting in temperate environments. Cold adaptability is crucial for cultivating bioenergy and grain sorghum at higher latitudes and elevations, or for extending the growing season. Identifying genes and alleles that enhance biomass accumulation under early cold stress can lead to improved sorghum varieties through breeding or genetic engineering. Methods We conducted image-based phenotyping on 369 accessions from the sorghum Bioenergy Association Panel (BAP) in a controlled environment with early cold treatment. The BAP includes diverse accessions with dense genotyping and varied racial, geographical, and phenotypic backgrounds. Daily, non-destructive imaging allowed temporal analysis of growth-related traits and water use efficiency (WUE). A genome-wide association study (GWAS) was performed to identify genomic intervals and genes associated with cold stress response. Results The GWAS identified transient quantitative trait loci (QTL) strongly associated with growth-related traits, enabling an exploration of the genetic basis of cold stress response at different developmental stages. This analysis of daily growth traits, rather than endpoint traits, revealed early transient QTL predictive of final phenotypes. The study identified both known and novel candidate genes associated with growth-related traits and temporal responses to cold stress. Discussion The identified QTL and candidate genes contribute to understanding the genetic mechanisms underlying sorghum's response to cold stress. These findings can inform breeding and genetic engineering strategies to develop sorghum varieties with improved biomass yields and resilience to cold, facilitating earlier planting, extended growing seasons, and cultivation at higher latitudes and elevations.


Introduction
Sorghum bicolor (L.) Moench, a C 4 crop native to Africa and known for its drought and heat tolerance, is a promising bioenergy crop because of its ability to grow in marginal environments and a high potential for biomass yield (Brenton et al., 2016).Early season planting, which extends the growing season of bioenergy sorghum and potentially increases exposure to rainfall is a strategy to maximize biomass yield (Yu and Tuinstra, 2001;Upadhyaya et al., 2015).Cultivation in northern regions and at higher elevations could increase the land available for bioenergy sorghum production, without using land areas needed for food or feed crop production.
However, early planting and cultivation in temperate areas will not lead to higher biomass yields if plant development stalls.Because of its cold-sensitivity, planting sorghum in temperatures below 12-15°C will diminish the yield (Burow et al., 2011;Salas Fernandez et al., 2014;Chopra et al., 2015).Therefore, identification of sorghum accessions with enhanced tolerance to cold is needed to lengthen the growing season, expand growing regions, and to achieve higher total biomass yields.
Increased cold tolerance is only one of several traits necessary for improved biomass yield.Identification of 'ideotype-positive' and 'ideotype-negative' accessions can prioritize germplasm with multiple positive characteristics that may contribute to increased yield and other desirable traits through breeding or engineering.For example, accessions with high biomass, reduced height, and high WUE may be more desirable for breeding bioenergy sorghum.Alternatively, accessions that are tall but have low biomass and low WUE may be the least desirable for multiple reasons, including a higher propensity for lodging (Esechie et al., 1977).Crosses of lines exhibiting these different ideotypes can also be used to create resources such as nested association mapping (NAM) populations to elucidate the genetic architectures underlying the desired characteristics.
The genetic basis of sorghum's response to cold stress and its potential for cold adaptability must be better understood in order to breed or engineer sorghum lines that can thrive in the lower temperatures of early spring and in colder climates.Cold tolerance is a complex quantitative trait, and there is phenotypic variability for cold tolerance among sorghum accessions (Yu et al., 2004;Chopra et al., 2017).Natural genetic variation in sorghum's response to cold stress has also been identified (Ortiz et al., 2017).Genome-wide association mapping of cold sensitivity traits in diverse germplasm is a promising approach to identify allelic variation that may be harnessed to improve the cold tolerance of sorghum.
The aim of this study was to analyze the phenotypic variability and genetic architecture of bioenergy sorghum's response to early cold stress conditions in a controlled environment.We used a high throughput imaging-based system to collect daily phenotypic measurements from a set of 369 diverse accessions from the sorghum Bioenergy Association Panel (BAP) (Brenton et al., 2016) genotyped with 232,303 high-quality single nucleotide polymorphism (SNP) markers (Supplementary File S1).Daily phenotypic measurements are more beneficial than endpoint measurements because the response to cold stress can change over time and as plants develop.Therefore, regular imaging and phenotyping allowed for the comparison of the cold stress response in the BAP at different developmental stages.
High throughput image-based phenotyping was performed to track growth over time and under early cold stress in this sorghum panel.GWAS on the daily phenotypic data revealed genetic variation underlying the response to early season cold stress and identified transient QTLs related to biomass, height, hull area, water use efficiency (WUE) and relative growth rate (RGR).Longitudinal GWAS was beneficial because it detected early transient QTL that predicted final phenotypes which would not have been possible with GWAS of endpoint data alone.Candidate genes associated with these phenotypes were identified and may prove useful for crop improvement through breeding or targeted genetic modifications.

Population structure and kinship
In this study, we tested 369 accessions of the BAP grown under early cold stress conditions.The BAP (Brenton et al., 2016) comprises six subpopulations that represent a racial, geographic, and phenotypically diverse selection of sorghum accessions.The BAP was designed to limit variation in key bioenergy traits, such as height and sensitivity to photoperiod, within a range of desirable values (Brenton et al., 2016) (Supplementary File S1).We generated a kinship matrix showing the genetic correlations in the BAP (Supplementary File S2).The kinship matrix revealed distinct groups of highly correlated accessions.In many of the highly correlated clusters, the individual accessions shared characteristics such as the country of origin and similar photoperiod sensitivity.For example, in one cluster, 47 of 48 accessions originated from Ethiopia.Another cluster comprised 29 accessions, and 25 of those originated from Ethiopia, while 26 were photoperiod-sensitive.One cluster was comprised entirely of photoperiod-insensitive accessions, while 29 of 30 accessions in another cluster were of the cellulosic type.These clusters reflect an uneven population distribution and illustrate the necessity for controlling for population structure in this study.

Germination and seedling vigor
Seeds were sown in pots at 15°C and grown at that temperature for 31 days.The temperature was then increased to 24°C for seven days and then increased to 32°C for the final 18 days of the experiment.Plant growth and development were observed using daily, nondestructive imaging throughout the eight-week experiment.Germination dates were determined for each plant by image analysis and manual validation.For 369 accessions, at least two plants germinated.Of these 369 accessions, 41 germinated within ten days after planting (DAP) during the 15°C cold treatment, and the remaining accessions germinated between 11 DAP and 46 DAP (Supplementary File S3).Accessions that failed to germinate or died by 46 DAP were excluded, leaving 309 accessions for subsequent phenotypic and genetic analyses (Supplementary File S3).

Trait heritability and phenotypic variation
To analyze the phenotypic response to early cold stress in the BAP, we evaluated five bioenergy-related traits of interest: biomass, height, hull area, water use efficiency (WUE), and relative growth rate (RGR).Images captured daily from three angles and the amount of water each plant received daily enabled calculation of phenotypic measurements for each plant for each day of the experiment.
We used the three replicates for each accession to calculate broad-sense heritability within the experiment for each of the traits at each DAP to determine the proportion of phenotypic variation explained by genetic variation in the BAP.The heritability for each trait changed over time and as the temperature increased (Figure 1A), generally increasing through the early growth period and decreasing as the plants became larger and more root bound.Heritability estimates ranged from 0.544 to 0 (Supplementary File S4).Height at 44-54 DAP was the most heritable trait with a heritability range of 0.544 to 0.503.The heritability values were moderate (0.493 to 0.3) for height at 23-43 DAP, height at 55-56 DAP, WUE at 36-45 DAP, and hull area at 42-56 DAP.The remaining traits had low heritability (0.297-0).
Heritability of WUE increased between 10-20 DAP, remained relatively flat until 42 DAP, then decreased.The decreased heritability of WUE was expected and coincides with an increased uncertainty of biomass calculations and water use estimates during the final two weeks of the experiment.There was a significant increase in the overlapping of leaves after 42 DAP, resulting in the underestimation of plant biomass from the imaging data.The heritability for RGR had no discernable trend.The range of variation observed for the phenotypes of interest is depicted in the density histogram (Figure 1B).The density histogram shows the frequency of the mean of the analyzed phenotypes at the DAP at which each trait had the highest heritability.

Trait correlations
The correlations among biomass, height, hull area, RGR, and WUE were evaluated at 51 DAP (Figure 2).The plants began to overlap and grow outside of the camera view after that date, making area measurements based on pixel numbers less indicative of the actual plant size.Biomass was most highly correlated to WUE with an r 2 of 0.95 but also showed high correlations to hull area and height.The correlation between biomass and RGR was negative, with an r 2 of -0.52.

Phenotypic ranking of accessions
The daily plant images were analyzed to quantify the biomass, hull area, height, and WUE phenotypes for each DAP.Based on the resulting values, for each trait, we ranked the accessions based on their performance.Some accessions were ranked high or low for multiple growth-related traits.The Venn diagram (Figure 3) illustrates overlaps among the top or bottom 5% of accessions at 51 DAP.Supplementary Table S5 provides a complete ranking of accessions for each trait at 51 DAP.Five accessions (PI152771, PI55064, PI514456, PI569454, and PI534165) are represented in the bottom-ranked accessions for all four traits.PI154988, PI535794, PI569148, PI535792, PI329584, and PI535793 are among the lowest-ranked accessions for biomass, hull area, and WUE, while PI569457 is among the lowest-ranked accessions for biomass, hull area, and height.Six additional accessions ranked as poor performers for two traits of interest: PI452544, PI552851, PI570038, and PI570254 for biomass and WUE; PI157035 for hull area and height; and PI540518 for hull area and WUE.Among the top-ranked accessions, PI19770, PI329299, PI455301, Traits of interest show distinct trends in heritability over time and temperature changes.and PI453106 were in the top 5% for biomass, height, hull area, and WUE.PI329403 and PI452619 were among the top accessions for biomass, hull area, and WUE, while PI92270 was a top performer for hull area, height, and WUE and PI563222 was among the top lines for biomass, height, and WUE.Finally, eight additional lines ranked among the top accessions for two traits.

Ideotype-positive and ideotypenegative accessions
Comparisons of accession rankings among multiple traits facilitated the identification of 'ideotype-positive' and 'ideotypenegative' accessions.Accessions ranked within the top 10% for WUE but with a height-to-biomass ratio in the bottom 5% of accessions were designated as ideotype-positive (Figure 4A).There were five ideotype-positive accessions in the early period of the experiment, from 15 DAP to 31 DAP at 15°C: PI329286, PI329517, PI329711, PI455221, and PI562730.There were three ideotypepositive accessions in the late period of the experiment, from 39 DAP to 56 DAP after the temperature had been increased to 32°C: PI329471, PI329517, and PI585461.Accession PI329517 was ideotype-positive throughout the experiment.

Temporal WUE profiles and clustering of accessions
The daily quantification of traits provided the opportunity to assess temporal changes in trait values as the temperature increased and plants developed.For brevity, we only refer to the five top and bottom-ranked accessions for WUE, a key trait.Supplementary Table S6 provides the complete ranked list of accessions for each trait.Figure 5 shows the temporal WUE profiles for the five highestand lowest-ranked accessions at each temperature stage.The topranked accessions for WUE at each temperature exhibited a WUE profile over time distinct from the profiles of the accessions ranked lowest for WUE accessions, and the separation became more apparent as the plants developed and the temperature increased (Figure 5).Even at 15°C, the highest-ranked accessions for WUE are distinguishable with a steeper increase in WUE over time compared to the bottom-ranked accessions.Some accessions were ranked high or low for WUE at one temperature stage only, whereas others were Overlaps of high and low performing accessions at 51 DAP.Venn diagrams depict numbers of accessions represented in the lowest-performing accessions (A) or highest-performing accessions (B) that overlap for biomass, WUE, hull area, or height at 51 DAP.Ideotype-positive and ideotype-negative accessions for WUE and biomass-associated traits.(A) Early and late ideotype-positive radar plots show accessions within both the top 10% of WUE and lowest 5% of height:biomass ratios.Each dashed square represents a percentile range: 0% is the center-most square, and 100% is the outermost square.Accessions whose phenotypes fall within a given percentile range for a given trait (height, height:biomass ratio, or WUE) reside on the inside of the corresponding square.(B) Early and late ideotype-negative radar plots show accessions within both the lowest 10% of WUE and the highest 5% of height:biomass ratios.(C, D) Heatmaps show the change in rank of early and late ideotype-positive and negative accessions over time and at each temperature stage for WUE (C) and height:biomass ratios (D).Dark blue indicates the lowest WUE (an ideotype-negative trait) and lowest height:biomass ratio (an ideotype-positive trait).Dark orange indicates the highest WUE (an ideotype-positive trait) and highest height:biomass ratio (an ideotype-negative trait).
ranked high or low for WUE across temperatures.For example, PI550604 was among the bottom-ranked accessions for WUE at each temperature, and PI534165 was among the bottom-ranked at 24°C and 32°C.PI455221 was a top-ranked accession for WUE at both 15°C and 24°C, while PI329403 was among the highestranked accessions at 15°C and 32°C.
For each accession, we analyzed the profile of each phenotype over time and at different temperatures.Again, for brevity, here we only refer to WUE.The temporal WUE profiles for the BAP accessions separate into seven clusters, each represented by a growth curve represented by the average of all accessions in that cluster (Figure 6A).The temporal WUE profile for each of the seven groups is distinct.The mean WUE for accessions that comprise Cluster 1 is low at 15°C, shows a minimal increase at 24°C, and increases steadily at 32°C, but never reaches levels as high as several other clusters.Cluster 2 shows a slightly higher WUE at 15°C compared to Cluster 1 and an increase in WUE at 24°C and 32°C, but never exceeds moderate levels of WUE.In Cluster 3 there is little increase in WUE until the temperature increases to 32°C and this cluster eventually reaches moderate WUE.WUE for Cluster 4 is low at 15°C, increases at 24°C and 32°C.Cluster 5 exhibits a modest increase in WUE at 15°C, a rapid increase at 24°C, but WUE increases slowly at 32°C.Cluster 6 has an early and robust increase in WUE across all three temperatures in the experiment and reaches the maximum peak WUE of all seven clusters.Finally, WUE for Cluster 7 is moderate at 15°C and increases steadily at 24°C and 32°C , but ultimately peaks at a moderate level.Cluster analysis for each additional trait is in Supplementary File S7.The dendrogram tree (Figure 6B) depicts the hierarchical clustering of the BAP accessions into seven groups based on their temporal WUE profiles, with each cluster represented by a different color.The clustergram of the PCA-weighted mean (Figure 6C) illustrates the divergence of the temporal WUE profiles of the BAP accessions into the seven stable clusters which captured the maximum phenotypic difference for this trait.

Longitudinal genome-wide association study
A longitudinal genome-wide association study (GWAS) was performed to characterize the genetic architecture underlying response to early cold stress conditions in the BAP.GWAS was performed using the multi-locus mixed model (MLMM) algorithm (Segura et al., 2012).MLMM uses multiple loci in the model to yield a higher detection power and lower the potential of false discoveries.Figure 7 illustrates the detection of highly significant signals related to each trait of interest on all ten sorghum chromosomes.Each significant SNP has a maximum corrected p-value of 0.05.This analysis identified 2,305 highly significant SNPs associated with biomass, height, hull area, WUE, and RGR over the course of the 56-day experiment.
To identify SNPs significantly correlated with phenotypes at a specific DAP and correlated with temperature changes, we also used the values for each phenotype at each DAP (Figure 8) for the GWAS.
This analysis revealed several QTL that displayed a transient response, appearing at specific times and developmental stages, and as temperatures increased.(Figure 8).For example, the highly significant hull area SNP 7:2,934,702 was transiently detected from 27 to 30 DAP and again from 34 to 35 DAP (Figure 9).This SNP was also detected at 17, 18, 26, 31, and 33 DAP using a lower confidence model (Figure 10).
For the transient QTL that were detected, we examined the phenotypic differences for the alternative alleles at a particular SNP.For the hull area SNP 7:2,934,702 shown in Figure 9, allele boxplots (Figure 10A; Supplementary File S8) shows that a "G" at the position corresponds to a significantly higher value for hull area than a "C" for 27 to 30 DAP and 34 to 35 DAP. Figure 10B shows that a "G" at this position also contributes to higher hull area through the final day of the experiment, essentially predicting the endpoint phenotype from detection of the earlier transient QTL.SNP 7:2,934,702, within gene Sobic.007G033300,encodes a plastocyanin-like protein.Plastocyanins are copper-containing chloroplast-localized proteins involved in photosynthesis and function in electron transfer between photosystem II and photosystem I (Sigfridsson, 1998).This putative plastocyanin, Sobic.007G033300, is highly expressed in sorghum shoots, leaves, and early inflorescences (https://phytozome.jgi.doe.gov/pz/portal.html)(Makita et al., 2014).We examined other sequence variation existing within this gene using public re-sequencing data and found seven SNP variants resulting in nonsynonymous (NS) amino acid substitutions and an additional inframe three basepair deletion segregating within individuals in the BAP (Supplementary File S9A).Eighty of the phenotyped BAP accessions contain no variants in Sobic.007G033300relative to the reference sequence, while the remaining accessions contain between one to eight of the variants summarized in Supplementary File S9A.When grouped, individuals carrying at least one variant in Sobic.007G033300have a significantly larger hull area from 27 DAP through the end of the experiment compared to accessions containing no variants (Supplementary File S9B).
The GWAS also identified pleiotropic QTL.A transient QTL on chromosome two (SNP 2:1,199,928) was significant for WUE and biomass from 14 to 19 DAP.SNP 4:64256235 on chromosome four was a significant transient QTL for WUE and biomass from 14 to 19 DAP, and height at 34 DAP.Transient QTLs on chromosomes six and nine (SNP 6:42,116,590 and SNP 9:55067080, respectively) were significant for WUE and biomass at 35 to 36 DAP.S10 shows the traits and DAP associated with each SNP.

FIGURE 8
Temporal manhattan plot of longitudinal GWAS on biomass, growth, and WUE traits reveals transient QTL.The plot depicts highly significant SNPs corresponding to the traits of interest at specific DAP.The y-axis denotes the DAP, and the x-axis represents the sorghum chromosome.Black boxes highlight transient QTL that "turn on/off" with developmental stage and as the temperature increases.Manhattan plot of GWAS on biomass, growth, and WUE traits.The Manhattan plot depicts significant SNPs for each trait of interest and the corresponding -log 10 (p) of the GWAS mixed linear model with no other SNPs included in the model.

Candidate gene identification
We conducted a genome scan of regions 15 kb upstream and 15 kb downstream of each significant SNP to identify candidate genes.Functional annotations, putative homologs, and polymorphisms within the candidate genes based on public genomic re-sequencing data were analyzed using Phytozyme (https://phytozome.jgi.doe.gov/pz/portal.html).GWAS identified highly significant QTL near 72 candidate genes (Supplementary File S11) with putative functions potentially related to biomass, height, hull area, WUE, and RGR, and in response to cold.Notably, many of these candidate genes are situated near pleiotropic QTLs, suggesting their involvement in cold stress response alongside biomass-associated phenotypes.
Among these, we spotlight an a priori candidate gene, Sobic.006G057866,encoding PSEUDO-RESPONSE REGULATOR 7 (PRR37)/MATURITY 1 (Ma1), which plays a crucial role in the regulation of flowering time in sorghum (Murphy et al., 2011).A specific SNP within this gene, SNP 6:40,312,463, significantly associated with RGR at 33 DAP, results in a non-synonymous amino acid substitution (Asn184Lys) in Ma1.Although accessions carrying this polymorphism show a slightly reduced mean RGR at 33 DAP, they were not significantly different (Supplementary File S12A).This observation suggests that allelic variation at SNP 6:40,312,463 is not likely to be a robust contributor to phenotypic variation in RGR.However, interestingly, most accessions (~90%) including both BAP accessions and lines from other diversity panels carry a polymorphism in Ma1 resulting in loss of the annotated reference stop codon, which could drastically increase the protein size and affect its function (Supplementary File S12B).This observation underscores the complexity of genetic factors influencing plant growth and development, highlighting the importance of considering a wide array of genetic variations in understanding the phenotypic outcomes of key agronomic traits.

Discussion
Sorghum is attractive as a bioenergy feedstock crop because it is heat and drought tolerant and can thrive in marginal environments.It is an ideal target for accelerated improvement through breeding and engineering because it has extensive genetic and phenotypic diversity but has not yet benefited from genomics or genetic

A B
Phenotypic differences in hull area associated with a SNP in a plastocyanin-like gene.(A) The allele boxplot shows the correlation of a "G" at SNP 7:2,934,702 with a significantly larger hull area at 27 DAP.(B) Temporal allele plots show the correlation of a "G" at SNP 7:2,934,702 maintains larger hull area than "C" at each DAP.Two asterisks represents SNP identification with the stringent GWAS models while one asterisk represents identification with the lower confidence model.modification like some other crops such as maize.Early season planting of sorghum provides the opportunity for an extended growing season with higher potential accumulation of biomass.However, as a tropical crop, sorghum is sensitive to cold stress.Identification of accessions that exhibit the highest WUE, RGR, height, hull area, and biomass under early cold stress conditions could facilitate genetic improvement of bioenergy sorghum for early season planting and cultivation in colder temperatures.
Our results identified the top performing BAP accessions for bioenergy-related traits under early cold stress.We identified the accessions with the highest and lowest rankings for each trait and multiple bioenergy-related traits.Accession PI329299 was among the top accessions for hull area and height, PI452619 was a top accession for hull area and WUE, PI329403 was a top accession for hull area and WUE, and PI585461 was among the top accessions for both biomass and WUE.These top-performing accessions, particularly the accessions that possess multiple advantageous traits, are promising candidates for sorghum bioenergy breeding programs and the development of additional genetic resources such as mapping populations (e.g., NAMs).
We defined ideotype-positive and ideotype-negative accessions as accessions that ranked high or low, respectively, for multiple traits or exhibited beneficial trait combinations.The ideotype-positive phenotype is the combination of high WUE with a low height-tobiomass ratio, i.e., a plant that achieves high biomass but does not grow tall.This ideotype is desirable for bioenergy sorghum because the required biomass accumulation is attained with reduced water use and without the increased risk of lodging associated with taller plants.Five accessions were ideotype-positive in the early phase of the experiment when the temperature was 15°C.These accessions may be beneficial for breeding bioenergy sorghum that can be planted early in the season or grown in colder environments.Accession PI329517, which attained the ideotype-positive phenotype early under cold stress, and maintained it after the temperature increased, may be beneficial for early season planting, for long growth periods, or under conditions of reduced water availability.In contrast, ideotype-negative accessions that have relatively low biomass accumulation, ranked low for WUE, and are tall may be undesirable for bioenergy sorghum because more water would be required and with a higher risk of lodging and decreased biomass accumulation.Even though these accessions may not be desirable for breeding, identification of these extreme ideotypes may be useful for the development of structured populations for further genetic analysis (Donald, 1968;Rötter et al., 2015).
The profiles of the traits over time and in response to temperature changes could be discerned, ranked, and clustered using the phenotypic data collected with single-day resolution.The BAP clustered into 6 to 8 distinct temporal profiles for each trait, identifying accessions that perform best under early cold stress, and at different temperatures and development stages.For example, accessions in WUE clusters two, six, and seven have the highest WUE at 15°C.Similarly, cluster six not only attains high WUE under 15°C but is the cluster with the highest WUE attained under increased temperatures and therefore may perform well under drought conditions.Accession PI455221 was among the top-performing accessions for WUE at both 15°C and 24°C, and PI329403 was among the top-performing accessions for WUE at 15°C and 32°C, demonstrating that high WUE can be maintained as the plants develop and as the temperature increases.
GWAS revealed 2,305 highly significant SNPs associated with biomass, hull area, WUE, RGR, and height phenotypes.We also identified significant SNPs that colocalized for multiple traits, suggesting that these polymorphisms may be near tightly linked genes or genes with pleiotropic effects.By using the daily values for each phenotype in the GWAS, we determined the temporal correlations of SNPs to traits.With this approach, we identified transient QTL where significance will "turn on/off" at a specific developmental stage or with a change in temperature.
Our GWAS analysis identified 72 candidate genes with functions pertinent to bioenergy production and cold-stress response.Notably, we discovered a transient QTL for Relative Growth Rate (RGR) at 33 days after planting (DAP) that maps to the gene Sobic.006G057866.This gene encodes PRR37/Ma1, a known repressor of flowering in long days (Murphy et al., 2011), Notably, the "activation" or significance of this QTL at 33 DAPimmediately following a temperature increase from 15°C and during a pivotal phase for sorghum's shift to floweringunderscores the significance of PRR37/Ma1.This observation is consistent with this gene playing a crucial role in promoting growth and biomass accumulation by regulating the plant's transition from vegetative growth to reproductive development.
Additionally, a separate transient QTL affecting hull area showed significance between 27 to 30 DAP and again from 34 to 35 DAP, mapping to Sobic.007G033300.This gene encodes a putative plastocyanin, implicated in photosynthetic electron transfer, with sequence variation correlating with increased hull area during these periods.The consistent observation of higher hull area values associated with the "G" allele at this locus emphasizes its significant role in shaping the final phenotype.It also distinctly illustrates that, had GWAS analysis been limited solely to endpoint phenotypes, such a QTL might have remained undetected, thereby obscuring its predictive value for the final phenotype.This case strongly supports the utility of conducting longitudinal GWAS to identify early-season genetic markers of vital traits like yield and biomass.Integrating this approach, which emphasizes the temporal dynamics of trait development, we can see how genes like Sobic.007G033300 can be fundamental in biomass accumulation, consistent with known deleterious effects of cold on photosynthesis, in particular, the capacity for electron transport and photosystem protection (Liu et al., 2018).The significance of plastocyanins, also supported by studies in other crops such as rice and Camelina (Zhang et al., 2017;Okooboh et al., 2023), underscores the complex genetic underpinnings of traits essential for yield and stress resilience.
In summary, we tracked growth over time and under early season cold temperature stress in a genetically diverse sorghum population using daily image-based phenotyping.We analyzed the changes in traits over time, which would not have been possible with the collection of endpoint phenotypes alone.The BAP segregated into distinct clusters for each trait, identifying accessions that performed best under early cold stress, and at different temperatures and development stages.GWAS on the daily phenotypes revealed transient QTLs for highly heritable bioenergy-relevant traits that could predict final phenotypes.In some cases, we could identify putative causative alleles in candidate genes that correlated to a significant difference in phenotype.These findings may facilitate targeted genetic modifications or genomicsdriven breeding efforts for the improvement of bioenergy sorghum.The transient nature of significance for certain QTLs, particularly those that reach significance under specific environmental conditions, underscores the importance of developing breeding strategies that consider the entire growth cycle.This approach enables the selection of sorghum lines that not only start strong but also maintain desirable traits throughout their development, optimizing for both yield and resource efficiency.

Plant materials and genetic data
This study used the sorghum Bioenergy Association Panel (BAP) (Brenton et al., 2016).Details about the panel design, GBS genotyping, marker distribution, population structure, and linkage disequilibrium (LD) decay have been previously described (Brenton et al., 2016).We selected 369 BAP accessions genotyped at 232,303 SNPs.The BAP accessions represent a racially, geographic, and phenotypically diverse selection of sorghum accessions, but is limited to accessions exhibiting key bioenergy traits, such as height, sensitivity to photoperiod, and delayed flowering.

Experimental conditions
Three replicates of 369 BAP accessions (1131 plants) were planted in 600 g of Turface ® in 8-inch tall tree pots, with +14-14-14 Osmocote (1.5 lb/cubic yard) fertilizer.The potted seeds were held overnight in a growth chamber at 32°C (day)/22°C (night).The next day, the pots were loaded into carriers on an automated phenotyping system within a controlled-environment plant growth chamber (Fahlgren et al., 2015).The phenotyping system moved the plants on a closed-loop conveyor path to stations for daily watering, weighing, and imaging.The plants were positioned in a randomized block design and were rotated one-half lane each day to reduce edge effects.To study the impact of cold stress on these accessions, plants were grown at 15°C (day)/15°C (night) for 31 days, 24°C (day)/19°C (night) for 7 days, and 32°C (day)/22°C (night) for 18 days.The soil was maintained at 100% field water capacity by watering the plants twice daily to a target mass of 1192 grams.The target mass was calculated by adding the mass of the plant carrier (342g), the water mass at saturation (250g), and the mass of the Turface ® -filled pot (600g).Water was added after each potted plant in its carrier was weighed, and the volume of water added to reach the target mass of 1192 grams was recorded.In our irrigation protocol, the calculation of the target mass for watering did not account for the increase in plant biomass over time.This decision was based on the practical challenges associated with nondestructive biomass measurements and their potential impact on the statistical integrity of our study.Although this methodology meant that the gradual increase in biomass was not directly factored into the water added, the target mass of 1192 grams was established to ensure that the soil moisture levels were optimally maintained for plant growth throughout the experiment.

Image collection and phenotyping
Each plant passed through a visible light imaging chamber daily while on the automated phenotyping system.The imaging cameras recorded two side views (sv) and one top view (tv) image.As the plants grew, the fields of view on the cameras were adjusted so that the entire plant could be captured in each image.The optical zoom level was reduced for the top view and side view images at 19 DAP and 40 DAP.Scaling factors were calculated for both area and height using a reference object of known size, so that pixel areas across zoom levels were comparable.After eight weeks (56 DAP), the plants were removed from the phenotyping system.The shoot of each plant was cut at the base of the stem.Fresh weight measurements of the shoots were collected immediately.
The images were analyzed using Plant Computer Vision (PlantCV), an image-processing tool coded in Python (Fahlgren et al., 2015).The pixel areas from the daily top view and two side view images of each plant were analyzed with PlantCV to generate measurements of the area, hull area, height, RGR, and WUE.
The area was calculated by adding the pixel count from the top view image to the two side view images.Endpoint fresh weight measurements collected at 56 DAP were correlated with area calculations at 51 DAP to estimate biomass (Fahlgren et al., 2015) (Supplementary File S13).The values at 51 DAP were used for the correlation because, after that date, plants began to overlap and grow outside of the camera's field of view, and the resulting pixel counts were less indicative of the actual plant size.
The hull area is the convex hull calculated from the pixel count in the smallest area that includes a set of given points in a plane.WUE was calculated by dividing the derived area by the cumulative water added to each plant.WUE was calculated by dividing the derived area by the cumulative water added to each plant.Here, it is crucial to recognize the limitations of our methodology, which reflect broader challenges in the precise measurement of complex traits such as WUE.Specifically, our estimation of WUE relied on non-invasive imaging techniques for biomass measurement, without direct assessments of soil moisture content, transpiration rates, or an analysis of root biomass.Height was determined from the side view images The relative growth rate (RGR) was calculated as described in Hoffmann and Poorter, using estimator 2 to determine the RGR for each distinct time point (Hoffmann and Poorter, 2002).Briefly, we calculated the natural log of all replicate area values in the experiment, then calculated the mean of the log-transformed values for two time points, t 1 and t 2 .The mean value of the logtransformed areas for t 2 , W2, was subtracted from the mean value of the log-transformed areas for t 1 , W 1 , and then divided by the difference between t 2 and t 1 as shown: Germination dates were determined for each plant based on the top view area measurements.The earliest DAP with a top view area measure was considered the germination date.The resulting germination date for the three replicates of each accession was averaged to determine the germination date for each accession.

Data processing and analysis
Image-derived phenotypic data were generated for 309 of the BAP accessions (Supplementary File S1).Phenotypic analysis was performed on the image-derived data using the R statistical software (R).Plants that never germinated or that died by 46 DAP were excluded from the analysis.A conservative initial outlier removal step was performed on each phenotype to remove likely artifacts from the image analysis.A data point was considered an outlier if its value was greater than 40 median absolute deviations (MAD) from the median (Davies and Gather, 1993).Raw measurements from PlantCV were smoothed using predicted values from a loess smoothing fit.These fitted values were used as the phenotypes for further analyses (Feldman et al., 2017).
The 'lmer' function in the lme4 R package was used to estimate variance components for broad-sense heritability (Bates et al., 2015).Heritability was calculated based on the subset of accessions that had three replicates that germinated.The model used to calculate variance components was as follows: where Y i is the phenotypic observation; m is the grand mean; genotype i is the effect of the ith genotype; and e i is the random error term.All terms were fit as random effects.Using the variance estimates from the model, broad-sense heritability was calculated as: After the heritability calculation, the median value of the replicates for each accession was used for further analysis.Traits were tested for normality and transformed as necessary using the Box-Cox procedure as implemented in R with the 'boxcox' function in the MASS package (Venables and Ripley, 2002).

Hierarchical clustering
Hierarchical clustering of each trait from 12 to 56 DAP was performed using the 'hclust' R function using a Euclidean distance matrix and the Ward agglomeration method (Murtagh and Legendre, 2014).Clustering results were visualized both as a dendrogram using the R package 'dendextend' and as a clustergram using the function 'clustergram.R' (Schonlau, 2002;Galili, 2010).

Ideotype selection
Height-to-biomass ratios were calculated for each day and converted to a percentile with values between 0 and 100.The average percentile of the accessions for height-to-biomass ratio and WUE over the growing period, DAP 15 to 31 or DAP 39 to 56, for the early and late periods, respectively, was used to select ideotype accessions.Ideotype-positive accessions were defined as those with a height-to-biomass ratio in the bottom 5 percent of the population and a WUE in the top 10 percent of the population.Ideotype-negative accessions have a height-to-biomass ratio in the top 5 percent of the population and a WUE in the bottom 10 percent.Radar plots of the ideotype accessions were generated with the R package 'fmsb' (Nakazawa and Nakazawa, 2018).Heatmaps were generated using the R package 'pheatmap' (Kolde and Kolde, 2018).

Genome-wide association study
Genotyping-by-sequencing (GBS) SNP markers for the BAP have been previously described (Brenton et al., 2016).GWAS was performed using a multi-locus mixed linear model (MLMM) in R to identify loci associated with each trait of interest.The first three principal components of the genotype matrix were included as covariates in the mixed model to control for population structure (Price et al., 2006).A kinship matrix, calculated from the genotype matrix using the Astle-Balding method in the 'synbreed' package, was included as a random effect to control for familial and cryptic relatedness between accessions (Astle and Balding, 2009;Wimmer et al., 2012) (Supplementary File S2).
MLMM tests for association with the phenotype using a stepwise mixed model regression.In each step, the SNP with the most significant association to the phenotype is added to the model as a covariate.Stepwise addition of SNPs continues until the heritable variance estimate (pseudo-heritability) reaches 0 (Segura et al., 2012).A final set of high-confidence SNPs were selected as those that were included as covariates in either of the two optimal models chosen by the MLMM software: the multiple-Bonferonni model or the extended BIC model (Segura et al., 2012).In the multiple-Bonferroni model, all cofactors with a p-value below a Bonferroni corrected threshold are selected.Multiple-Bonferroni was the more stringent of the two models (Segura et al., 2012).The extended BIC model selects a model based upon BIC penalized by the model complexity (Chen and Chen, 2008;Brenton et al., 2016).
The temporal GWAS results were analyzed using ZBrowse, an interactive browser that runs using R (Ziegler et al., 2015).Using ZBrowse, we were able to view the SNPs for multiple traits simultaneously and plot those traits over time to determine the DAP on which each SNP was significant and thus identify transient QTL peaks that turn "on/off" at specific DAP.
The BTx623 sorghum reference genome version 3.1 was used to identify genes colocalizing with or adjacent to the associated SNPs.A genome scan of 15 kb upstream and downstream of each significant SNP was performed to identify candidate genes.Phytozome (https://phytozome.jgi.doe.gov/pz/portal.html)was used to analyze the functional annotation of candidate genes and to identify putative homologs in other species.Analysis of polymorphisms within candidate genes was done using three diversity panels of re-sequencing data available for sorghum via Phytozome (Paterson et al., 2009;Mace et al., 2013;McCormick et al., 2018).VCFtools generated final variant calls after merging the three diversity VCFs, based on a minor allele frequency cutoff of >0.05, and to create distinct VCFs for the accessions used in this study (Danecek et al., 2011).Variant effects were estimated using the SNPeff pipeline and the v3.1 sorghum reference to assess the potential impact of sequence variants on annotated gene models (Cingolani et al., 2012).
(A) Heritability of traits over time at each temperature (15°C, 24°C, and 32°C).(B) Density histograms of each trait analyzed show the day at which each corresponding trait is most heritable (height at 48 DAP, WUE at 42 DAP, hull area at 50 DAP, biomass at 50 DAP, and RGR at 24 DAP).Each trait was transformed into an approximately normal distribution using the Box-Cox method.

FIGURE 2
FIGURE 2Correlations among biomass, height, hull area, RGR, and WUE at 51 DAP.Dark blue represents the strongest correlation, and dark orange represents the weakest correlation; the size of the circles are proportional to the correlation coefficients.The r 2 for each correlation is in the corresponding box below the colored circles.
FIGURE 6 BAP accessions separate into seven distinct temporal WUE clusters.(A) The temporal WUE profiles of BAP accessions from 12 DAP to 56 DAP separate into seven clusters.The growth curves represent the average WUE profile over time for each cluster.(B) Hierarchical dendrogram tree showing the segregation of the BAP accessions into 7 clusters, each represented with a different color.(C) Clustergram of the PCA-weighted mean of each cluster showing the divergence of the BAP accessions into the seven clusters, and the distances between the clusters.The cluster numbers and colors correspond to the dendrogram (B) and temporal WUE profiles (A).

FIGURE 9 A
FIGURE 9A transient QTL associated with the hull area trait.The boxed region represents a transient hull area QTL at SNP 7:2,934,702 that "turns on/off" with developmental stage and as temperature increases.There are five genes within 15kb of the transient QTL.SNP 7:2,934,702 is directly on candidate gene Sobic.007G033300.