The 245 sorghum accessions were used for forage quality characters as described previously (Li et al., 2018). The sorghum accessions were grown in four environments (2 locations × 2 years), i.e., Fengyang campus of Anhui Province (Fengyang, China, 32°52′N, 177°33′E) in 2015 and 2016, and Tengqiao town of Hainan Province (Tengqiao, China, 18°24′N, 109°45′E) in 2016 and 2017. All experiments were performed using a completely randomized block design with three replicates each. Sorghum cultivar Tx430 was used to perform RNA-seq. It was grown at Fengyang campus.

The four yield traits (PH, TN, SD, and FW) were measured when all accessions were at the heading stage. The middle stem of each plant was used to measure SD, and only aerial plant parts were used to determine FW.

Total DNA was extracted using a DNAsecure Plant Kit (Cat. No. DP320, Qiagen, Hilden, N.W, Germany). Library construction, restriction site-associated DNA (RAD) sequencing, and SNP analysis were performed as described previously (Li et al., 2018).

Linkage disequilibrium (LD) analysis was performed using PopLDdecay, with a MaxDist of 1,000 kb. All SNPs were filtered for population structure (Q), and relative kinship analysis (K) was performed using Plink v1.07 (MAF < 0.05, r² = 0.2; Purcell et al., 2007). Number of clusters in the population (k) was set from 1 to 10, with five independent runs (Pritchard et al., 2000).

GWAS was performed using TASSEL 5.2.70 (Bradbury et al., 2007), with a mixed linear model (MLM) to calculate associations and the incorporation of Q matrix/PCA and kinship data (K; Zhao et al., 2011). The MLM was applied using default settings (P3D for variance component analysis and compression set to the optimum level). For MLM (Q + K), the significance threshold for significantly associated markers was set to p ≤ 4.06 × 10^–4 or [-log10 (p) = 3.39], as described previously (Li et al., 2018).

Four weeks after planting, Tx430 leaves were sampled every 2 weeks until 10 weeks after planting, representing stages 1–4, respectively. The samples were flash-frozen in liquid nitrogen and stored at –80°C before RNA extraction. Each sample had three biological replicates. Total RNA was extracted using an RNAprep Pure Plant Kit (Tiangen, Beijing, China). Gel electrophoresis and a BioDrop (Biochrom, Cambridge, London, United Kingdom) were used to measure the quality and quantity of total RNA. Libraries were constructed and sequenced at the Beijing Genomics Institute.

Raw data were initially filtered using SOAPnuk v1.5.2 (Chen et al., 2018), and then histat2 (Kim et al., 2019) was used to map clean reads to a sorghum reference genome of BTx623 version 3.1.1 (McCormick et al., 2018).² Differential expressed genes (DEGs) were identified using R package DESeq2 with a padj < 0.05-Benjamini-Hochberg multiple test correction (FDR) and the absolute value of a log2 (FC) > 1 (Love et al., 2014). Both gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were performed using the “clusterProfiler” package in R (Wu et al., 2021). The enrichment results of GO and KEGG pathways were obtained using P < 0.05 as the significance threshold. The top 5 KEGG pathways and top 5 terms of each GO domain were identified.

Weighted gene correlation network analysis (WGCNA) was performed using the R package “WGCNA” (Langfelder and Horvath, 2008). Firstly, the genes were ranked by median absolute deviation from large to small, and the top 50% genes were selected for WGCNA using the “goodSamplesGenes” function in package “WGCNA.” Subsequently, the power parameter ranging from 1 to 20 was screened out using the “pickSoftThreshold” function in package “WGCNA.” A suitable soft threshold of 8 was selected, as it met the degree of independence of 0.85 with the minimum power value. Finally, modules were obtained following dynamic branch cutting with a merging threshold of 0.25. The modules were visualized by the “plotDendroAndColors” function in package “WGCNA.” The correlation map between modules and traits was visualized using the R package “ggcor.”

Stable QTLs were those detected across two environments. Genes within 50kb of QTL-associated SNPs were considered for further analysis based on the LD results. The candidate genes in the QTLs were obtained according to the reference genome (Sorghum bicolor v3.1.1, McCormick et al., 2018) and annotation GFF3 file (Sorghum bicolor v3.1.1)³ using BEDTools (Quinlan and Hall, 2010).

The phenotypic mean data of PH, TN, SD, and FW were calculated with “Descriptive Statistics” pack in “Data Analysis” tool using Excel 2010. Correlation analysis and histogram construction were performed for the traits using the “PerformanceAnalytics” package in R.

Extensive variation in PH, TN, SD, and FW was observed in all four environments in the 245 accessions (Table 1). The extent of variation for the traits ranged from 1.5- to 5.6-fold. PH ranged from 90.0 to 476.7 cm, with 3.3–5.3-fold variation in the different environments, whereas TN ranged from 0 to 6.3, with 3.5–5.6-fold variation, SD ranged from 5.10 to 29.55 mm, with 1.5–1.8-fold variation, and FW ranged from 0.073 to 5.830 kg, with 2.4–4.9- fold. Mean PH, SD, and FW were significantly lower at Tengqiao (Tq) than at Fengyang (Fy), but no significant difference was observed for TN. It suggests that PH, SD, and FW, but not TN, were affected by photoperiod which was shorter at Tengqiao than at Fengyang.

