A core set of 272 accessions from the NPGS sweet sorghum germplasm collection that included 99 improved germplasm accessions, 132 landraces from 22 countries, and 41 accessions from unknown origin were selected based on sugar concentration [≥9.0°Brix; Germplasm Resources Information Network (accessed August 23,2023)], sorghum race and place of origin ( Supplementary Table S1 ). This core set [named after herein as Sweet Sorghum Diversity Panel (SWDP)] included the five major sorghum races (Bicolor, Caudatum, Kafir, Durra and Guinea), represented ~10% of the NPGS sweet sorghum germplasm collection, and only shared 15 and 5 accessions with the BAP (Brenton et al., 2016) and SAP (Casa et al., 2008), respectively.

The SWDP and five reference lines that included both anthracnose resistant (SC112-14 and SC784-5) and susceptible accessions (RTx430, BTx623, and Early Honey) were planted in April 2017 and December 2018 at the research farm of the USDA-ARS Tropical Agriculture Research Station in Isabela, Puerto Rico, U.S.A., for evaluation and seed production. In summer 2017 and 2018, the SWDP and reference lines were planted at the Black Shank Farm of the University of Georgia in Tifton, Georgia, U.S.A. (31°29’58.9”N, 83°32’53.8”W), at the Research Farm of Texas A&M University at College Station, Texas, U.S.A. (30°31’55.5”N, 96°25’23.7”W), and at the University of Florida North Florida Research and Education Center-Suwannee Valley near Live Oak, Florida, U.S.A. (30°18’47.8”N, 82°54’07.8”W). The accessions at the four locations were planted in a randomized complete block design with two blocks in single-row plots of 3.1 m in length with 20-30 plants and 0.9 m between rows.

Even though the high relative humidity at these four locations makes anthracnose an endemic sorghum disease, several plants per plot were manually inoculated to have a uniform disease distribution in the field. The inoculations were similar to those described by Prom et al. (2009), where fungal cultures were prepared with different isolates of C. sublineolum that represent the pathotypes present at each location. These isolates were used to colonize autoclaved sorghum seeds that were later placed into the leaf whorl of 30-45 day- old plants. The anthracnose resistance response of each plot was determined approximately 30-45 days after flowering (hard-dough state to physiological maturity). The plots were visually inspected and the most infected plant was used to determine the diseases rating using the 1-to-5 scale as follows: 1 = no symptoms or chlorotic flecks on leaves; 2 = hypersensitive reaction on inoculated leaves, but no acervuli in the center; 3 = infected bottom leaves with acervuli formation in at least one plant within the plot; 4 = necrotic lesions with acervuli on bottom leaves showing spreading to middle leaves in at least one plant within the plot; and 5 = most leaves dead because of infections and flag leaf infected in at least one plant within the plot (Prom et al., 2009). An analysis of variance (ANOVA) for anthracnose resistance response was performed using the proc mixed covtest method type3 procedure of SAS (version 9.4, SAS Institute, Cary, NC). In the mixed linear model for such ANOVA the locations and years were treated as fixed effects and accessions as a random effect with the following model: A = μ + Y + L + L×Y + G + G×L + ϵ; where A is anthracnose resistance response, μ is the common effect, Y represents the effect of the year, L is the location effect, Y×L is location × year interaction, G is the effect of the accessions, L×G is the location × accession interaction and ϵ is the error term (residual).

The anthracnose resistance response of each accession was also categorized as resistant (score ≤ 2.0) or susceptible (score > 2.0), and the frequency of resistant accessions was determined for each population and the SWDP as a whole in each location. To examine if any of the populations within the SWDP had been selected for anthracnose resistance, we performed a χ² test using the observed number of resistant genotypes within each population and the expected number based on the frequencies in the entire SWDP. Segregation distortion against the expected ratio [χ² P(value) < 0.01] was considered a signature of selection. To determine the correlation between photoperiod-sensitivity and anthracnose resistance, the flowering time of each plot at Tifton, GA was recorded as the number of days post planting until 50% of the plants within a plot reached anthesis. The Pearson correlation coefficient was calculated in SAS (Cary, NC) using an assigned value of 300 days to photoperiod-sensitive accessions (i.e. plants that did not flower).

