Phenotypic and genotypic variation in Australian native Sorghum species along aridity clines

Plant material and growth conditions

Six S. plumosum accessions, five S. stipoideum accessions and seven S. timorense accessions were germinated and grown for 10 weeks. Accessions were chosen on the basis of their collection sites occurring along aridity clines in northern Australia (Table 1). Within species, accessions were ordered numerically (P1–P6, S1–S5 and T1–T7) on the basis of increasing aridity of their collection sites (Fig. 1). The S. plumosum accessions ranged in aridity index from 0.109 to 0.920, the S. stipoideum accessions from 0.358 to 0.589, and the S. timorense accessions from 0.128 to 0.296.

Fig. 1.

(a) Mean annual aridity index (Williams et al. 2010) map of northern Australia, with collection sites of the 18 study accessions pinpointed with white circles. Accessions are labelled according to Table 1. High aridity indices (magenta) indicate less arid areas, whereas low aridity indices (cyan) indicate more arid areas. (b–d) Collection sites of the study accessions (yellow circles) are superimposed on the respective modelled distributions (blue) of (b) S. plumosum, (c) S. stipoideum and (d) S. timorense (Myrans et al. 2020).

Seeds were provided by the Australian Grains Genebank (Horsham, Vic., Australia) and were collected by the Australian Tropical Crops Genetic Resource Centre. All plants (n = 251) were grown under controlled greenhouse conditions at Monash University (−37.91, 145.14) between September and December 2020, with a mean temperature of 31.35 ± 0.09°C/21.46 ± 0.03°C, day/night, and a mean photoperiod of 13.7 h. Supplementary light was provided using sodium lamps for 13 h per day (MK-1 Just-a-shade, Ablite, NSW, Australia).

Seeds were germinated in a constant-temperature cabinet, following Cowan et al. (2020). Once the coleoptile and radicle had emerged, plants were transferred to the greenhouse and planted in seedling trays containing Debco seed-raising substrate (Evergreen Garden Care, NSW, Australia) and perlite (1:1, v/v ratio) for 1.5 weeks. Seedlings were then transplanted into 1.9 L pots containing Debco seed-raising substrate, perlite and coarse sand (2:2:1, v/v ratio). Plants were grown for a further 8 weeks and watered three times per week, namely twice with 250 mL water, once with 250 mL Thrive all-purpose soluble fertiliser (Yates, Vic., Australia) with an N:P:K ratio of 25:5:8.8 (v/v ratio). Pot positions were randomised and rotated weekly to minimise microenvironment effects. At 10 weeks post-germination, plants were harvested and separated into roots, sheaths and leaves and dried in ovens at 50°C for 72 h.

So as to confirm that stress was not a factor in any phenotypic differences, F_v/F_m, a useful indicator of plant stress, was measured on the youngest fully unfurled leaf using a MINI-PAM Photosynthesis Yield Analyser (Heinz Walz GmbH, Effeltrich, Germany) 7 weeks post-germination (n = 6). Mean F_v/F_m values were similar across all accessions, ranging between 0.78 and 0.82, indicating that all accessions were healthy and approaching the optimal F_v/F_m of 0.83 (Maxwell and Johnson 2000) (Supplementary Table S1).

Morphological parameters

Plant morphology was quantified by counting tillers, measuring total leaf area, and determining the biomass of leaves, sheaths and roots of each plant. Total leaf area was measured using a LI-3000 portable area meter and a LI-3050A belt conveyor (LI-COR Biosciences, NE, USA). Senesced leaves were included in total biomass measurements, but not in leaf-area measurements or chemical analysis. Specific leaf area (SLA; cm² g⁻¹) and root:shoot ratio were calculated using the following equations:

Relative growth rate (RGR; g g⁻¹ day⁻¹) was calculated using the equation:

where W is whole-plant biomass, and t is time. W₁ was measured at the start of the experiment (t₁), at the seed stage, with the mean W₁ for the accession used in every RGR calculation for that accession. W₂ was obtained by measuring the biomass of the entire plant at the end of the treatment period (t₂).

