A S. bicolor diversity set consisting of n = 330 inbred lines of different origin, type of use (grain, dual-purpose, and forage) and subspecies, comprising all five sorghum races (bicolor, caudatum, durra, guinea, and kafir) and most of their respective intermediaries, was utilized for the present study. All these lines are photo-insensitive and show a similar, early maturity, allowing for an adequate scoring of seed set traits under cool environments with a short growing season. A representative selection of early-maturing sorghum conversion lines, obtained in the 1960s by repeated backcrossing of genetically diverse tropical accessions to a short, photoperiod-insensitive cultivar (Stephens et al., 1967) accounts for 35% (n = 117) of the set. A further 22% (n = 73) consist of other publicly available accessions from temperate countries, mainly China, Russia, the United States, Turkey, and Hungary. The original germplasm of these accessions and the conversion lines was received from the United States Department of Agriculture Agricultural Research Service (USDA-ARS). The remainder of the diversity set consists of 140 diverse breeding lines (grain and dual-purpose sorghum) from a joint breeding program of Justus Liebig University Giessen, Norddeutsche Pflanzenzucht Hans-Georg Lembke KG (Hohenlieth, Germany) and Deutsche Saatveredelung AG (Lippstadt, Germany). The composition of the diversity set is shown in Supplementary Table 1.

Field trials were conducted at five locations (three in Germany and two in Mexico, Table 1). At Texcoco (TEX17, TEX18), one of the Mexican locations, the diversity set was scored in two subsequent years (2017 and 2018), while at all other locations, the experiments were carried out in 1 year (2017 or 2018). Hence, data from six environments (location × year combination) are available for analyses. The locations of this study represent different mega-environments contrasting in climate, radiation and day length. Among the German locations, Asendorf (AS; located in NW-Germany) and Poel (PL; a small island in the Baltic Sea) have a cool, maritime climate, usually providing harsh conditions for sorghum. In contrast, Gross Gerau (GG) is located in the Upper Rhine Valley and characterized by a warm and sunny climate, being a suitable control environment without cold stress. The locations in Mexico differ strongly from Germany, having shorter days during the growing season but much stronger radiation. While San Juan del Río (SJR; 1,920 m, federal state Querétaro) is considered to be at the altitude limit for commercial sorghum cultivation in Mexico, Texcoco (2,250 m, federal state México) is a tropical high-altitude stress environment for sorghum, providing the lowest minimum temperatures of all locations.

To avoid shading effects by neighbors of different plant height, the inbred lines were split into three subgroups, based on previous scorings of plant height. These subgroups were planted in adjacent but separate blocks at all sites. Within these subgroups, an un-replicated randomized complete block design was used. Entries were grown in microplots, consisting of single rows (2.5 × 0.7 m) at Gross Gerau and double rows (2.5 × 1.4 m) at all other locations, with 0.7 m row spacing and a plant density of approx. 20 plants/m². Plant protection and fertilizer application were executed according to site specific best practice. Per entry, the primary panicles of five healthy plants were covered before anthesis with a transparent Cryovac^® bag (330 mm × 750 mm, 15 μm) to avoid cross pollination. These five self-pollinated panicles were considered as biological replications for further analyses. However, due to lodging and insufficient maturity, some plants and entries had to be excluded site-specifically. At maturity, the panicles were harvested with secateurs and dried. The peduncles of each panicle were cut just below the first branches before determining the panicle dry weight. Seed yield (SY) per panicle was measured after threshing. Subsequently, PHI was calculated according to Krishnamurthy et al. (2014) as follows:

Consequently, a PHI value of 0 implies absolutely no seed set, while values close to 1 indicate a high seed set. However, even assuming complete spikelet fertility, PHI will be <1 due to the panicle raw weight. Moreover, seed number (SN) was measured using a Contador seed-counter (Pfeuffer, Germany). In addition to these seed set traits, plant height and start of flowering (in all environments except SJR) were scored on a plot basis.

For statistical analyses of SY, seed number (SN) and PHI, a general linear model was used (IBM SPSS Statistics version 27, IBM Software, Armonk, NY, United States), in which entries (genotypes) and environments (combination of location and year) were considered as fixed effects and replicates (individual plants) as random effects:

where μ represents the population mean, G_i is the genotypic effect, E_j is the environmental effect, GE_ij is the genotype-by-environment interaction, R_kj is the replicate effect, and e is the residual effect.

To compare the levels of genotypic variance obtained in the different environments, ANOVA was also computed separately for each environment, using the following general linear model, where genotypes were considered as fixed and replicates as random effects:

To determine which environments were significantly different from one another, Student-Newman-Keuls post hoc-test was applied.

The heritability was calculated as proposed by Piepho and Möhring (2007) using the following formula:

where h² represents broad-sense heritability, σ²_G is the genotypic variance calculated by a random effect model considering genotype and environment as random factors, and v⁢d¯ is the average variance of the difference between two means.

The phenotypic stability of the best performing inbred lines with a complete data set in all stress environments was analyzed using the coefficient of variation (CV) across the mean trait values of a line in each environment.

The sorghum diversity set was genotyped using DArTseq reduced-representation sequencing¹ to identify genome-wide single-nucleotide polymorphism (SNP) markers at high resolution. SNPs with more than 20% missing data or a minor allele frequency lower than 5% were removed from the final dataset. After filtering, a total of 21,520 high-quality SNPs remained and were used for downstream analyses.

Analysis of phylogenetic relatedness was conducted with TASSEL version 5.0 (Bradbury et al., 2007) using the neighbor-joining method (Saitou and Nei, 1987). Dendroscope 3.7.3 (Huson et al., 2007) was utilized to visualize the genetic relatedness among accessions in a phylogram. The R package GenABEL (Aulchenko et al., 2007) was used to perform a GWAS for the target traits. To adjust population stratification, a mixed linear model approach was implemented by using the kinship matrix as covariates (Stich et al., 2008). We used a threshold of −log10(p) ≥ 3 to minimize type II error and identify SNP-trait association (Gabur et al., 2020).

Candidate genes were identified based on the S. bicolor reference genome v3.1.1 hosted by Phytozome 12², the same reference genome used by Diversity Array Technology to call SNPs. Linkage disequilibrium and haplotype blocks across the entire genome were calculated using the squared allele frequency correlations (r2) between each pair of SNPs, using an LD threshold of r² > 0.7 block as described by Gabriel et al. (2002) and implemented using the tool LDBlockShow (Dong et al., 2021).

Genes within haploblocks surrounding trait-associated SNPs with the maximum −log10(p) values from GWAS analysis were selected and annotation was checked using the sorghum reference assembly. The gene sequences were further blasted against maize and rice for validation and identification of potential candidates. Additionally, we used the high resolution, open access sorghum QTL Atlas (Mace et al., 2019) to identify and characterize overlapping abiotic stress QTL.

We compared the sequence of S. bicolor reference genome v3.1.1 and the sequence based on the alternate allele discovered at the SNP position Sb10-5394955 with the recently published sorghum pangenome data (Tao et al., 2021). The gene sequence was selected from the reference genome (v3.1.1) and blasted (blastn: 2.2.31; Camacho and Madden, 2013) against the individual thirteen assembled genomes from the sorghum pangenome dataset. The fasta sequence of the results were individually trimmed from the respective genome sequences using bedtools (v2.30.0; Quinlan and Hall, 2010) and aligned using MUSCLE (v3.8.31; Edgar, 2004). The alignment was visualized using SeaView5 (Galtier et al., 1996).

The temperature conditions and consequently the level of cold stress differed strongly among the phenotyping environments (Table 1), which is also reflected in a high degree of environmental variance on the traits (Table 2). Environments AS, PL, TEX17, and TEX18 showed a significant reduction of PHI, SY, and SN compared to the environments without thermal stress (GG and SJR; Figure 1 and Table 2). In consequence, GG and SJR were regarded as control environments for further analyses, whereas all other environments were considered and pooled as cold stress environments. Highly significant differences among the entries for the scored seed set traits were observed in both environmental groups. However, as expected, the CV was higher for the group “stress environments,” especially for the trait PHI. While genotype × environment interaction between the two control environments GG and SJR was high and of a similar magnitude as genotypic variance, it was lower for the stress environments (depending on the trait, 29–40% of genotypic variance). Heritability estimates for the groups “stress environments” and “all environments” were high for both PHI (h² = 0.72 and 0.69, respectively) and SN (h² = 0.68 and 0.63, respectively), and somewhat lower for SY (h² = 0.60 and 0.51, respectively).