Genomic DNA was isolated from seedlings of the SWDP using the method of Guillemaut and Marechal-Drouard (1992) and purified using the ZR96 DNA Clean & Concentrator-5 kit (Zymo Research, Irvine, CA, USA). A genotyping-by-sequencing (GBS) library was prepared using the restriction enzyme ApeK1 for digestion (Elshire et al., 2011) and sequenced in three lanes of an Illumina HiSeq 2500 platform at the Biotechnology Center of the University of Wisconsin, Madison, WI, USA.

The Tassel 5.0 v2 GBS pipeline (Glaubitz et al., 2014) was used to call SNPs using the most recent version of the BTx623 sorghum genome [version v.3; www.sorghumbase.com; (Paterson et al., 2009)] ( Supplementary Figure S1 ). Raw genotypes with minor allele frequencies (MAF) > 0.01 and missing data < 0.25 were filtered to retain 171,954 SNPs that were imputed using Beagle 4.1 (Browning and Browning, 2016) with a probability call of >0.80. The number of SNPs within genes, in proximity of genes (i.e. within 5 kb upstream or downstream of a gene), and in intergenic regions, as well as the number of synonymous and non-synonymous SNPs were determined using the Variant Effect Predictor (McLaren et al., 2016), as implemented in EnsemblPlants (http://plants.ensembl.org/).

Population structure of the SWDP was assessed using the model-based clustering method implemented in STRUCTURE 2.1 (Pritchard et al., 2000). The genotypic data were thinned to 2,345 unlinked SNPs (r² <0.10) using PLINK (Purcell et al., 2007) and used in the analysis. For each k value set from 1 to 15, three independent analyses were run using an admixture model with correlated frequencies, 25,000 burn-in periods, and 125,000 Monte Carlo Markov Chain (MCMC). The actual number of populations was determined using the ad hoc statistic Δk based on the rate of change in the log probability of data between successive k values (Evano et al., 2005), as implemented by the Structure Harvester software (Earl and Vonholdt, 2012). The ancestry membership coefficient of each accession was obtained by matching the three independent runs of the selected k values in CLUMPP (Jakobsson and Rosenberg, 2007). Accessions with an ancestry membership coefficient >0.60 were assigned to their corresponding population. Principal component analysis (PCA) was conducted using Tassel 5.0.

The identical-by-state (IBS) genetic distances among SWDP accessions were calculated in Tassel 5.0 using the 2,345 unlinked SNPs. The resulting matrix was subjected to a clustering analysis using neighbor-joining and visualized using Interactive Tree of Life (Letunic and Bork, 2011). The pairwise fixation indexes (F_ST) among populations were estimated in the R package HIERFSTAT (Goudet, 2005).

The genetic relationship between SWDP and the BAP was investigated through clustering and PCA analysis. Genotyping information of the BAP was obtained from the Dryad Data Repository (Hu et al., 2019) and merged with SWDP resulting in a total of 94,999 common SNPs. These genotyping data were thinned to 7,792 unlinked SNPs (r² <0.10) using PLINK (Purcell et al., 2007) and used in the analysis. The IBS pairwise genetic distance and PCA analysis among the 609 accessions from SWDP and BAP were calculated in Tassel 5.0 (Glaubitz et al., 2014). The resulting matrix was subjected to a clustering analysis using neighbor-joining and visualized using Interactive Tree of Life (Letunic and Bork, 2011).

Genome-wide association analyses were performed using 132,399 SNPs with MAF ≥ 0.05 and the Bayesian-information and Linkage-disequilibrium Iteratively Nested Keyway (BLINK) model implemented in the R-package GAPIT (Lipka et al., 2012). Log quantile-quantile (Q-Q) p-value plots were examined to determine that the inclusion of the first three ancestry coefficients adequately controlled for the population structure and family relatedness. Association analyses were performed for photoperiod-sensitivity (i.e. non-flowering under short days) and anthracnose resistance response for each location. The empirical significance threshold for BLINK based on the false discovery rate (Benjamini and Hochberg, 1995) as implemented in GAPIT using a cutoff p-value ≤0.01. Manhattan and Q-Q plots were visualized using the R package qqman (Turner, 2014).

The genomic regions in linkage disequilibrium with associated SNPs were delimited using the block function of PLINK (Purcell et al., 2007) and plots were visualized using the R package LD heatmap (Shin et al., 2006). The candidate genes within these genomic regions were identified based on the annotation of the sorghum reference genome BTx623 v.3.1.1 (https://phytozome-next.jgi.doe.gov/). Functional annotations included protein domains encoded by putative genes and associated PlantFAMs identified via hmmsearch (Finn et al., 2011) that evaluate remote protein homology as implemented in Phytozome 13.

Analysis of the anthracnose resistance response of the SWDP across the four locations revealed that the resistance response of multiple accessions is determined by C. sublineolum pathotypes ( Table 1 ; Supplementary Table S1 ). The ANOVA indicated year, location, accessions, and the interactions year × location and accessions × location have statistically highly significant effects (P ≤0.001; Supplementary Table S2 ) on the anthracnose resistance response. We observed that the frequency of resistant accessions in Puerto Rico (130 accessions) was higher than in Texas (90 accessions), Florida (80 accessions) and Georgia (67 accessions). Twenty-six accessions were consistently resistant (≤2.0; no acervuli observed) across locations, and this group of accessions includes both landraces and advanced breeding germplasm.

The GBS analysis of the SWDP identified 171,954 SNPs with an average of one SNP per 4,071 kb ( Table 2 ). These SNPs tagged 24,892 annotated genes (within 5 kb upstream or downstream) and included 22,975 non-synonymous mutations, suggesting the presence of novel alleles. After removing SNPs with a frequency lower than 0.05 (i.e. rare alleles), 132,399 SNPs were retained with an average LD block of 11,020 kb across the genome. Moreover, we observed that 16,178 SNPs were contained within 490 LD blocks ranging from 50-200 kb. The genomic comparison of this panel against the BAP identified 94,999 common SNPs. Therefore, the SWDP could be used to expand the genetic diversity of the BAP and to increase the statistical power to detect SNPs associated with under-represented traits of interest.

Population structure analysis revealed that the SWDP could be stratified into four populations and one admixed group ( Figure 1 ). In total, 188 accessions were assigned to one of these four populations, while 84 accessions were an admixture. This population structure was consistent in the neighbor-joining tree analysis where the four populations were dispersed into different clades separated by admixed clades. We observed that the population structure could be associated with the sorghum domestication process. Racial classifications based on panicle architecture showed that most of the accessions in population 1 belong to the Guinea race. In population 2, Caudataum, Guinea and their intermediates (Caudatum/Guinea) constituted 77% of the accessions. Most of the accessions in population 3 (81%) belonged to Caudatum, Durra and their intermediate (Caudatum/Durra). The Ethiopian accessions classified into the Durra race were grouped within population 4. The F_ST values showed that population 2 is highly genetically differentiated from populations 3 and 4 (F_ST > 0.27; Figure 2 ). Similar analyses of populations 1 and 4, as well as 3 and 4 also identified high levels of genetic differentiation (F_ST = 0.26 and 0.25, respectively). Most of the improved sweet sorghum germplasm was distributed within population 2 (47%) and the admixed group (27%). Remarkably, population 1 consisted only of landraces.

Phylogenetic analysis indicated that the SWDP expands the genetic diversity present in the BAP ( Figure 2 ). Multiple accessions from the SWDP clustered between the sweet sorghum accessions present in the BAP, suggesting that both panels could be combined to increase the statistical power to study sweet sorghum-related traits. Indeed, principal component analysis showed the genetic diversity within populations 2 and 3 is greater than the genetic diversity within sweet sorghum present in the BAP. The level of genetic differentiation between the high-biomass sorghums in the BAP and the SWDP was very low (F_ST = 0.02). Nevertheless, most of the SWDP germplasm forms a distinct cluster from the BAP in the phylogenetic analysis. Therefore, the SWDP contains additional genetic variation that could also be used to increase the statistical power of detection of under-represented alleles in the BAP relevant for the improvement of high-biomass sorghums.

The anthracnose resistance response was associated with SWDP population structure. When comparing the ratio of resistant to susceptible accessions at each field location where the SWDP was grown, we found evidence of selection against pathotypes from Georgia and Florida ( Table 1 ). The overall frequency of resistant accessions in Georgia (0.30) was lower than that observed in population 1 (0.59) and 3 (0.54). Likewise, in Florida the overall frequency of resistant accessions (0.36) was lower than that observed in population 1 (0.63) and 3 (0.66). Remarkably, an excess of susceptible accessions in population 2 was observed in Florida (0.88) and Georgia (0.90). The ratio observed in Puerto Rico, Texas and across the four locations showed no evidence of selection for anthracnose resistance. The twenty-six accessions that were consistently resistant across locations represent the four populations present in the SWDP. Likewise, the photoperiod-sensitivity of accessions was associated with the SWDP population structure and anthracnose resistance response ( Supplementary Table 1 ). We observed that 98 accessions did not reach flowering stage in Georgia, 46 of which (47%) were resistant to anthracnose. In contrast, only 19 accessions (16%) of the 123 that flowered were resistant to anthracnose. We observed selection for photoperiod-sensitive germplasm in populations 1, 3 and 4, and selection against photoperiod-sensitivity in population 2. The admixture group showed no evidence for selection for flowering time.

The photoperiod-sensitivity of the SWDP was first analyzed to validate the accuracy of the SNP data and their potential use in GWAS. The analysis identified four genomic regions in chromosomes 4, 6 and 7 associated with flowering time variation ( Supplementary Table S3 ). The positive effects of the loci S4_53762801 (p-value = 2.67 ×10^-8) and S6_39401690 (p-value = 3.44 ×10^-10)] indicate that both are associated with photoperiod response, while the negative effects of the loci S4_66999721 (p-value = 3.97 ×10^-9) and S7_16432373 (p-value = 6.39 ×10^-9)] are likely to be related to the flowering time variation among photoperiod-insensitive accessions. Indeed, the two loci associated with photoperiod-sensitive responses colocalized with the flowering genes Ma1 (Chr.6; Sobic.006G057866) and Hd1 (Chr.4; Sobic.004G216700) (Murphy et al. (2011); Supplementary Figure S2 ), which explained 21.6% of the phenotypic variance. The analysis of Q-Q plots indicated that the BLINK model including the first three ancestry coefficients from STRUCTURE (K=4) provided an adequate control of the population structure and family relatedness.