Chemical analysis

HCN potential was measured through colorimetric assay in 96-well microtiter plates, by using the evolved cyanide method, following Gleadow et al. (2012). Briefly, 10 mg of dried, ground tissue was mixed with 300 μL of β-glucosidase solution to catabolise all its CNglcs and release HCN. The HCN reacted with 200 μL of 1 M sodium hydroxide solution to form sodium cyanide, the concentration of which was then measured through mixing 50 μL with 125 μL of 2 M acetic acid, 125 μL of succinimide reagent and 50 μL of barbituric acid dissolved in pyridine before incubating at room temperature for 15 min. Absorbance was then measured at 590 nm (FLUOstar Galaxy, BMG Labtechnologies, Vic., Australia).

Phenolic concentration was also measured through colorimetric assay in 96-well microtiter plates, by using a Folin–Ciocalteu method adapted from Ainsworth and Gillespie (2007) and Bärlocher and Graça (2020). Phenolics were extracted from 10 mg of dried, ground tissue into 1 mL of 70% acetone solution. 20 μL of sample was then mixed with 20 μL distilled water, 20 μL dilute Folin–Ciocalteu reagent and 200 μL of 2% sodium carbonate in 0.1 M sodium hydroxide solution, before being incubated at room temperature for 2 h. Absorbance was then measured at 760 nm (FLUOstar Galaxy, BMG Labtechnologies, Vic., Australia). Gallic acid (91215, Sigma-Aldrich, MO, USA) was used as the standard.

Silicon concentration ([Si]) was measured in leaf, sheath and root tissue by using 80 mg dried, ground tissue through X-ray fluorescence spectrometry (Epsilon-3x, PANalytical, Netherlands), following Reidinger et al. (2012) and Hiltpold et al. (2016), at the Hawkesbury Institute for the Environment, Western Sydney University (NSW, Australia). Silicon was calibrated against a certified plant reference material of known silicon concentration (NCS ZC73018 Citrus leaves, China National Institute for Iron and Steel). Plant silicon concentration (plant [Si], %) was calculated using the following equation:

Statistical analysis

Within species, variation in phenotypic traits among accessions was analysed using one-way ANOVAs. Where significant variation was found, post hoc comparisons were made using Tukey’s tests. Covariance of traits with aridity was analysed using linear regression, by using aridity index as the explanatory variable and the mean trait value for each accession as the response variable. Slope and R² were calculated to give the direction and strength of relationships. Linear regression results are presented in Table S2. For all analyses, P < 0.05 signified statistical significance. Where raw data were not normally distributed, data were log-transformed to reduce skewness. ANOVA and linear regression analyses were performed using GraphPad Prism ver. 9.0.2 (GraphPad Software, San Diego, CA, USA).

Similarity between suites of phenotypic traits among accessions was analysed using principal-component analysis followed by K-means cluster analysis in R ver. 4.0.0 (R Core Team 2020), using accession means for specific leaf area, herbivory tolerance traits, and the most biologically relevant chemical defence metrics collected (plant silicon concentration, root HCN potential and leaf phenolic concentration). Mean trait values were transformed to z-scores (mean = 0, s.d. = 1) and the silhouette method was used to determine the optimal number of clusters (six) on the basis of having the highest average silhouette width (Kaufman and Rousseeuw 2009).

Genetic analysis