Furthermore, PH, FW, and SD were normally distributed, whereas the distribution of TN was relatively skewed in all four environments. According to Pearson’s correlation coefficients, FW was significantly and positively correlated with PH, SD and TN in all four environments, whereas TN was significantly and negatively correlated with SD (Figure 1). It suggests that the traits are genetically linked or that the traits are affected by genes with pleiotropic effects.

LD in the 245 accessions was calculated using parameter r² with 3,026 SNPs. The LD in the 245 accessions decayed after 25 kb on average (Figure 2A), which suggests that the QTLs detected in multiple environments were less than 25 kb from the causal marker(s). The best K-value in the population structure was 7 (Figure 2B) and was used for GWAS analysis.

GWAS was performed using a MLM and 85,585 SNP markers. A total of 338 SNPs, or quantitative trait nucleotides (QTNs), were associated with the four yield traits with phenotypic variation explained (PVE) ranged from 4.1 to 57.07% (Supplementary Table 1). Seventy-four SNPs were associated with PH, and 97, 83, and 84 were associated with TN, SD and FW. The association p-values for the QTNs ranged from 9.38E-04 to 7.90E-12. The QTNs were distributed relatively evenly across the 10 chromosomes with the highest of 38 QTNs on chromosome 3 and the lowest of 29 on chromosome 6. The numbers of QTNs detected in each of the four environments were not significantly different: 90 QTNs in 2015Fy, 85 in 2016Tq, 87 in 2016Fy, and 76 QTNs in 2017Tq.

In addition, 21 stable QTNs were detected in two environments at least. Four PH QTNs were detected on chromosomes 1, 3, 4, and 8 (one each) in two environments. Ten TN QTNs were detected on chromosomes 3, 4, 5, 7, and 10 in all four environments. Six SD QTNs were detected on chromosomes 3, 7, 8, and 10. One FW QTN was detected on chromosome 8. In these stable QTNs, the PVE of eight QTNs was greater than 10% in both environments, and the PVE of S4_1261758 and S3_69018585 were greater than 20% (Table 2).

Number of DEGs decreased with development stage from 5,456 (2,417 up- and 3,039 down-regulated) between stages 1 and 2 (stage1_stage2) to 1,246 (684 up- and 562 down-regulated) between stages 2 and stage 3 (stage2_stage3), and 1,379 (543 up- and 836 down-regulated) between stages 3 and stage 4 (stage3_stage4, Figure 3). In addition, the three sets of DEGs included 4,443, 573, and 711 unique DEGs, respectively. However, 106 DEGs were shared among all three-stage comparisons, which suggested that the DEGs play important roles in vegetative development.

Furthermore, KEGG enrichment analysis indicated that the stage1_stage2, stage2_stage3, and stage3_stage4 DEGs were associated with pathways related to ribosome and amino acid biosynthesis, stress and glutathione metabolism, and circadian rhythm and photosynthesis, respectively (Figure 4). From stage3 to stage4, sorghum transitioned from vegetative to reproductive development. Four Flowering Locus T-like (FTL) genes were significantly up-regulated at stage 4, and the expression of two phytochrome biosynthesis-related genes were also affected in the stage3_stage4 DEGs (Table 3). The FTL genes may play a role in the transition between vegetative to reproductive stages. The 106 DEGs shared among the four developmental stages encoded proteins involved in plant hormone signal transduction, including three jasmonate-zim-domain proteins that were down-regulated in stage1_stage2 and stage3_stage4 but up-regulated in stage2_stage3 (Table 3). In Arabidopsis, JAZ10/JAZ11 regulates root growth (Liu et al., 2021).

Weighted gene correlation network analysis (WGCNA) indicated that 45 separate modules were correlated with PH and SD (Figure 5A and Supplementary Table 4). Eight of the modules were significantly correlated with PH, and three with SD (p < 0.05; Figure 5B). The blue module had the highest positive correlation with PH and SD (r = 0.914 and r = 0.915, respectively), and the turquoise module had the highest negative correlation with PH and SD (r = –0.717 and r = –0.958, respectively). There were 4,733 and 8,635 genes in the blue and turquoise modules, respectively, which indicates that the development of PH and SD is complex and involves a large number of genes.

The DEGs in the blue module were up-regulated in all three-stage comparisons, whereas the most of DEGs in the turquoise module were down-regulated in the stage1_stage2 and stage3_stage4 (Supplementary Table 5). KEGG enrichment analysis indicated that the genes in the blue module were mainly enriched in pathways related to spliceosome and protein processing in the endoplasmic reticulum (Supplementary Figure 1A), whereas the genes in the turquoise module were mainly enriched in ribosome and purine metabolism pathways (Supplementary Figure 1B).

The sorghum annotation file was used to annotate genes associated with the 21 stable QTNs. Of the 86 candidate genes associated with the 21 stable QTNs (Supplementary Table 2), PH, TN, and SD were associated with 14, 40, and 32, respectively.

Further analysis reduced the number of candidate genes to eight that were associated with seven QTNs by RNA-seq (Table 4). A gene for indole-3-glycerol phosphate synthase (Zhao, 2010), which is involved in the biosynthesis of indole-3-acetic acid (IAA), was associated with PH and was down-regulated in the stage1_stage2 comparison. Of the four candidate genes associated with TN, two were up-regulated, and two were down-regulated. Of the three candidate genes associated with SD, two (Sobic.003G047700 and Sobic.003G047800) encoded cytokinin-O-glucosyltransferase 3 and both genes were up-regulated. Another (Sobic.003G375100) encoded a mitochondrial DNA repair RAD52-like protein 1, which plays a very important role in plant development, especially vegetative development (Table 4).