The race (i.e., panicle architecture) of the genotypes did not influence SY, SN, or PHI (p ≥ 0.063), regardless of whether one-way ANOVA was applied for the pooled groups of stress environments, control environments or all environments. Comparing the public material (conversion lines and other public accessions) and the breeding lines developed under Central European conditions as groups, the latter were superior (p < 0.001) for all considered traits under the pooled stress environments, while under the control environments a difference between the groups was only detected for SY (p = 0.004). The phenotype data of all entries are provided in Supplementary Table 1.

As expected, SY, SN, and PHI were highly correlated with each other (Table 3). However, the correlation of PHI to SY and SN was stronger for the group of stress environments (r = 0.880^∗∗∗ and 0.780^∗∗∗, respectively) than for the control environments (r = 0.647^∗∗∗ and 0.578^∗∗∗, respectively). In contrast, TKW was only weakly related to these traits, implying SN as the predominant factor for SY and PHI in all environmental groups. Plant height showed a significant and positive correlation to SY, PHI and – at a lower level especially in control environments- also SN. Whereas TKW correlated with plant height in the control environments, no relation between these traits was found in the stress environments.

Trait correlations were inconsistent among the different single environments (Supplementary Table 2). While moderate correlations were observed between AS, PL, TEX18 and SJR, TEX17 showed only weak correlations to the other environments. Surprisingly, GG showed no correlation to any of the other environments (p ≥ 0.103). Pearson’s correlation coefficients of r = 0.343^∗∗∗ for SY, r = 0.442^∗∗∗ for SN and r = 0.385^∗∗∗ for PHI were observed between the pooled groups of stress and control environments.

The effect of flowering time on the yield-related traits was generally weak and inconsistent over the environments (Supplementary Table 3), as expected due to the vast ecogeographical variation between the testing locations. For the control environment GG, there was no effect of flowering time at all (p ≥ 0.053). Whereas a significant and moderately negative correlation between days to flowering and TKW was observed for the temperate, short-season stress environments of AS and PL, this effect was much weaker for the tropical stress environments (insignificant for TEX17 and weak for TEX18). Also for PHI, later flowering genotypes showed a slight disadvantage at three out of four stress environments. In contrast, later flowering resulted in no penalty for SN nor SY (with the exception of location PL).

A substantial proportion (n = 110) of the diversity set had previously been extensively scored for juvenile cold tolerance and early vigor, in independent field trials and climate chamber experiments (Schaffasz et al., 2019). However, only few, weak and inconsistent correlations between juvenile and reproductive cold tolerance could be found. Seedling emergence of an early-sown field experiment in Giessen, Germany in 2014 showed a positive, but weak correlation to mean SN (r = 0.197^∗) and PHI (r = 0.200^∗) of the pooled stress environments in the present study. Similarly, seedling emergence under controlled cold stress (13/10°C) was also positively related (r = 0.204^∗) to mean PHI of stress environments. In contrast, mean juvenile biomass of four field trials (Schaffasz et al., 2019) correlated negatively with mean SN (r = −0.221^∗) and mean PHI (r = −0.198^∗) of the pooled stress environments in the present study. Furthermore, leaf greenness, a score for photosynthetic activity under controlled and prolonged cold stress, was negatively related to PHI as well (r = −0.238^∗).