Genotypes were compared using SNPs obtained through DArTseq™ (Jaccoud et al. 2001; Kilian et al. 2012). Leaf samples were collected from five plants of each study accession and freeze dried. Samples were sent to Diversity Arrays Technology Pty Ltd (Bruce, ACT, Australia) for DNA extraction and genotyping. The raw data comprised 169 228 codominant, genome-wide, biallelic SNPs. Data filtering was performed using the R package ‘DARTR’ (ver. 2.9.7) (Gruber et al. 2018), following a method adapted from Gruber et al. (2019). First, we removed loci with less than 99% reproducibility. Then we removed monomorphic loci and loci with call rates below 95%. Last, we removed all but one SNP locus from each sequenced fragment. Our final dataset, which was used for analysis, contained 10 159 SNP loci with 1.64% missing data. We analysed population structure by running a principal-coordinate analysis (PCoA) in DARTR, as well as estimating pairwise F_ST (Weir and Cockerham 1984) in the R package ‘HIERFSTAT’ (Goudet 2005).

Biochemical variation among accessions

In all study accessions, excepting P1, T5, T6 and T7, mean HCN potential was higher in the roots than in the leaves or sheaths. Within accessions, mean HCN potential often varied among different organs across several orders of magnitude. For example, in S. plumosum accession P4, mean leaf, sheath and root HCN potentials were 1.20 μg g⁻¹, 2.59 μg g⁻¹ and 119.91 μg g⁻¹ respectively. Mean leaf HCN potential ranged from 0.57 μg g⁻¹ in T5 to 5.03 μg g⁻¹ in S3; mean sheath HCN potential ranged from 0.48 μg g⁻¹ in S1 to 53.23 μg g⁻¹ in P1; and mean root HCN potential ranged from 11.96 μg g⁻¹ in T5 to 119.91 μg g⁻¹ in P4 (Fig. 2). Root HCN potential significantly varied among S. plumosum accessions, but there was no significant co-variation between this trait and aridity. There was no significant variation in root HCN potential among accessions of S. stipoideum or S. timorense (Table S2).

Fig. 2.

HCN potential (HCNp) in the (a–c) leaves, (d–f) sheaths, and (g–i) roots in six S. plumosum accessions, five S. stipoideum accessions and seven S. timorense accessions. Plants were grown for 10 weeks in standardised greenhouse conditions. Data were analysed by one-way ANOVA using log-transformed data. F- and P-values are presented, with significant P-values in bold. Columns with different letters are significantly (P < 0.05) different. Values are means ± s.e.

Leaves had higher concentrations of phenolics than did the sheaths or roots in all study accessions, ranging from 13.72 mg g⁻¹ in S1 to 32.38 mg g⁻¹ in P1 (Table S3). Leaf phenolic concentration significantly varied among accessions in all three study species (Fig. 3). However, the relationship between aridity and phenolic concentration was inconsistent across organs and species, with significant positive correlations found in the roots of S. plumosum and leaves of S. stipoideum, and a significant negative correlation found in the sheaths of S. timorense (Table S2).

Fig. 3.

(a–c) Leaf phenolic concentrations and (d–f) plant silicon concentrations in 18 accessions of three wild Sorghum species grown for 10 weeks in a common garden. Accessions are labelled according to Table 1. Data were analysed by one-way ANOVA. F- and P-values are presented, with significant P-values in bold. Columns with different letters are significantly (P < 0.05) different. Values are means ± s.e.

Mean plant silicon concentration ranged from 0.35% in T6 to 0.64% in P1 (Fig. 3, Table S4). Plant silicon concentration significantly varied among S. plumosum accessions. Low-aridity accession P1 had the highest concentration, but there was no significant variation among the remaining accessions or significant linear correlation with aridity. In S. stipoideum, plant silicon concentration significantly decreased as aridity increased (Table S2). There was no significant difference in plant silicon concentration among accessions of S. timorense.

Morphological variation among accessions

Mean biomass ranged from 6.5 g in P1 to 18.8 g in P3. Mean relative growth rate ranged from 0.094 g g⁻¹ day⁻¹ in S2 to 0.117 g g⁻¹ day⁻¹ in P4. Both biomass and relative growth rate significantly varied among accessions in S. plumosum and S. stipoideum, but neither significantly co-varied with aridity. No significant variation in either trait was found in S. timorense (Fig. 4).

Fig. 4.