The genome-wide association scan for anthracnose resistance detected 16 loci across locations ( Figure 3 ; Table 3 ). The largest number of associations were identified in Puerto Rico, with six resistance loci on chromosomes 2 (S2_2625340), 3 (S3_4470322), 6 (S6_57256182), 8 (S8_575263 and S8_3139785) and 9 (S9_8900776), which together explained 30.5% of the phenotypic variation. We detected four loci effective against pathotypes from Texas on chromosomes 2 (S2_262340), 3 (S3_51933644) and 8 (S8_575263, S8_616316), which together explained 23.7% of the phenotypic variation. The loci on chromosomes 2 (S2_2625340) and 8 (S8_575263) were identified in Puerto Rico and Texas, and thus provide protection against pathotypes from both locations. In Georgia, the genome-wide association scan identified four loci on chromosomes 4 (S4_16285663), 5 (S5_66523668), 6 (S6_45707450) and 8 (S8_585667), which together explained 39.0% of the phenotypic variation. The locus on chromosome 8 is 10.40 kb upstream of the locus S8_57263 identified in Texas and Puerto Rico. The lowest phenotypic variation explained (11.6%) was observed against the pathotypes of Florida. The GWAS detected four loci on chromosomes 4 (S4_55304023), 5 (S5_1971654 and S5_65730857) and 7 (S7_7905713). The loci S5_65730857 and S5_66523668 detected in Florida and Georgia, respectively, are separated by only 0.8 Mb.