To identify accessions with a superior and stable reproductive cold tolerance, the 5% (n = 12) best performing genotypes for each trait, based on their mean over all four stress environments [for which only entries with a complete dataset (n = 240) were considered], were further dissected for their stability, expressed by the CV across the mean trait values of a line in each environment. Figure 2 shows that breeding lines SB14109 and SB14122 attained the highest PHI scores. Among the public material, accessions PI53578101 (from United States), PI550586 and PI550599 (both from Russia) showed the best results. Regarding SN, the conversion line SC37 (originating from Ethiopia) reached the highest mean value, but a rather low environmental stability due to a poor performance in Texcoco2017. In contrast, SC20 (also originating from Ethiopia) showed a very high environmental stability at the expense of just a slight reduction in mean performance, hence appearing to be the preferable tolerance source. Further, public lines SC701 (from Sudan) and PI176767 (from Turkey) ranked among the top 5% with regard to reproductive cold tolerance. As expected due to the high trait correlations, SC37, SC701, and SC20 (again at a superior environmental stability) were among the top-performers for SY. Breeding line SB14022 reached the highest mean SY, but showed a rather high environmental fluctuation. Moreover, public lines SC1069 and SC566 (both from Nigeria) were among the top performers for SY under cold stress. Chinese kaoliang types (n = 15), which are known for their superior juvenile cold tolerance, showed only moderate levels of reproductive cold tolerance (Supplementary Table 1). This finding is concordant with the low and inconsistent correlation between juvenile and reproductive cold tolerance (see section “Trait Correlations”).

As presented in Table 4, ten lines ranked among the top 5% of performers for more than one trait. Of these, line SB14022 excelled in all three traits simultaneously and is a strong candidate for further focus in breeding.

The phylogenetic relatedness of the diversity panel, and the localization of genotypes with superior reproductive cold tolerance is depicted as a phylogram in Figure 3, and as a principal component analysis in Supplementary Figure 1. Genotypes with superior reproductive cold tolerance tend to cluster within related branches. However, these branches are rather evenly distributed across the phylogram, indicating that sources of reproductive cold tolerance are not confined to a particular origin.

Using pooled phenotype data from the stress environments AS, PL, TEX17 and TEX18, GWAS identified a total of 55 marker-trait associations for PHI (5 associations), SN (21), SY (9), TKW (9), DTF (1), and PH (10), respectively, (Figure 4 and Supplementary Figures 2, 3). Full details of chromosome positions of identified MTAs (marker trait associations) for each trait, including explained variance (R²) values and p-values, are provided in Supplementary Table 4. Due to the high phenotypic correlations among the SY traits, some overlapping associations were detected between the traits under different environments (Supplementary Table 5). The average phenotypic variance explained (R²) was 4.34 for the stress environments.

Considering the most relevant traits in regard to reproductive cold tolerance, namely PHI and SY, marker-trait associations were identified on chromosomes Sb01, Sb05, Sb06, Sb09, and Sb10 (Supplementary Figure 2 and Supplementary Table 4). SNP marker Sb10-5394955 on chromosome Sb10 was associated to both PHI and SY (and also associated with TKW) in the cold-stress environments (Figure 4).