(a–c) Whole plant biomass and (d–f) relative growth rate (RGR) in 18 accessions of three wild Sorghum species grown for 10 weeks in a common garden. Accessions are labelled according to Table 1. Data were analysed by one-way ANOVA. F- and P-values are presented, with significant P-values in bold. Columns or points with different letters are significantly (P < 0.05) different. Values are means ± s.e.

Mean specific leaf area ranged from 83.5 cm² g⁻¹ P4 to 197.4 cm² g⁻¹ in T7. Specific leaf area significantly varied among S. plumosum accessions, but there was no significant co-variation with aridity. There was no significant variation in specific leaf area among accessions of S. stipoideum or S. timorense. Mean number of tillers was generally high, ranging from 23.6 in P3 to 95.4 in S1. Significant variation in number of tillers was found among accessions in all three study species. However, this trait did not significantly co-vary with aridity in S. plumosum or S. stipoideum. In S. timorense, significant positive correlation was found (Table S2). Mean root:shoot ratio ranged from 0.11 in T7 to 0.32 in P6. Significant variation in root:shoot ratio was found among accessions in S. plumosum and S. timorense, but not S. stipoideum. In S. plumosum, root:shoot ratio had a significant positive correlation with aridity, whereas in S. timorense a significant negative correlation was found (Table 2, Table S2).

Variation in suites of phenotypic traits among accessions

The K-means cluster analysis separated the 18 study accessions into six clusters according to phenotype. No species had all its accessions grouped into a single cluster, and three clusters contained accessions of multiple species (Fig. 5). Clusters 2 and 5 collectively contained every S. timorense accession and were both characterised by not varying greatly from the mean for any parameter, compared with other clusters. Cluster 6 contained the high-aridity S. plumosum accessions P4 and P6 and was characterised by high root HCN potentials and root:shoot ratios, and low specific leaf areas and numbers of tillers. Low-aridity accessions P1 and S2 formed Cluster 1, characterised by high silicon and low relative growth rates. Cluster 4 contained only accession P3, characterised by a high root HCN potential, and a low number of tillers and specific leaf area. Cluster 3 contained the remainder of the S. plumosum and S. stipoideum accessions. Like in Clusters 2 and 5, the mean z-scores in Cluster 3 did not vary greatly from the overall mean for any parameter (Fig. 5).

Fig. 5.

K-means cluster analysis of seven phenotypic traits from 18 accessions of three wild Sorghum species grown for 10 weeks in a common garden. (a) Cluster analysis results plotted on the first two principal-component (PC) axes and coloured according to cluster. Accessions are labelled according to Table 1. (b) Bar graph showing clusters’ mean z-scores for each trait.

Genotypic variation among accessions

The genotype PCoA separated the 18 accessions into four major clusters across the first two principal components; all S. stipoideum accessions formed one cluster; all S. plumosum accessions formed one cluster; accessions T1–4 formed a cluster; and accessions T5–7 formed the final cluster (Fig. 6a). Across the third and fourth principal components, the S. timorense accessions split into the same clusters, but S1 and S2 formed a separate cluster from the other S. stipoideum accessions, and P1 separated from the S. plumosum cluster (Supplementary Fig. S1). The pairwise F_ST comparisons showed similar patterns. Every S. plumosum accession’s highest F_ST comparison arose when compared with P1 (Fig. 6b). Similarly, all S. stipoideum accessions other than S1 gave their highest F_ST value when compared with S2 (Fig. 6c). In S. timorense, there were low F_ST values when comparing any two accessions within the T1–4 and T5–7 clusters, and high F_ST values when comparing any two accessions from different clusters (Fig. 6d).

Fig. 6.

(a) Principal-coordinates analysis of single nucleotide polymorphisms from 18 accessions of three wild Sorghum species, plotted on the first two principal-component (PC) axes. Accessions are colour coded and encompassed by ellipses. (b–d) Pairwise F_ST comparisons among study accessions within (b) S. plumosum, (c) S. stipoideum and (d) S. timorense. Numbers are in colour-coded cells, with greener and yellower shades signalling lower and higher F_ST values respectively. Accessions are labelled according to Table 1.