The size of the LD blocks of these 16 loci varied across the genome. For instance, the five loci on chromosomes 4, 5 and 6 contained LD blocks larger than 25 kb. However, the LD block for the loci S8_575263 and S8_585667 detected in Georgia, Texas and Puerto Rico were smaller than 12.5 kb. We did not detect a correlation between the size of the LD blocks and the effect of these loci.

The resistance response against pathotypes from Texas and Puerto Rico mapped to the distal end of the short arm of chromosome 8. This region contained multiple R genes ( Figure 4 ). For instance, the LD block containing the SNPs S8_575263 and S8_585667 is 1.6 kb downstream of gene Sobic.008G006800, which encodes leucine-rich repeats (LRR) and protein kinase receptor domains. Similarly, the LD block containing the SNP S8_616316 harbors the genes Sobic.008G007900 and Sobic.008G008000, both encoding LRR domains. In fact, six of the ten genes predicted in this 41 kb region have features related to disease resistance.

We identified two LD blocks on chromosome 5 associated with resistance against pathotypes from Florida and Georgia that harbor disease resistance genes. The SNP S5_65730857 and its 800 bp LD block on chromosome 5 are located in gene Sobic.005G175500, which is annotated to encode a protein tyrosine kinase domain implicated in antifungal activity. Likewise, the 27.12 kb block that contains SNP S5_66523668 harbors only one annotated gene (Sobic.005G182500.2) that encodes protein kinase and Jacalin-like lectin domains, both often associated with disease resistance responses.