The SNP marker Sb10-5394955 on chromosome Sb10, which showed multiple cold-related trait associations to PHI, SY, and TKW, showed relatively weak LD to the adjacent SNPs on either side (R² < 0.8; Figure 4), indicating recombination within the haploblock region. Based on the S. bicolor reference genome assembly (S. bicolor reference genome v3.1.1; https://phytozome.jgi.doe.gov/pz/portal.html#!info?alias=Org_Sbicolor), the haplotype block containing Sb10-5394955 overlapped with the gene Sobic.010G067300 and no other genes were present in the region. Strong synteny to the very well-annotated public reference genomes of rice³ and maize (B73 RefGen_v5; https://www.maizegdb.org/genome/assembly/Zm-B73-REFERENCE-NAM-5.0) was exploited to obtain putative functional information for this gene. Protein sequence BLAST of the genes present in syntenic rice and maize chromosome regions revealed that Sobic.010G067300 is a leucine-rich repeat receptor-like protein kinase gene. It has been found to be critically involved in cold response in rice and soybean (Gao and Xue, 2012; Yang et al., 2014).

Haploblocks surrounding each significant SNP were further investigated using the sorghum QTL Atlas (Mace et al., 2019) for overlaps to known QTL for relevant traits in other sorghum populations. QTL overlapping with the associated markers were identified for all traits under both stress and stress+control conditions (Supplementary Table 6). On each of chromosomes Sb03 and Sb10 we identified four overlapping, previously identified QTL related to cold stress overlapping the haploblock regions of traits PHI and SY (Supplementary Figure 4).

The high amount of phenotypic variation for reproductive cold tolerance observed in the present study, along with satisfying heritability estimates and identified tolerance sources with high environmental stability, suggests that a robust breeding progress for this crucial adaptation trait is feasible, even though marker-assisted selection will be challenging due to its quantitative inheritance. The cold stress environments chosen for this study represent contrasting mega-environments (temperate vs. tropical). While the lowest minimum temperatures occurred at the Mexican highland station Texcoco, the most severe stress reaction of the plant material was observed in the maritime, high-latitude environment of Asendorf (Germany). As already outlined by Schaffasz et al. (2019), this shows that other factors besides minimum temperatures also play a role. For locations Asendorf and Poel, lack of radiation, high moisture and suboptimal daily temperatures of frequently <20°C implied additional stressors, while at Texcoco radiation was not limiting and the daily temperatures were also more favorable. Nevertheless, despite these climatic differences, satisfactory heritability estimates under stress were observed for all traits (h² = 0.72 for PHI, 0.68 for SN and 0.60 for SY). Indeed, Schaffasz et al. (2019) reported even higher heritability estimates (h² > 0.80 for PHI, SN, SY) for F₁-hybrids across more or less the same set of environments as those used in the present work.

Panicle harvest index showed the highest heritability and proved to be the most suitable and reliable trait for scoring reproductive cold tolerance in sorghum, representing an efficient proxy for spikelet fertility by reducing the effects of different panicle sizes (i.e., spikelet numbers) considering panicle raw weight. In contrast to our expectations, the results showed that PHI was not influenced by race, which in sorghum is expressed in terms of contrasting panicle architecture. However, the use of PHI as a proxy for spikelet fertility might at least theoretically be distorted to a degree by compensating effects of increased TKW in panicles with poor seed set and, hence, less sink competition. The lack of negative correlation between TKW and SN for the stress environments seems to support this hypothesis, in contrast to the control environments, where a weak negative correlation between these traits was observed as expected (Supplementary Table 3). However, comparing the correlation between PHI and TKW (r = 0.321^∗∗∗) vs. PHI and SN (r = 0.780^∗∗∗) reveals that the influence of SN on PHI is much higher. In a study on sorghum heat stress, which can reduce seed set in the same way as cold stress, Singh et al. (2015) also found that the effect of stress on seed-set percentage was much higher than on TKW, excluding relevant compensation effects of an increased seed mass. Furthermore, there are reports suggesting a determination of maximum seed size of sorghum prior to anthesis based on meristem size (Yang et al., 2009).

In consequence, SN itself represents also an efficient score for reproductive cold tolerance. While it is not affected by possible fluctuations of TKW, it is obviously strongly determined by spikelet number, so that for actual spikelet fertility, it seems a less accurate score than PHI. Though, from a breeder’s point of view, a high spikelet number is desirable, reflecting satisfying yield potential, so that SN might be preferred over PHI for selection. Altogether, both PHI and SN showed a very high correlation to SY (r = 0.880^∗∗∗ and 0.914^∗∗∗, respectively) in stress environments. Since SY itself had a somewhat lower heritability, PHI and SN may represent preferable surrogate traits for selection.

Differences in flowering time can potentially complicate comparisons of reproductive cold tolerance levels across divergent genotypes, since they can lead to possible escape of earlier or later flowering genotypes from brief cold-weather conditions at the most critical growth stages. In this regard, tropical highland locations like Texcoco are particularly suitable selection environments due to their consistently low and generally stable minimum temperatures during reproductive phases, whereas temperatures tend to be more fluctuating in temperate areas. However, also for Asendorf and Poel, only weak effects of flowering time were observed on stress symptoms. We observed that later flowering entries tended to have a lower TKW in these environments, due to their shorter grain filling period. This also implies a minor disadvantage in terms of PHI, whereas SN was not influenced by flowering time.

The observed positive correlations between plant height and SY traits are in line with earlier reports (e.g., Jordan et al., 2003; George-Jaeggli et al., 2011). Except for TKW, which was probably distorted by the aforementioned compensation effects for reduced spikelet fertility under stress, the correlations with plant height were somewhat higher for the group of stress environments. The increased availability of stem reserves in taller sorghum, which is especially beneficial under stress (Blum et al., 1997), is believed to be the main reason for this phenomenon.

We identified several inbred lines with superior and environmentally stable reproductive cold tolerance. The diversity panel consisted of both publicly available material (sorghum conversion lines and other public accessions) and breeding lines developed under temperate conditions of Central Europe. As expected due to their breeding history with selection in temperate environments, the latter group performed significantly better under stress conditions. However, some public accessions also ranked among the top 5% for the stress-related traits we evaluated. These superior public accessions were of diverse origin, spanning accessions from temperate (Russia, Turkey, United States) to tropical (Ethiopia, Nigeria, Sudan) countries. While the Ethiopian highland is a known diversity center for cold tolerance (Singh, 1985), the identification of high reproductive cold tolerance in Nigerian sorghum materials may seem surprising.

However, the fact that tolerance sources from both groups (public accessions and private breeding lines) were found among genetically diverse branches over the entire phylogram (Figure 3) suggests a polyphyletic origin of reproductive cold tolerance.

The results of our study indicate that reproductive cold tolerance in sorghum is a quantitative and complex trait. This finding is consistent with the description of reproductive cold tolerance in sorghum as a heterotic trait (Schaffasz et al., 2019) and also in line with studies on the reproductive cold tolerance of rice (Jiang et al., 2010; Mori et al., 2011). Since rice diverged from maize and sorghum 50–70 Mya, it has been used as a model for sorghum using a comparative genomic approach (Roulin et al., 2009). The candidate genes discussed hereinafter have been well studied in rice, maize and other crops, indicating that a similar function in sorghum can be expected.

To detect regions associated with reproductive cold tolerance traits, we applied a combination of GWAS and comparative combined-QTL analysis, utilizing a published Sorghum QTL Atlas (Mace et al., 2019) as a reference. Combined-QTL analysis including results from different environments and populations yields a better overview of the genetic control of a trait than any single study (Rong et al., 2007). Since our study comprised environments with and without cold stress, separate analyses were performed to compare QTL detected only in cold stress environments with those detected in all (stress and control) environments. Considering cold stress environments alone, the QTL Atlas (Mace et al., 2019) revealed QTL related to abiotic stress tolerance (Supplementary Table 7 and Supplementary Figure 4) for traits such as emergence rate, seedling vigor and survival that overlapped with our identified associations. However, when focussing on only the associations consistent between the traits PHI and SY, identified in the present study as the two trait complexes most impacted by reproductive cold tolerance, we found overlaps with known QTL for dry matter growth rate, chlorophyll content and chlorophyll fluorescence. Obviously the relatively low resolution of the QTL Atlas means these findings are preliminary, however, they provide a starting point for future functional studies to dissect the genetic and physiological basis of the identified QTL.

As expected due to the high phenotypic trait correlations, GWAS revealed multiple overlapping loci among the considered SY traits. This strengthens our confidence of associations not being false positives (Shen and Carlborg, 2013; Boyles et al., 2016). While some common associations between SY traits and DTF or plant height were also observed, in line with the significant phenotypic correlations among these traits, the majority of associations for reproductive cold tolerance were not linked to a taller stature or earlier flowering. Hence, the breeding of short-statured grain sorghum with enhanced cold tolerance and maturity adapted to the temperate target environment appears to be feasible using the stress tolerance sources we identified. Presence of SNP marker highly associated to cold tolerance within the gene Sobic.010G067300, a leucine-rich repeat receptor-like protein kinase further enhances our hypothesis of this genomic region being involved in reproductive cold tolerance. A novel allele of the gene was found to be enriched in breeding lines selected for adaption to temperate climates of central Europe. Based on the results of this study, it can be speculated that the haploblock region around Sb10-5394955 may be involved in temperature adaptation of S. bicolor and could be of interest for future breeding efforts.