Morphological variation

Several morphological traits varied intraspecifically, but it was difficult to identify consistent trends between these traits and aridity. In general, plants had large numbers of tillers (ranging from 23.6 to 95.4), consistent with the finding of Lazarides et al. (1965) where unburnt S. plumosum populations had up to 60 tillers per plant in the wild. Number of tillers significantly increased with aridity in S. timorense, but not in the other two species. Root:shoot ratio decreased as aridity increased in S. timorense, but in S. plumosum it increased with aridity.

Although morphological variation has been found among accessions from different provenances in other grass species (Guo et al. 2017), this variation can become more visible when accessions are grown in a range of environmental conditions. For example, Quiroga et al. (2013) grew two accessions of Trichloris crinita in drought and well-watered conditions, finding that morphology did not significantly differ between accessions when grown under well-watered conditions, but did differ when grown under drought conditions. The next steps in quantifying local adaptations in wild sorghums should therefore include assessments of plasticity, such as surveys of wild plants, reciprocal transplant experiments, and growth of accessions in multiple conditions.

Cyanogenic glucoside storage consistent with role in nitrogen metabolism

In 14 of the 18 study accessions, mean HCN potential was higher in roots than in aboveground organs, mirroring patterns found in other wild sorghum species (Cowan et al. 2020; Myrans et al. 2021). CNglcs being stored away from the organs most vulnerable to herbivore attack, and at low concentrations, suggests that defence might not be the primary function of CNglcs in these species, at least when mature. The consistency of low aboveground HCN potentials in Australian wild sorghums, regardless of treatment or provenance (Cowan et al. 2020; Myrans et al. 2021), highlights their potential utility in development of safer forage sorghum suited to a range of environments.

It has been proposed that CNglcs could serve as a store of reduced nitrogen on top of their defensive role (Møller 2010), and storage of CNglcs in the roots would reduce nitrogen loss in the event of herbivory or fire. The nitrogen stored in these CNglcs could then be remobilised for use in compensatory growth following biomass losses (Selmar et al. 1988; Zagrobelny et al. 2008). It is also possible that root CNglcs could defend these plants against root herbivores such as nematodes (Curto et al. 2012; Moore and Johnson 2017). However, it is unclear whether their HCN potentials are high enough to hinder root herbivores, especially given that many non-ephemeral herbivores in Australia’s arid zones are highly specialised and can tolerate the chemical defences of their prey plants (Smith and Morton 1990).

While root HCN potential significantly varied among S. plumosum accessions, there was again no clear relationship between this trait and aridity. There was no significant variation in root HCN potential among accessions of S. stipoideum or S. timorense. These results contrast those of Kooyers et al. (2014), who found that the likelihood of CNglc production in T. repens grown in standardised conditions was significantly higher in accessions from more arid provenances. In T. repens, CNglc production is a binary trait controlled by one gene, whereas in wild sorghums it is a complex trait controlled by many genes and the environment. HCN potential is therefore a more plastic trait in wild sorghums, making relationships between HCN potential and environmental conditions harder to characterise.

More targeted research into CNglc functions in wild sorghums is needed before we can outline the selection pressures that are likely to be acting on HCN potential (Gleadow et al. 2021). One possibility is to measure HCN potential in plants grown under a range of temperatures or watering regimes so as to quantify the interactive effects of genotype and the environment and assess the plasticity of accessions grown within and outside their optimal environmental conditions.

Concentrations of phenolics and silicon

It was unexpected that the relationship between phenolic concentrations and aridity would vary among species. In the roots of S. plumosum and leaves of S. stipoideum, phenolic concentrations increased with aridity, whereas in the sheaths of S. timorense, phenolic concentration decreased as aridity increased, with no significant relationships detected in other organs. In contrast, Carvalho et al. (2022) found that arid-adapted Arctotheca calendula populations had consistently higher leaf phenolic concentrations, and proposed that this could mitigate heat and drought stress. Our result further highlighted the difficulty of predicting plant defence traits on the basis of environmental parameters alone. Like CNglcs, phenolics are thought to play multiple roles in plant metabolism, including defence against herbivores and stress mitigation (Lattanzio et al. 2012). These multiple roles can confound relationships with the environment, making trends difficult to identify.

Plant silicon concentration significantly varied among accessions in S. plumosum, with its low-aridity accession, P1, having the highest concentration. In S. stipoideum, plant silicon concentration significantly decreased with an increasing aridity. High silicon concentrations might be especially beneficial as a defence in tropical, rainy environments, where there are typically higher concentrations of insect herbivores (Ryalls et al. 2018). Plant silicon concentrations in this study were comparable to those found previously in both domesticated and wild Sorghum taxa (Simpson et al. 2017; Myrans et al. 2021).

Aridity, phenotype and genotype did not consistently covary

Of the 10 phenotypic traits we tested, S. plumosum showed significant intraspecific variation in eight, S. stipoideum in four, and S. timorense in four, supporting our expectation that S. plumosum would display the most phenotypic variation of the study species. This result could have stemmed from S. plumosum covering a broader distribution than the other study species (Myrans et al. 2020). Despite this, significant linear relationships between traits and aridity were not common in S. plumosum, and were no more common than in the other two study species, implying that other abiotic and biotic factors might be more important than aridity for S. plumosum’s high phenotypic diversity. Further research into the study accessions would also allow us to better estimate the aridity of provenances during their growing seasons, rather than using an annual average.

No species had all its accessions grouped into a single cluster in the K-means cluster analysis, and the three largest clusters all contained accessions of multiple species. This suggests that not all accessions within species had the same suites of phenotypic traits and could indicate some intraspecific local adaptation. Moreover, there were occasionally more similarities between accessions from different species than within a species. Evidence of phenotypic similarity not necessarily reflecting genetic relatedness has also been found in Asclepias L., implying the occurrence of convergent evolution (Agrawal and Fishbein 2006).

Sorghum timorense was the only study species to be split over two clusters across the first two principal components of the genotype PCoA, and it had the greatest differences in F_ST values among accessions. These results suggest that there is greater genotypic variation among study accessions in S. timorense than among those of the other two species, and implies the existence of a partial reproductive barrier between accessions T1–4 and T5–7. The finding was surprising considering that S. timorense is considered the most habitat-specific of the study species (Lazarides et al. 1991), and that T6 was separated from T5 and T7 in the phenotype cluster analysis despite their similar genotypes.

The cause of the restricted gene flow between accessions T1–4 and T5–7 is not yet known, with no obvious geographical barriers between the groups, and no immediately evident differences in soil type (Table 1), geology or hydrology (Boone et al. 2022). More extensive genotyping across the species’ range may help us better understand its genetic differentiation. The relative phenotypic similarity of T6 to T1–4 could have resulted from the habitat and soil similarity across S. timorense’s range, preventing phenotypic differentiation from being necessary.

The disequilibrium between genotypic and niche diversity in our study accessions challenges broader assumptions often used in conservation, such as that made by Myrans et al. (2020), that broader geographic distributions confer greater phenotypic and genotypic diversity in wild species. One potential explanation is that restricted gene flow is crucial for the creation of local adaptation (Kawecki and Ebert 2004), meaning that the largely contiguous distributions of S. plumosum and S. stipoideum might allow sufficient gene flow to prevent sizeable local adaptation, as evidenced by their relative genetic similarity among accessions. Going forward, investigation into phenotypic and genotypic differentiation of S. plumosum and S. timorense populations in Indonesian islands, compared with those in mainland Australia, could provide further information on the role of gene flow restriction in differentiation in wild sorghums.