Genome-wide association study reveals the genetic complexity of fructan accumulation patterns in barley grain

Abstract

We profiled the grain oligosaccharide content of 154 two-row spring barley genotypes and quantified 27 compounds, mainly inulin- and neoseries-type fructans, showing differential abundance. Clustering revealed two profile groups where the ‘high’ set contained greater amounts of sugar monomers, sucrose, and overall fructans, but lower fructosylraffinose. A genome-wide association study (GWAS) identified a significant association for the variability of two fructan types: neoseries-DP7 and inulin-DP9, which showed increased strength when applying a novel compound ratio-GWAS approach. Gene models within this region included three known fructan biosynthesis genes (fructan:fructan 1-fructosyltransferase, sucrose:sucrose 1-fructosyltransferase, and sucrose:fructan 6-fructosyltransferase). Two other genes in this region, 6(G)-fructosyltransferase and vacuolar invertase1, have not previously been linked to fructan biosynthesis and showed expression patterns distinct from those of the other three genes, including exclusive expression of 6(G)-fructosyltransferase in outer grain tissues at the storage phase. From exome capture data, several single nucleotide polymorphisms related to inulin- and neoseries-type fructan variability were identified in fructan:fructan 1-fructosyltransferase and 6(G)-fructosyltransferase genes. Co-expression analyses uncovered potential regulators of fructan biosynthesis including transcription factors. Our results provide the first scientific evidence for the distinct biosynthesis of neoseries-type fructans during barley grain maturation and reveal novel gene candidates likely to be involved in the differential biosynthesis of various types of fructan in barley.

Barley, expression analysis, fructans, 6G-FFT, grain, GWAS, neoseries, oligosaccharides, ratio-GWAS

Introduction

Starch, fructans, and (1,3;1,4)-β-glucans represent the major plant reserve carbohydrates (Vijn and Smeekens, 1999; Burton and Fincher, 2009). Among them, fructan biosynthesis has evolved polyphyletically in ~15% of higher plants, including species of the orders Asterales, Buxales, Asparagales, and Poales (Hendry and Wallace, 1993; Cairns et al., 2000; Van den Ende, 2013). In cereals, fructans accumulate in all plant organs (Pollock and Cairns, 1991).

Fructans consist of repeating fructose residues linked to a sucrose unit. The classification relates to the position of the sucrose, the linkage type between the fructose residues [i.e. β(2,1), inulin; β(2,6), levan; or containing both β(2,1)- and β(2,6)-d-fructosyl units referred to as graminan-type], and the chain lengths (Cochrane, 2000; Matros et al., 2019). Fructans can form oligomers with a degree of polymerization (DP) of 3–9 or polymers with a DP ≥10. Here, the term fructans is used to indicate either fructooligosaccharides (FOS) or fructan polymers. Fructans are typically discussed in the literature without differentiation of the DP but, since they have become more important in a dietary context (Dwivedi et al., 2014; Verspreet et al., 2015a; Liu et al., 2017), more attention has recently been paid to the role of fructans according to their DP level.

Fructan biosynthesis involves four fructosyltransferase enzymes (reviewed in del Pozo et al., 2019). First, sucrose:sucrose 1-fructosyltransferase (1-SST) mediates the transfer of fructose from one donor sucrose to the C1 position of one receptor sucrose, forming 1-kestose and glucose. Next, sucrose:fructan 6-fructosyltransferase (6-SFT) mediates the transfer of fructose from a donor sucrose to the C6 position of 1-kestose or sucrose, forming 1&6-kestotetraose or 6-kestose and glucose, thus facilitating the formation of fructans with β(2,6) linkages (levan- and graminan-type). Then, fructan:fructan 1-fructosyltransferase (1-FFT) catalyses the transfer of d-fructosyl units from one fructan to another, forming β(2,1) linkages and leading to chain elongation (inulin-type). Finally, fructan:fructan 6G-fructosyltransferase (6G-FFT) transfers a d-fructosyl unit from one fructan to the glucosyl unit of another fructan or sucrose molecule with a β(2,6) linkage and facilitates the subsequent elongation of the chain with β(2,1)- or β(2,6)-linked d-fructosyl units, allowing the respective formation of so-called neoseries-inulin-type or neoseries-levan-type fructans. Fructan hydrolysis is usually catalysed by different isoforms of fructan exohydrolyases, which release terminal fructose units from the various linkage types, namely fructan 1-exohydrolase (1-FEH), fructan 6-exohydrolase (6-FEH), and fructan 6G&1-exohydrolase (6G&1FEH).

All types of fructans are known to occur in the Poaceae (Carpita et al., 1991; Pollock and Cairns, 1991; Bonnett et al., 1997). However, Triticum, Secale, and Hordeum are believed to mainly contain branched-type fructans (graminan-type) whereas, the Poeae tribe mostly comprises levan-type fructans (Bonnett et al., 1997; Huynh et al., 2008a). Recently, the presence of graminan- and neoseries-type fructans was reported in wheat (Verspreet et al., 2015b). Neoseries-type fructans, in contrast to other fructan types, are characterized by an internal glucose unit (Matros et al., 2019; Witzel and Matros, 2020). Additional structural variations are likely to occur between different plant organs.

New developments in fructan analysis based on MS detection revealed the fine structure of cereal grain fructans with DP3–5 (Verspreet et al., 2017). Variations in fructan composition pattern and abundance were observed in oat, barley, rye, spelt, and wheat flour, suggesting a putative link between accumulation of certain fructan types and cereal phylogeny (Verspreet et al., 2017).

Reports of the beneficial health effects of fructans (Verspreet et al., 2015a; Liu et al., 2017; Andrade et al., 2019) have prompted screens for variation in their natural abundance and composition, and biotechnological approaches to increase FOS content in classical non-fructan cereals, such as maize (Dwivedi et al., 2014). However, most studies on grain fructan content still focus on wheat (Huynh et al., 2008a, b; Veenstra et al., 2017, 2019). Investigation of two doubled haploid (DH) populations (Berkut×Krichauff and Sokoll×Krichauff) revealed several quantitative trait loci (QTLs) for high fructan content in wheat grain (Huynh et al., 2008b). Winter wheat grain fructan content was found to be significantly influenced by either the genotype or the environment, as well as by genotype×environment interactions (Veenstra et al., 2019). Fructan content in developing barley grain was compared between seven genotypes, demonstrating peak accumulation between 6 and 17 days after pollination (DAP) (De Arcangelis et al., 2019) as previously reported (Peukert et al., 2014). Notably, a comparative mapping approach involving wheat and barley revealed clusters of genes encoding fructan biosynthesis enzymes (Huynh et al., 2012) on 7AS in wheat and 7HL in barley. These clusters included 1-SST, 1-FFT, 6-SFT, and several vacuolar invertase genes. Similar gene structures and physical positions of these clusters of functionally related genes in both genomes indicate that they may have evolved in parallel and that the genes within a cluster may be linked functionally in controlling fructan accumulation.

Due to its increasing potential as a health-promoting functional cereal, there is considerable interest in identifying factors that influence barley grain quality (Meints et al., 2016; Langridge and Waugh, 2019). Here we report on an analysis of natural variation in fructan content and composition across a diversity panel of two-row spring barley. We identified significant associations between fructan composition/content and fructan biosynthesis genes. We obtained support for the involvement of some of these in underpinning the observed variation from transcriptomic analysis. Additionally, potential regulators of fructan biosynthesis were assigned by co-expression analyses.

Materials and methods

Plant material

We used 154 two-row spring barley genotypes sourced from The James Hutton Institute, complemented by three Australian elite barley varieties and the wheat line Piccolo as checks (Supplementary Table S1). The germplasm was selected for minimum population structure while maintaining as much genomic diversity as possible based on principal components analysis of a much larger set of genotypes (>800). Three plants per genotype (biological replicates) were grown in a randomized main-unit design in a glasshouse compartment in a mix of clay–loam and cocopeat (50:50 v/v) and day/night temperatures of 22 °C/15 °C between July and December 2014 in The Plant Accelerator, Adelaide, Australia. Mature grains were harvested and stored until oligosaccharide analysis. For each sample, five grains were ground together to a fine powder with a PowerLyzer™24 Homogenizer (QIAGEN) and used immediately for oligosaccharide analysis.

Oligosaccharide extraction and profiling

A ‘mixed sample’ was assembled composed of equal amounts from each individual sample (154 genotypes×3 biological replicates) to capture systematic shifts during extraction and measurement. Soluble sugars were extracted following a method adapted from Verspreet et al. (2012) by incubation in 80% ethanol at 85 °C for 30 min followed by Milli-Q water at 85 °C for 30 min on a mixer (700 rpm) in a final dilution of 1:40 (w/v, mg µl^–1), and supernatants combined. Extracts were diluted with water to 1:1000 (w/v, mg µl^–1) and 25 µl per sample were analysed by high pH anion exchange chromatography with pulsed amperometric detection (HPAEC-PAD) on a Dionex ICS-5000 system using a DionexCarboPAC™PA-20 column (3×150 mm) with a guard column (3×50 mm) kept at 30 °C and operated at a flow rate of 0.5 ml min^–1. The eluents used were (A) 0.1 M sodium hydroxide and (B) 0.1 M sodium hydroxide with 1 M sodium acetate. The gradient used was: 0% (B) from 0 min to 2 min, 20% (B) from 2 min to 35 min, 100% (B) from 35 min to 36.5 min, and 0% (B) from 37.5 min to 38.5 min. Detector temperature was maintained at 20 °C, data collection was at 2 Hz, and the Gold Standard PAD waveform (standard quadratic potential) was used.

Data acquisition, processing, and peak integration were performed using the Chromeleon™ version 7.1.3.2425 software (Thermo Scientific). Compounds were annotated based on available analytical standards. Glucose, fructose, sucrose, raffinose, 1-kestose, maltose, maltodextrin, nystose, and mixtures of inulin from chicory (DP2–60) and levan from Erwinia herbicola were purchased from Sigma-Aldrich, while 1,1,1-kestopentaose was obtained from Megazyme. Additional inulin- and neoseries-type fructans were isolated from onions and barley grain and analysed by MS. Fructan-related chromatographic peaks were identified based on fructanase digestion and mild acid hydrolysis. A total of 27 peaks were annotated (Table 1; Supplementary Table S2).

Table 1.

List of the 27 annotated metabolites

Peak no.	Compound	RT (min)	DP	Molecular family	Fructan-type
1	Glucose	3.43	1	Monosaccharide
2	Fructose	3.80	1	Monosaccharide
3	Melibiose	4.36	2	Disaccharide
4	Sucrose	5.95	2	Disaccharide
5	Raffinose	9.20	3	Raffinose family oligosaccharides
6	1-Kestose/6-kestose	9.74	3	Fructan	Inulin/levan
7	Maltose	10.47	2	Maltose-type oligosaccharides
8	Fructosylraffinose	10.73	4	Raffinose family oligosaccharides
9	B-type procyanidin 1	11.53	2	Flavonoids
10	6G-Kestose	12.13	3	Fructan	NS-inulin
11	Nystose	12.76	4	Fructan	Inulin
12	B-type procyanidin 2	13.08	2	Flavonoids
13	1&6-KT (bifurcose)/6G&1-KT (NS-DP4)	13.52	4	Fructan	Graminan, NS-inulin
14	Maltotriose	14.11	3	Maltose-type oligosaccharides
15	UK-Fructan1	14.68	–	Fructan	Unknown
16	1,1,1-Kestopentaose	15.68	5	Fructan	Inulin
17	6G,1-KP (NS-DP5)	16.38	5	Fructan	NS-inulin
18	6G&1-KP (NS-DP5)	16.61	5	Fructan	NS-inulin
19	Maltotetraose	17.51	4	Maltose-type oligosaccharides
20	Inulin-DP6	18.47	6	Fructan	Inulin
21	Neoseries-DP6	19.28	6	Fructan	NS-inulin
22	Inulin-DP7	21.10	7	Fructan	Inulin
23	Neoseries-DP7	22.73	7	Fructan	NS-inulin
24	Inulin-DP8	23.48	8	Fructan	Inulin
25	Neoseries-DP8	25.15	8	Fructan	NS-inulin
26	Inulin-DP9	25.94	9	Fructan	Inulin
27	Inulin-DP10	28.10	10	Fructan	Inulin

Peaks were annotated by comparison with analytical standards and isolated fractions from barley grain and onion bulb samples, as well as based on fructanase digestion and mild acid hydrolysis. Further details of compound annotation are provided in Supplementary Table S2.

Abbreviations: DP, degree of polymerization; KP, kestopentaose; KT, kestotetraose; NS, neoseries-type fructan; RT, retention time.

Fructanase assay and acid hydrolysis

Grain extract (30 µl) was diluted with water (335 µl) and treated with 10 µl of either water (control) or fructanase (1 U µl^–1, Megazyme fructan assay kit, Deltagen). The mixture was incubated at 40 °C for 2 h for complete hydrolysis, and heated to 90 °C for 5 min. After centrifugation, 25 µl of the mixture were injected to the chromatographic system and oligosaccharide profiles analysed as described above.

For mild acid hydrolysis, grain extract (10 µl) was either treated with water (250 µl, control) or hydrolysed with 250 µl of 0.2 M trifluoroacetic acid (TFA) for 1 h at 85 °C and the acid removed by evaporation. Specifically, samples were evaporated to dryness, treated with 50 µl of methanol, and then evaporated to dryness again. For HPAEC-PAD analysis, samples were dissolved in 375 µl of water and 25 µl was injected for every sample.

Oligosaccharide isolation

Soluble sugars were extracted from 300 mg of the ‘mixed sample’ of barley grain as described above, resulting in a 12 ml volume extract. Oligosaccharides were then fractionated by solid phase extraction (SPE) on graphitized carbon. In detail, 3 ml of extract were loaded per water-conditioned column (carbon SPE column, 38–125 µm, 150 mg, 4 ml, GRACE). The flowthrough was collected and the columns were each washed with 3 ml of water. The oligosaccharide fractions were eluted with a stepwise increasing concentration of organic solvent, namely 2, 10, 20, and 98% of acetonitrile, 3 ml for each. The resulting fractions were combined, and oligosaccharide profiles analysed as described above. As considerable amounts of fructans were still present in the flowthrough and water-wash fractions, these fractions were subjected to a repeat round of SPE fractionation as described and fractions were combined with those previously obtained. As fructans were mainly eluted with 10% acetonitrile, this fraction was further subjected to preparative HPLC on an Agilent 1200 series HPLC system operated with the OpenLAB software (Agilent) and coupled to evaporative light scattering detection (ELSD) with an Alltech^® ELSD 800 (Buchi). The total 10% acetonitrile fraction was evaporated to dryness and dissolved in 1700 µl of water. Portions of 50 µl were injected onto a Thermo Scientific™ Hypercarb Porous Graphitic Carbon LC Column (250 Å, 100 mm×4.6 mm, 5 µm; Thermo Fisher Scientific) at a flow rate of 1 ml min^–1. Fructan separation was performed with water (A) and 90% acetonitrile (B) using the following gradient: 5% (B) initially, 23% (B) from 0 min to 14 min, 45% (B) from 14 min to 15 min, 80% (B) from 15 min to 16 min, and 5% (B) from 17 min to 18 min. The column temperature was set to 20 °C. The column effluent was connected to a splitter directing one part to the ELSD and nine parts to the collection tube. ELSD settings were 3 bar nitrogen gas pressure, 70 °C, and laser attenuation 8, while the detector response in the OpenLAB software was set to 1000. Fractions of 30 s each were collected from 0 min to 18 min, combined from the repetitive separations, and oligosaccharide profiles analysed as described above. After comparison with barley grain extracts, fractions containing isolated fructans were further subjected to MS analysis.

Onion bulbs contain high amounts of inulin- and inulin-neoseries-type fructans with a DP range of ~3–13 (Yamamori et al., 2015), thus representing a suitable source for the isolation of reference compounds for fructan annotation. Soluble sugars were extracted from 115 mg of freeze-dried onion bulb tissue as described above, resulting in a 4.5 ml extract. Oligosaccharides were then isolated from this extract by means of SPE and preparative HPLC in the same way as described for the barley grain sample, and further subjected to MS analysis.

Mass spectrometry analysis

Samples were prepared 1:1 (v/v) using 2,5-dihydroxybenzoic acid as a matrix (diluted to 30 mg ml^–1 in 50% methanol and 0.2% TFA) and applied onto an MTP384 target plate polished steel BC (Bruker Daltonics). MS measurements were performed using an ultrafleXtreme MALDI time-of-flight (TOF)/TOF device run with the flexControl v3.4 software (Bruker Daltonics) in positive ionization mode following Peukert et al. (2014). Instrument calibration was performed with a polyethylene glycol (PEG) mixture (1:1 mixture of PEG200 and 600, diluted 1:300 in 30% v/v acetonitrile and 0.1% w/v TFA). The settings for fragmentation analyses were: m/z range 80–2000 and sample rate 0.63 GS s^–1 for MS spectra acquisition as well as m/z range 40–1000, sample rate 1.25 GS s^–1, and precursor ion selector range 2 Da. The laser power and resolution were adjusted as required. The acquired MS data were processed with flexAnalysis v3.4 (Bruker Daltonics).

Metabolic data analyses

Peak area entry means and variances with respective SDs were calculated in Excel 2007 (Microsoft) from the ‘mixed sample’-normalized integrated peak area values of the individual biological replicates for each two-row spring barley line and the check lines (Supplementary Table S3). Data from at least three biological replicates were available for 143 lines. For 10 lines (Agenda, Alliot, Appaloosa, Cellar, Drought, Goldie, Scarlett, Tankard, Tartan, and Turnberry), data from two biological replicates were available. The mean values for the two lines with just one entry (Calgary and Saana) were replaced by the only available data. Therefore, descriptive statistics could only be achieved for 153 genotypes.

Bonferroni outlier test was performed and pairwise correlations between the abundances of the 27 metabolites were revealed by applying the average linkage clustering method, based on Pearson correlation coefficients implemented in the MVApp (Julkowska et al., 2019; https://mvapp.kaust.edu.sa/). The metabolite abundances were analysed with the software package MATLAB (The MathWorks, Inc.) with a log-logistic distribution applied. The neural gas (NG) algorithm, implemented in MATLAB, was applied for cluster analysis following Kaspar-Schoenefeld et al. (2016) and Peukert et al. (2016). Analyses were performed for the biological replicates individually and the number of NG clusters was set to four.

GWAS

GWAS was carried out by combing the phenotypic data for the 154 barley accessions with genotypic data generated using the Barley 50K iSelect genotyping platform (Bayer et al., 2017). We focused on two-row spring barley accessions to reduce the confounding effects of population structure (Comadran et al., 2012) that could have been introduced by including other row types and growth habits (Darrier et al., 2019). Prior to analysis, any single nucleotide polymorphism (SNP) with a minimum allele frequency (MAF) of <0.05 was removed, which left 24 925 polymorphic markers for our analysis. Marker–trait association analysis was carried out using R 2.15.3 (http://www.R-project.org) and performed with a compressed mixed linear model (Zhang et al., 2010) implemented in the GAPIT R package (Lipka et al., 2012). Linkage disequilibrium (LD) was calculated across the genome between pairs of markers using a sliding window of 500 markers and a threshold of R²<0.2 using Tassel v 5 (Bradbury et al., 2007) to identify local blocks of LD, facilitating a more precise delimitation of QTL regions. We anchored regions of the genome containing markers that had passed the Benjamini–Hochberg threshold (P<0.05) as implemented in GAPIT to the barley physical map (Mascher et al., 2017) using marker positions provided in Bayer et al. (2017) and then expanded this region using local LD derived from genome-wide LD analysis as described above. Putative QTLs represented by <5 SNPs with –log10(P) values <3 were not considered to be robust given the marker density and extensive LD present in the barley genome (Mascher et al., 2017). The SNP with the highest logarithm of odds (LOD) score which passed the false discovery rate (FDR) threshold of –log10(P) 6.02 was used to represent a significant QTL. We investigated significantly associated regions using BARLEX (https://apex.ipk-gatersleben.de/apex/f?p=284:39) to identify putative candidate genes. Gene annotations refer to entries in the UniProt database (https://www.uniprot.org/uniprot/, June 2019). Unknown genes were searched against the non-redundant entries for plants (June 2019) in the NCBI database using the BLASTX 2.9.0+ software (https://blast.ncbi.nlm.nih.gov/Blast.cgi).

For compounds showing an association that passes the FDR calculated in GAPIT, ratios between these and all other compounds quantified were generated. The ratios were log transformed and then used to carry out further GWAS. The ‘p-gain’, defined as the ratio of the lowest P-value of the two individual metabolites and the P-value of the metabolite ratio (Petersen et al., 2012), was then calculated. A critical value for the p-gain was derived using B/(2×α), where α is the level of significance (0.05) and B the number of tested metabolite pairs. As we tested 52 pairs of compounds, our critical value threshold was 5.2×10².

Publicly available exome capture datasets (Mascher et al., 2017) were used to identify potential causal polymorphisms in candidate genes. We only considered non-synonymous SNPs with <10% missing data across the set of germplasm to be informative.

Gene transcript expression analyses of various developmental stages and tissues

Gene transcript data were obtained from previously published studies or public datasets as follows. The transcript abundance of genes of interest was measured in whole germinated grain (mean data from genotypes Navigator and Admiral) and isolated Navigator grain tissues from 0 to 96 hours after imbibition (hai). Aleurone tissues were divided into approximately thirds, with the proximal aleurone closest to the embryo (Betts et al., 2020).

Data for seedling tissues (germinated embryo, root, and shoot) were obtained from the Expression Atlas organ dataset (https://www.ebi.ac.uk/gxa/home).

Data from epidermal strips (4 weeks after sowing, W4), roots (W4), inflorescences-rachis (W5), inflorescences-lemma, inflorescences-lodicule, dissected inflorescences-palea (W6), inflorescence (10 mm), and internode, as well as for senescing leaves were obtained from BARLEX (https://apex.ipk-gatersleben.de/apex/f?p=284:10) (Colmsee et al., 2015).

A developing anther dataset was obtained from Barakate et al. (2020, Preprint) covering four anther stages (pre-meiosis, leptotene/zygotene, metaphase I to tetrad, and pachytene/diplotene) and two meiocyte stages (leptotene/zygotene and pachytene/diplotene). Raw expression data were mapped against the transcriptome of barley (https://webblast.ipk-gatersleben.de/barley_ibsc/downloads/; merging high-confidence and low-confidence transcripts as well as isoforms) using Salmon v14.0 (Patro et al., 2017).

RNA-sequencing data were obtained from developing pistils at Waddington (W) stages W8, W8.5, W9, W9.5, and W10 (Wilkinson, 2019) and are shown as a mean value for five genotypes including Golden Promise (1× replicate per stage; Aubert, 2018), Salka, Wren, Forum, and Gant (2× replicates per stage; Wilkinson et al., 2019). In addition, RNA-sequencing data from individual pistil tissues including the nucellus, integuments, ovary wall, embryo sac, egg apparatus and central cell, antipodal cells, and chalaza were analysed from the Sloop genotype (Aubert et al., 2018; Wilkinson et al., 2019).

Gene transcript expression data for whole developing grain (from 7 to 20 DAP) minus the embryo were generated by RNA-sequencing and are shown as a mean value from six genotypes including Sloop (1×replicate per time point; Aubert, 2018), Alabama, Pewter, Extract, Taphouse, and Hopper (1× replicate per time point; Aubert et al., 2018), while isolated developing grain tissues of interest including the pericarp, aleurone, subaleurone, and starchy endosperm were generated from medial sections at 7–25 DAP for the genotype Sloop (1× replicate per time point; Aubert et al., 2018).

Data obtained from BARLEX were available as fragments per kilobase, per million mapped reads (FPKM), while all other data represent (TPM).

Correlation analyses of gene transcript expression

Correlations among transcript abundance of fructan metabolism genes with other gene models from the QTL interval detected in the GWAS were evaluated for each of the RNA-sequencing datasets listed above, individually. Pairwise correlations between the gene transcript expression levels were revealed by applying the average linkage clustering method, based on Pearson correlation coefficients implemented in the MVApp (Julkowska et al., 2019).

Results

Grain oligosaccharide profiling revealed the abundance of fructans

HPAEC-PAD chromatograms of non-structural soluble carbohydrates from mature barley grain allowed separation of monosaccharides, disaccharides, and oligosaccharides with a DP <15 (Table 1; Supplementary Fig. S1). Among the latter, we identified two raffinose family oligosaccharides (RFOs) and three maltose-type oligosaccharides. Most compounds were found to be fructans (Supplementary Table S2), including levan-, inulin-, graminan-, and neoseries-inulin types. A high abundance of fructans with DP3 and DP4 was observed in mature grain extracts. Oligosaccharide profiles were obtained and evaluated from all 154 two-row spring barley lines and four checks (Supplementary Table S1) and integrated peak areas extracted for 27 compounds (Table 1). The resulting data matrix was used for analyses of oligosaccharide distribution, abundance variation, metabolite correlations, and GWAS.

Large variations detected in oligosaccharide profiles

The abundance of most compounds followed a log-logistic distribution. Only the two most abundant compounds, sucrose and raffinose, followed a normal distribution (Supplementary Fig. S2).

Oligosaccharide profiles were grouped separately for each of the biological replicates. Applying NG clustering to the data identified four statistically significant patterns of abundance (clusters). Each cluster can be interpreted as a prototypic abundance profile of the underlying metabolite values (Supplementary Fig. S3). They mainly differed in height of the normalized peak areas. Overlaps between cluster 1 and 4 and cluster 2 and 3 were detected. Samples in the latter clusters were characterized by significantly higher levels of sugar monomers and sucrose, lower fructosylraffinose, and higher overall fructan values compared with cluster 1 and 4 samples (Supplementary Fig. S3). Accordingly, we rationalized the four clusters into two profile groups (Fig. 1); cluster 1 and 4 forming profile group 1 (‘low’) and cluster 2 and 3 forming profile group 2 (‘high’). The largest peak in each sample was sucrose, whilst among the oligosaccharides, the highest values were detected for raffinose, the co-eluting fructans 1-kestose/6-kestose, nystose, and the co-eluting 1&6-kestotetraose (KT, bifurcose)/6G&1-KT (neoseries-DP4). Generally, a higher abundance of fructans with DP3 and DP4 was observed for all accessions, and differentiation between individuals was mainly attributed to the overall abundance of all fructan types.

Fig. 1.

Oligosaccharide profile groups as obtained by neural gas clustering. The mean profiles of the two major profile groups are shown. Value (y-axis) represents the peak area [nC×min] of the individual compounds (x-axis) listed in Table 1. Profile group 2 was characterized by higher levels of sugar monomers and sucrose, lower fructosylraffinose, and higher overall fructan values. Abbreviations: G, glucose; F, fructose; Me, melibiose; S, sucrose; R, raffinose; 1-K/6-K, 1-kestose/6-kestose; M, maltose; FR, fructosylraffinose; P1, procyanidin B1; 6G-K, 6G-kestose; N, nystose; P2, procyanidin B2; B-1&6-KT (bifurcose)/6G&1-KT (NS-DP4); Mtr, maltotriose; UK, unknown fructan; I-DP5, inulin DP5; NS-DP5, neoseries-DP5 (probably 6G;1-KP and 6G&1-KP); Mte, maltotetraose; I-DP6, inulin DP6; NS-DP6, neoseries-DP6; I-DP7, inulin DP7; NS-DP7, neoseries-DP7; I-DP8, inulin DP8; NS-DP8, neoseries-DP8; I-DP9, inulin DP9; I-DP10, inulin DP10.

We then assigned individual barley accessions to profile groups according to abundance profiles in each individual replicate (Supplementary Table S1). Accessions with only two biological replicates and mixed representation of their clusters in the two profile groups were assigned as ‘mixed’, as their oligosaccharide profile group was not distinct. In total 76, 77, and five accessions were assigned to the profile groups ‘low’, ‘high’, and ‘mixed’ respectively.

Significant metabolite correlations reflect underlying biochemical pathways

We have investigated metabolite correlations in order to prove the connection of the individual fructan types to related biosynthetic processes. Of the 349 pairwise correlations, 184 (52.72%) were highly significantly correlated (P<0.001), 204 (58.45%) moderately significantly (P<0.01), and 228 (65.33%) were just significant (P<0.05) (Supplementary Table S4). Several regions with highly correlated metabolites were identified in the results matrix, reflecting, in many cases, biochemical relationships (Fig. 2). The most significant positive correlations were observed between the various branched neoseries-type fructans, such as between neoseries-DP6, 6G&1-KP (neoseries-DP5), and neoseries-DP8, between 1&6-KT (bifurcose)/6G&1-KT (neoseries-DP4) and 6G&1-KP (neoseries-DP5), as well as between 6G,1-KP (neoseries-DP5), 6G&1-KP (neoseries-DP5), neoseries-DP6, and 1&6-KT (bifurcose)/6G&1 KT (neoseries-DP4). Similarly, significant positive correlations were observed between the linear inulin-type fructans, such as inulin-DP6 and inulin-DP7 and 1,1,1-kestopentaose (DP5). Notable significant positive correlations were also observed between the linear inulin-DP9 and the branched neoseries-DP7, as well as between nystose (DP4) and 6G&1-KP (neoseries-DP5). In contrast, inulin-DP6, inulin-DP7, and inulin-DP8 all showed slightly negative significant correlations with neoseries-DP7, which might hint at different enzyme kinetic particularities depending on the chain length. Significant positive correlations were detected between the monosaccharides and their related disaccharides, as well as between the maltose-type oligosaccharides. In contrast, glucose, fructose, and sucrose all showed weaker significant positive correlations with all fructans, except inulin-DP10. The most significant negative correlations were detected between fructosylraffinose and fructans, being highest with nystose (Fig. 2).

Fig. 2.

Correlation pattern among metabolites. Pairwise Pearson correlations are shown in a heat map representation, whereas metabolites are sorted according to correlation-based hierarchical cluster analysis. High positive correlations are represented by dark blue and negative ones by red circles, whereas the circle diameter is indicative for the strength of the correlation. X, not significantly correlated (P>0.05). The raw data are presented in Supplementary Table S4.

Differences in grain oligosaccharide profiles are genetically controlled in barley

GWAS for variation in mature grain oligosaccharides identified a single highly significant association for two compounds, neoseries-DP7 (LOD=8.65, P=2.25×10^–9) and inulin-DP9 (LOD=6.74, P=1.81×10^–7), with other less significant associations for both of these on chromosome 7H (Fig. 3A). Regression of these two traits showed a high level of correlation (R²=0.86, Fig. 2; Supplementary Table S4). Both QTLs on 7H overlapped, and the most significant marker from the analysis was the same, JHI-Hv50k-2016-438638 (Table 2A). This marker had adjusted P-values after FDR correction of P=0.00004 for neoseries-DP7 and P=0.003 for inulin-DP9. We anchored this QTL to the physical map (Mascher et al., 2017) which, based on local LD, spans 3.88 Mb from bases 174 327 (JHI-Hv50k-2016-435062) to 4 056 691 (JHI-Hv50k-2016-439312).

Table 2.

Significant GWA results for the two metabolites neoseries-DP7 and inulin-DP9

A. Information on the detected significant QTLs
Trait	Chromosome	Peak marker	Peak marker bp	MAF	LOD	QTL start and end bp
Neoseries-DP7	7H	JHI-Hv50k-2016-438638	3407292	0.14	8.65	174327..4056691
Inulin-DP9					6.74
B. Candidate gene models related to fructan biosynthesis
Gene model	HORVU7Hr1G000250.3	HORVU7Hr1G000260.2	HORVU7Hr1G000270.1	HORVU7Hr1G001040.6		HORVU7Hr1G001070.17
Gene model location	261065..262018	276441..280334	318543..322514	2257431..2260153		2423328..2427280
Description	Acid beta-fructofuranosidase, GH family 32 protein	Acid beta-fructofuranosidase, GH family 32 protein	Acid beta-fructofuranosidase, GH family 32 protein	Acid beta-fructofuranosidase, GH family 32 protein		Acid beta-fructofuranosidase, GH family 32 protein
PFAM	PF00251, PF08244, PF11837	PF00251, PF08244, PF11837	PF00251, PF08244	PF00251, PF08244		PF00251, PF08244
UniprotKB	A0A287VBC4	M0XA31	J7GIU6	M0YBF9		A0A287VC78
Similar proteins	J7GHS0	B6ECP1	M8C9V9	Q9AUH3		Q6PVN1
Similarity	0.9	1	0.9	1		0.9
Description	Fructan:fructan 1-fructosyltransferase	Sucrose:sucrose 1-fructosyltransferase	6(G)-fructosyltransferase	Sucrose:fructan 6-fructosyltransferase		Vacuolar invertase1
Species	Hordeum vulgare	Hordeum vulgare	Aegilops tauschii	Hordeum vulgare		Triticum monococcum
Abbreviation	1-FFT	1-SST	6G-FFT	6-SFT		VI-1

Abbreviations: DP, degree of polymerization; LOD, logarithm of odds; MAF, minimum allele frequency; QTL, quantitative trait locus.

Fig. 3.

Results of the GWA scan. (A) Manhattan plots are shown for the two fructans neoseries-DP7 and inulin-DP9. The –log₁₀ (P-value) is shown on the y-axis and the seven barley chromosomes are shown on the x-axis. Marker–trait association analysis was based on mean integrated peak areas. Integrated box plots show information for the top SNP (JHI-Hv50k-2016-438638). The variable alleles found (A and G) are shown on the x-axis. The y-axis shows the median of the metabolite amount for all lines with the respective allele variant; the width of the box is indicative of the number of lines with this particular allele. The false discovery rate significance threshold= –log 10(P) 6.02. (B) Assembling of the annotated gene models from the significant QTL region according to functional categories as obtained from the UniProt database (https://www.uniprot.org/uniprot/, June 2019). Numbers represent the count of gene models with the respective functional annotation.

In total, 194 gene models were detected within the QTLs, of which 65 are unannotated (Supplementary Table S5). The highest number of annotated gene models was involved in protein modification (32) and degradation (17) (Fig. 3B), with others involved in transcription/translation (27), lipid/sterol/terpenoid metabolism (12), transcription factors (TFs) (11), or carbohydrate metabolism (11). Among the latter category, we identified five candidates that could influence fructan content (Table 2B). These included HORVU7Hr1G000250.3, HORVU7Hr1G000260.2, and HORVU7Hr1G001040.6, which were genetically similar or identical to 1-FFT, 1-SST, and 6-SFT from Hordeum vulgare, respectively. Of the two others, HORVU7Hr1G000270.1 was similar to 6(G)-fructosyltransferase (6G-FFT) from Aegilops tauschii, and HORVU7Hr1G001070.17 to vacuolar invertase1 (VI-1) from Triticum monococcum.

To further explore relationships between compounds, we used a hypothesis-free analysis of metabolite ratios in a GWAS. This analysis generates a ‘p-gain’ statistic which is calculated from the significance of increases in –log10(P) values of the metabolite ratios compared with an estimated threshold derived from the P-values obtained in GWAS of the individual compounds (Petersen et al., 2012). Using the log-transformed ratios between neoseries-DP7 and inulin-DP9 with all other compounds, 17 pairs of compounds correlated with a QTL that passed the FDR threshold of –log 10(P) 6.02 in the same region of chromosome 7H as neoseries-DP7 and inulin-DP9 alone. The ratios of neoseries-DP7:inulin-DP10, inulin-DP9:neoseries-DP8, and inulin-DP9:inulin-DP10 passed the p-gain threshold of 5.2×10² (P<0.05) for markers with an MAF of >10% with 1.97×10⁵, 4.08×10⁹, and 3.71×10³, as well as 1.93×10⁵ and 6.43×10⁵, respectively (Table 3; Supplementary Fig. S4), indicating metabolic links between these compounds. The QTLs on 7H overlapped for all ratios. Significant markers identified were SCRI_RS_8079, JHI-Hv50k-2016-435510, and JHI-Hv50k-2016-438638, the latter being the same as identified with the metabolite concentrations for neoseries-DP7 and inulin-DP9 alone (Table 2A). GWAS for the ratio neoseries-DP7:inulin-DP9 did not identify any significant associations (Table 3).

Table 3.

Significant GWA results for the ratios of neoseries-DP7 and inulin-DP9 with all other metabolites

Ratio combination	Peak marker	Peak marker bp	MAF	LOD	P	P for non sig trait	FDR adjusted P-value	FDR adjusted P-value for non sig trait	p-gain	FDR-adjusted p-gain
Neoseries-DP7:fructosylraffinose	JHI-Hv50k-2016-438638	3407292	0.14	7.0	1.08E-07	6.26E-04	2.24E-03	6.76E-01	2.08E-02	2.08E-02
	JHI-Hv50k-2016-438980	3689408
Neoseries-DP7:B-type procyanidin1				6.9	1.39E-07	4.73E-05	2.87E-03	1.00E+00	1.62E-02	1.62E-02
Neoseries-DP7:1,1,1-kestopentaose				9.1	8.62E-10	8.22E-05	1.79E-05	2.84E-01	2.61E+00	2.61E+00
Neoseries-DP7:6G,1-KP				9.0	8.96E-10	9.06E-01	1.86E-05	1.00E+00	2.51E+00	2.51E+00
Neoseries-DP7:inulin-DP6				8.4	3.88E-09	9.78E-04	8.04E-05	1.00E+00	5.80E-01	5.80E-01
Neoseries-DP7:neoseries-DP6				7.9	1.28E-08	1.51E-01	2.66E-04	1.00E+00	1.75E-01	1.75E-01
Neoseries-DP7:inulin-DP7				9.1	8.45E-10	3.17E-03	1.75E-05	1.00E+00	2.66E+00	2.66E+00
Neoseries-DP7:neoseries-DP8				10.8	1.64E-11	7.89E-02	2.14E-07	1.00E+00	1.37E+02	2.17E+02
Neoseries-DP7:inulin-DP10	SCRI_RS_8079	2325562	0.11	7.4	4.24E-08	3.46E-01	1.72E-03	1.00E+00	1.97E+05*	5.80E+02
Inulin-DP9:1-kestose	JHI-Hv50k-2016-438638	3407292	0.14	8.0	1.03E-08	1.72E-02	2.14E-04	1.00E+00	1.75E+01	1.75E+01
	JHI-Hv50k-2016-438980	3689408
Inulin-DP9:1,1,1-kestopentaose				9.4	4.38E-10	8.22E-05	9.08E-06	2.84E-01	4.14E+02	4.14E+02
Inulin-DP9:6G,1-KP				7.3	6.29E-06	9.06E-01	9.40E-04	1.00E+00	2.88E-02	4.00E+00
Inulin-DP9:inulin-DP6				8.4	4.03E-09	9.78E-04	8.35E-05	1.00E+00	4.50E+01	4.50E+01
Inulin-DP9:neoseries-DP6				8.0	1.04E-08	1.51E-01	2.15E-04	1.00E+00	1.75E+01	1.75E+01
Inulin-DP9:inulin-DP7				9.0	9.04E-10	3.17E-03	1.87E-05	1.00E+00	2.01E+02	2.01E+02
Inulin-DP9:neoseries-DP8	JHI-Hv50k-2016-435510	280043	0.07	10.7	1.93E-11	3.71E-04	4.01E-07	1.00E+00	4.08E+09*	2.02E+06
	JHI-Hv50k-2016-437376	2350955
	JHI-Hv50k-2016-438638	3407292	0.14	10.3	4.89E-11	7.89E-02	5.06E-07	1.00E+00	3.71E+03*	7.42E+03
	JHI-Hv50k-2016-438980	3689408
Inulin-DP9:inulin-DP10	JHI-Hv50k-2016-435510	280043	0.07	6.9	1.36E-07	3.21E-01	1.93E-03	1.00E+00	1.93E+04*	4.18E+02
	JHI-Hv50k-2016-437376	2350955
	SCRI_RS_8079	2325562	0.11	6.9	1.40E-07	3.46E-01	1.93E-03	1.00E+00	6.43E+05*	5.17E+02
Neoseries-DP7:inulin-DP9	JHI-Hv50k-2016–438638	3407292	0.14	1.6	2.31E-02	1.80E-07	1.00E+00	3.74E-03	9.75E-08	4.66E-05
	JHI-Hv50k-2016–438980	3689408

*Significant results passing the p-gain threshold of 5.2×10², which are also highlighted in light grey. Corresponding Manhattan and box plots are shown in Supplementary Fig. S4

Abbreviations: DP, degree of polymerization; FDR, false discovery rate; LOD, logarithm of odds; MAF, minimum allele frequency; p-gain, ratio of the lowest P-value of the two individual metabolites and the P-value of the metabolite ratio; P-values and FDR-adjusted P-values from the initial GWAS for traits that did not identify significant associations but when included as a ratio do identify significant associations are included in columns ‘P for non sig trait’ and ‘FDR adjusted P-value for non sig trait’.

Evaluation of exome capture data revealed several non-synonymous SNPs

Mascher et al. (2017) presented exome capture data for 25 of the genotypes included in our study, which we evaluated to identify putative casual SNPs for our five regional candidates involved in fructan biosynthesis. Eight non-synonymous SNPs in 1-FFT, three in VI-1, two in 6G-FFT, and one in 1-SST were identified (Supplementary Table S6). All identified SNPs are located within functional protein-coding regions of the genes (Fig. 4). Changes in just one out of 25 genotypes were observed for three markers among the eight SNPs detected in 1-FFT. The other five SNPs in 1-FFT represent changes from methionine to leucine (position 7H_ 262685), alanine to threonine (7H_263547 and 7H_263700), isoleucine to threonine (7H_264127), and leucine to isoleucine (7H_264198). They showed significant effects (P<0.05) on 1-kestose and several neoseries-type fructans (Supplementary Fig. S5). The two SNPs in 6G-FFT were in LD and represent changes from glycine to glutamic acid (7H_321608), and alanine to threonine (7H_319284). Notably, they have a significant effect on 1-kestose and several inulin-type fructans (Supplementary Fig. S5). However, the SNP in 1-SST, representing a change from threonine to isoleucine (7H_279526), as well as the three SNPs in VI-1, representing changes from glutamic acid to aspartic acid (7H_2423349), tryptophan to arginine (7H_2425560), and arginine to cysteine (7H_2425578), did not show a significant effect on either trait. For 6-SFT no SNP was identified.

Fig. 4.

Location of non-synonymous SNPs in fructan biosynthesis candidate genes. All significant causal SNPs (position assigned in bold red) are located within functional protein-coding regions of the genes (black regions of the transcripts). All identified SNPs are presented in Supplementary Table S6, and box plots for significant effects are shown in Supplementary Fig. S5. For 6-SFT no SNP was identified.

Fructan biosynthesis genes show developmental stage- and tissue-specific expression patterns

We compared expression patterns of the five candidate genes (Table 2) and three known fructan hydrolyase-encoding genes in various tissues across barley plant development (Fig. 5). In the vegetative phase, the highest expression for 1-SST (facilitating the biosynthesis of 1-kestose, the precursor for production of inulin- and graminan-type fructans) was observed during early germination in embryo and all seedling tissues (Fig. 5A). During the reproductive phase, 1-SST is expressed in all vegetative tissues, with peak expression in the leaf epidermis, as well as in all reproductive tissues and stages, with a pronounced peak of expression in the ovary wall, the embryo sac (ES), the egg apparatus and central cell (EC+CC), and the antipodal cells (ANT) during late pistil development (stages W8–W10, Fig. 5B). In the grain development phase, 1-SST expression is highest in the early stages (7–9 DAP) in maternal grain tissues [pericarp, aleurone, subaleurone/outer starchy endosperm (SA)] while decreasing during the storage stage (from 11 DAP onwards) in all tissues (Fig. 5C). 1-FFT (mediating the biosynthesis of inulin-type fructans) showed tight co-expression with 1-SST during early germination in embryo and all seedling tissues (Fig. 5A) as well as all vegetative tissues (Fig. 5B), while in reproductive (Fig. 5B) and grain tissues (Fig. 5C) much lower expression levels were observed. In contrast, 6-SFT (mediating the biosynthesis of graminan-type fructans) showed very tight co-expression with 1-SST during meiosis and pistil development (Fig. 5B) as well as at early grain development (Fig. 5C). During germination, 6-SFT expression was extremely low while it was observed to be moderate in seedling (Fig. 5A) and all vegetative tissues (Fig. 5B). Notably, expression of 6G-FFT (mediating the biosynthesis of neoseries-type fructans) was restricted to the outer grain tissues (see aleurone tissues in Fig. 5A and aleurone, pericarp, and endosperm tissues in Fig. 5C) during late grain development (from 11 DAP onwards). VI-1, with an as yet unknown function, showed low expression levels in germinated grain tissues, all vegetative tissues, in pericarp at late grain development, and in senescing leaf, while higher levels were notable during late pistil development (Fig. 5). Among the fructan hydrolyases, 1-FEH (HORVU6Hr1g011260) and 6-FEH (HORVU2Hr1G109120) seem to be involved in balancing fructan biosynthesis, with 1-FEH tightly co-expressed with 1-SST in all tissues and stages, and pronounced 6-FEH expression during the reproductive phase in all tissues and in the pericarp at late grain development. In contrast, only marginal expression levels were observed for 6-FEH/CWI2 (HORVU2Hr1G118820) (Fig. 5).

Fig. 5.

Transcript expression of the fructan biosynthesis genes in barley across plant development and in various tissues. Datasets included in the analysis are detailed in the Materials and methods section. Genes included comprise the five fructan biosynthesis genes from the detected significant QTL region (1-FFT, HORVU7Hr1G000250; 1-SST, HORVU7Hr1G000260; 6G-FFT, HORVU7Hr1G000270; 6-SFT, HORVU7Hr1G001040; and VI-1, HORVU7Hr1G001070) as well as three known fructan hydrolyase-encoding genes (6-FEH, HORVU2Hr1G109120; 6-FEH/CWI2, HORVU2Hr1G118820; and 1-FEH, HORVU6Hr1g011260). Expression levels are colour coded, whereas different scales were used for TPM and FPKM values as indicated in the key. (A) Expression levels in the early vegetative phase for whole germinated grain tissues and for isolated germinated grain tissues from 0 to 96 hours after imbibition (hai). Also, expression levels in seedling are shown for germinated embryo (GE), root, and shoot. (B) Data from the reproductive phase. Vegetative tissues included are: EPI, epidermal strips (4 weeks after planting, W4); ROO2, roots (W4); RAC, inflorescences, rachis (W5); LEM, inflorescences, lemma (W6); LOD, inflorescences, lodicule (W6); PAL, dissected inflorescences, palea (W6); INF2, inflorescence (10 mm); and NOD, internode. Meiosis stages included are: A.Pre, pre-meiosis anthers; A.LepZyg, leptotene/zygotene anthers; M.Lep/Zyg, leptotene/zygotene meiocytes; A.MetTet, metaphaseI-tetrad anthers; A.PacDip, pachytene/diplotene anthers; M.PacDip, pachytene/diplotene meiocytes. Waddington (W) stages for pistil development are: W8, W8.5, W9, W9.5, and W10. Isolated pistil tissues are: nucellus (including nucellus and embryo sac for W8.5); integument; ovary wall; ES, embryo sac; EC+CC, egg apparatus and central cell; ANT, antipodal cells; and chalaza. (C) Data from the grain development phase for whole developing grain (from 7 to 20 days after pollination, DAP) and for isolated developing grain tissues from 7 to 25 DAP. Abbreviations: SA, subaleurone/outer starchy endosperm; SE, starchy endosperm/inner starchy endosperm. Also, data from senescing leaf (SN) are presented. Other abbreviations are: At, anthesis; CI, collar initiation; Em, emergence; Hv, harvest; and Sw, sowing.

Fructan biosynthesis genes show differential co-expression patterns in developing barley grain

Besides the five fructan biosynthesis genes, the association of differential oligosaccharide profiles with other candidates in the identified genomic region may be possible. We hypothesized similar expression patterns for fructan biosynthesis genes and other candidates influencing the fructan levels in barley. Therefore, transcript expression levels were evaluated for all gene models within the QTL interval, and co-expression of genes was assessed individually within the developmental phases and tissues (Supplementary Table S7). We have focused on developing barley grain and significant correlations for the expression of fructan metabolism genes with each other and with TFs (Table 4).

Table 4.

Significant correlations for the expression of fructan metabolism genes with each other and with potential regulatory gene models from the detected QTL interval for developing barley grain

Gene name	Annotation	Molecular function	1-FFT	1-SST	6G-FFT	6-SFT	VI1	1-FEH	6-FEH
HORVU2Hr1G109120.2	Fructan 6-exohydrolase (6-FEH)	Carbohydrate metabolism			0.81			0.62
HORVU6Hr1g011260.19	Fructan 1-exohydrolase (1-FEH)	Carbohydrate metabolism							0.62
HORVU7Hr1G000250.3	Fructan-fructan 1-fructosyltransferase (1-FFT)	Carbohydrate metabolism		0.92	–0.97	0.88
HORVU7Hr1G000260.2	Sucrose:sucrose 1-fructosyltransferase (1-SST)	Carbohydrate metabolism	0.92		–0.84	0.99
HORVU7Hr1G000270.1	6(G)-fructosyltransferase (6G-FFT)	Carbohydrate metabolism	–0.97	–0.84		–0.80			0.81
HORVU7Hr1G000820.1	WD_REPEATS_REGION domain-containing protein	Transcription factor	0.95	0.99	–0.87	0.98
HORVU7Hr1G000910.1	NAC domain-containing protein	Transcription factor						0.81
HORVU7Hr1G001040.6	Sucrose:fructan 6-fructosyltransferase (6-SFT)	Carbohydrate metabolism	0.73	0.98
HORVU7Hr1G001050.1	AP2/ERF domain-containing protein	Transcription factor						–0.83
HORVU7Hr1G001070.17	Vacuolar invertase1 (VI-1)	Carbohydrate metabolism
HORVU7Hr1G001120.1	Protein ALWAYS EARLY 3	Transcription factor	0.80	0.92		0.95	0.78
HORVU7Hr1G001300.3	GRAS domain-containing protein (Scarecrow-like protein 6)	Transcription factor	0.68	0.84		0.88
HORVU7Hr1G001310.1	Scarecrow-like protein 22	Transcription factor	0.99	0.96	–0.95	0.94	0.80
HORVU7Hr1G001320.1	GRAS domain-containing protein	Transcription factor				0.80
HORVU7Hr1G001830.3	HTH Myb-type domain-containing protein	Transcription factor	0.88		–0.91			0.68	0.88

Positive correlations are shown in blue and negative correlations in orange, whereas the colour code is indicative for the strength of the correlation (the darker, the stronger). The significance threshold was P>0.05. Numbered boxes without formatting indicate values just above significance (P<0.055). Bold framed boxes indicate correlations which were detected in various developing grain datasets. Datasets correspond to those presented in Fig. 5, and methods are detailed in the Materials and Methods. The raw data are presented in Supplementary Table S8.

Highly positive correlations among the fructan metabolism genes were observed between 1-FFT, 1-SST, and 6-SFT; between 6G-FFT and 6-FEH; and for 1-FEH with 1-FFT, and 6-FEH. Notable negative correlations were observed for 6G-FFT with 1-FFT, and 1-SST.

Co-expression patterns with TFs were highly similar for 1-FFT, 1-SST, and 6-SFT in the developing grain. Notable positive correlations for 1-FFT, 1-SST, and 6-SFT expression were identified with the WD_REPEATS_REGION domain-containing protein (HORVU7Hr1G000820.1), the ALWAYS EARLY 3 (HORVU7Hr1G001120.1), and the two scarecrow-like protein genes (HORVU7Hr1G001300.3 and HORVU7Hr1G001310.1). In contrast, 6G-FFT showed significant negative correlations with the WD_REPEATS_REGION domain-containing protein (HORVU7Hr1G000820.1), scarecrow-like protein 22 (HORVU7Hr1G001310.1), and HTH myb-type domain-containing protein (HORVU7Hr1G001830.3). Significant positive correlations for VI-1 were observed in developing grain with the protein ALWAYS EARLY 3 (HORVU7Hr1G001120.1) and scarecrow-like protein 22 (HORVU7Hr1G001310.1), as observed for 1-FFT, 1-SST, and 6-SFT. 1-FEH showed co-expression patterns partly like those observed for 1-FFT. Besides the positive correlation with a myb-type transcription factor (HORVU7Hr1G001830.3), a strong negative correlation with the AP2/ERF domain-containing protein (HORVU7Hr1G001050.1) and a positive interaction (not significant) with a NAC domain-containing protein gene (HORVU7Hr1G000910.1) were identified. The highest positive correlations were noted for 6-FEH with 6G-FFT and the HTH myb-type domain-containing protein (HORVU7Hr1G001830.3), which in contrast was negatively associated with 6G-FFT (Table 4).

Additional and partly different patterns for the co-expression of fructan metabolism genes with other genes were observed across a range of developmental phases and tissues (Supplementary Table S7).

Discussion

Neoseries-type fructans are abundant in mature barley grain

Profiling of DP3–DP10 oligosaccharides revealed the abundance of 6G-kestose and higher DP neoseries-type fructans in mature barley grain for the first time (Table 1; Supplementary Fig. S1; Supplementary Table S2). While the presence of 6G-kestose has been reported in wheat and barley grain (Nilsson et al., 1986; Henry and Saini, 1989), higher DP variants of this fructan type have not been previously identified. Recent studies revealed the presence of neoseries-type fructans in oat, rye, spelt, and wheat flour (Verspreet et al., 2015a, 2017), but claimed its absence in barley (Verspreet et al., 2017). These contrasting observations may be explained by both the plant materials used (flour versus whole grain) and the technical constraints of fructan profiling. Noticeable accumulation of unidentified higher DP fructans was reported for the outer pericarp of developing wheat grain (Schnyder et al., 1993). Thus, utilizing whole grain here may have facilitated the detection of neoseries-type fructans in barley, probably accumulating in outer grain parts, as suggested by expression analysis for 6G-FFT (Fig. 5). Additionally, electronic properties of PAD, typically used for fructan profiling, require higher concentrations for the detection of higher molecular weight fructans (Rocklin and Pohl, 1983), resulting in a pronounced log-logistic distribution for those compounds that is also observed in our study (Supplementary Fig. S2) and which may have led to the under-representation of fructans with >DP4 in other studies. In the future, comprehensive grain fructan profiling could be improved by employing recently established LC-MS methodologies, as reviewed in Matros et al. (2019).

Barley accessions group according to their oligosaccharide accumulation patterns

Reported genotypic variation in grain fructan content ranges from 0.9% to 4.2% of dry matter (DM) among 20 barley breeding lines (Nemeth et al., 2014) and from 1.1% to 1.6% of DM among seven barley cultivars (De Arcangelis et al., 2019). These results correspond well to the variability for total fructan values (0.02–1.94% of grain DW) determined here (Supplementary Table S9). When discriminating between different chain length, Henry and Saini (1989) measured varying amounts of FOS with 0.26% (DP3), 0.2% (DP4), 0.03% (DP5), and 0.23% of DM (>DP5) in mature barley grains, which was confirmed by results from Jenkins et al. (2011). In our study, the lowest abundance range showed FOS with DP5 (traces to 0.16% of DM), while FOS with DP3 and DP4 ranged from 0.02% to 0.53% and traces to 0.44% of DM, respectively. Nemeth et al. (2014) observed a positive correlation between fructan values and the content of long chain fructans (>DP9, r=0.54, P=0.021). However, such an association could not be found in our dataset. Clustering of the oligosaccharide profiles from the 154 lines revealed two major profile groups, one each of higher and lower sugar values (Fig. 1). We detected significant positive correlations between biosynthetically closely related metabolites (e.g. within and between the different fructan types) with negative associations for antagonistic compounds (e.g. fructosylraffinose with all fructans, sugar monomers, and dimers, Fig. 2). Co-occurrence of fructans and RFOs has been reported for many plant species including wheat (Haskå et al., 2008) and barley (Henry, 1988), with the proposal that strong RFO and fructan accumulation do not occur together in a single plant species (Van den Ende, 2013). Notably, fructosylraffinose was only speculated to occur in barley (Cerning and Guilbot, 1973), while its presence was described in wheat decades ago (White and Secor, 1953; Saunders, 1971). In our study, fructosylraffinose was clearly identified with medium abundance (Supplementary Fig. S1C).

Differences in oligosaccharide profiles are genetically controlled in barley

A significant QTL on chromosome 7H affecting barley grain fructan levels was identified (Fig. 3A) and five genes involved in fructan metabolism were detected in this region (Tables 2, 3; Supplementary Table S5). We increased the power of GWAS by analysing metabolite ratios using the novel p-gain approach. As the p-gain passed an appropriate threshold (defined by the data), using ratios provided more information about the traits, and the genomic locus underlying them, than looking at the traits individually. Using ratios reduces background ‘noise’ in datasets, increasing statistical power to detect significant associations between traits and genomic loci (Petersen et al., 2012). Previous studies have demonstrated that including ratios between pairs of traits can strengthen associations identified and uncover novel information about biochemical pathways (Gieger et al., 2008; Illig et al., 2010; Suhre et al., 2011). Thus, ratio-GWAS represents an innovative approach for the discovery of new biologically meaningful associations in plants, as we could also demonstrate for the linked oligosaccharide pathways described here. However, when we used the neoseries-DP7:inulin-DP9 ratio in the GWAS, we did not identify an association on 7H, indicating less information provided by the ratio than the individual values. This may relate to the high positive correlation of these two compounds across the barley lines (R²=0.86, Fig. 2; Supplementary Table S4) and the close genomic location of the related fructan biosynthesis genes (Table 2). In contrast, for the ratios inulin-DP7:inulin-DP10 (R²=0.011, P=0.80), inulin-DP9:neoseries-DP8 (R²=0.56, P<0.05), and inulin-DP9:inulin-DP10 (R²=0.028, P=0.54), a significant QTL was identified. This points towards a stronger association between the identified genomic locus and the molecular weight of the fructans than with the fructan structure.

In wheat, two loci for differential total fructan contents in grain were identified on chromosomes 7A and 6D, which did not show significant interactions (Huynh et al., 2008b). Subsequent physical mapping provided indications for clustering of fructan biosynthesis genes in the genomes of both dicots and monocots (Huynh et al., 2012). For wheat and barley, the formation of a functional cluster was shown containing 1-SST (provided are the IDs of the most probable barley gene product; J7GM45_HORVV), 1-FFT (J7GHS0_HORVV), and 6-SFT (Q96466_HORVU) (Huynh et al., 2012), which were also identified here. Additionally, the authors found two vacuolar invertase (VI) genes (J7GIU6_HORVV and J7GR98_HORVV) in this cluster, of which we identified one, which is similar to 6G-FFT (J7GIU6_HORVV). The identification of 6G-FFT matches the detection of neoseries-type fructans in our study, which we reported for the first time in barley grain. Among the five candidates we identified was also a gene coding for an uncharacterized gene product (M0X3V0_HORVV) which is similar to a VI-1 from T. monococcum (Q6PVN1_TRIMO) that has not been described or annotated in barley before.

The evaluation of exome capture data (Mascher et al., 2017) led to the identification of several significant SNPs in the five fructan biosynthesis genes. SNPs in 1-FFT were associated with grain neoseries-type fructan content, while SNPs in 6G-FFT were associated with inulin-type fructan content (Fig. 4; Supplementary Table S6). Ideally, the influence of these SNPs would be validated in the complete set of germplasm used to quantify fructan content. This analysis would be likely to reveal additional SNPs that have not been identified in this subset of lines.

Developmental- and tissue-specific nature of barley fructan biosynthesis

Throughout plant development, 1-SST, 1-FFT, and 1-FEH showed strong co-expression, starting in embryo tissue during germination, accompanied later by 6-SFT expression in root, leaf, and stem, probably leading to the biosynthesis of inulin- and graminan-type fructans in those tissues until senescence (Fig. 5). These observations matched the consensus of inulin- and graminan-type fructans being the predominant polymers in barley tissues (Pollock and Cairns, 1991; Bonnett et al., 1997; Huynh et al., 2008a). Accordingly, 1-SST, 1-FFT, 6-SFT, and 1-FEH are the best studied fructan biosynthesis genes (Duchateau et al., 1995; Henson, 2000; Lüscher et al., 2000; Huynh et al., 2012). A key role was assigned to 1-SST (Wagner et al., 1983), and correlated transcription and activity was reported for 1-SST and 6-SFT in barley leaves (Nagaraj et al., 2004). A role for fructans as a temporal carbohydrate reserve has been widely accepted in vegetative tissues and roots (Pollock et al., 1996; Vijn and Smeekens, 1999; Housley, 2000), and can be assumed for the inulin- and graminan-type fructans in barley.

We observed co-expression of 1-SST, 6-SFT, 6-FEH, and 1-FEH in reproductive tissues, with a pronounced peak during late pistil development in ovary tissues (Fig. 5B) probably leading to specific accumulation of graminan-type fructans. In Campanula rapunculoides, the largest inulin-type fructan concentrations were found in petals and ovaries (Vergauwen et al., 2000). Based on the observation that petals in daylily (Hemerocallis) (Bieleski, 1993) and C. rapunculoides (Vergauwen et al., 2000) and leaves of Phippsia algida (Solhaug and Aares, 1994) rapidly degrade fructans upon flower opening, a role for them in flower expansion was suggested. However, the function of the different fructan types accumulating in C. rapunculoides ovary (inulin-type) and barley ovary tissues (graminan-type) remains unresolved at present. Also, the newly identified VI-1 showed peak expression specific to ovary tissues at late pistil development, while its function in fructan biosynthesis remains unclear.

In accordance with the detection of neoseries-type fructans and the identification of 6G-FFT in the significant QTL region, we observed 6G-FFT expression in barley grain (Fig. 5A, C). Notably, its expression was restricted to developing grain from 11 DAP onwards and confined to the outer endosperm and maternal tissues. Reports on this fructan type in developing grains of other cereals do not exist to our knowledge. Accumulation of neoseries-type fructans in the aleurone of mature grain may be related to favourable structural characteristics when compared with inulins and graminans and with the function of this tissue during germination. Some reports showed that fructan branching architecture is critical to physicochemical properties, such as water solubility or formation of aggregates at high concentration (Eigner et al., 1988; Wolff et al., 2000; Ponce et al., 2008). The more compact shape of neoseries-type fructans would allow higher concentrations to be stored in the desiccated aleurone. Better water solubility and pH stability of neoseries-type fructans would be advantageous during germination, when the aleurone hydrates and enzymes must be quickly activated and reach their substrates. Proving these hypotheses will require the comparative evaluation of physicochemical properties of neoseries-inulin- and graminan-type fructans in the future. Indeed, several studies have provided strong evidence for a positive relationship between enhanced fructan concentrations and better malting characteristics in barley varieties (Smith et al., 1980; Cozzolino et al., 2016, and references therein)

Potential regulators of barley grain fructan biosynthesis

Despite the increasing evidence of tissue specificity, there is limited knowledge of how fructan metabolism is orchestrated to adjust the storage and use of photosynthates during grain development. Within the significant QTL interval, we found several genes differentially co-expressed with the various fructan biosynthesis genes in developing grain. Among them were several TFs (Table 4; Supplementary Table S7; Fig. 6).

Fig. 6.

Fructan types, suggested biosynthesis routes, and potential regulators in developing barley grain. Specific spatiotemporal biosynthesis of oligofructans was observed for barley grains. Continuous bold arrows illustrate the major route of biosynthesis during the pre-storage phase (until 14 DAP) and the dashed bold arrows indicate the major route during the storage phase (until 20 DAP). During the pre-storage phase, high transcript levels for 1-SST and 6-SFT were observed for the endosperm leading to an accumulation of 6-kestose and bifurcose. With transition to the storage phase, a transcriptional switch was observed, resulting in high transcript levels of 1-SST in the nucellar projection (NP). 1-FFT was found to be exclusively expressed in the NP during the storage phase. Induction of the inulin-type fructan biosynthesis pathway led to high amounts of 1-kestose and nystose accumulating in the endosperm cavity (Peukert et al., 2014). The dotted bold arrow illustrates the major route of biosynthesis during the late storage phase with 6G-FFT transcripts detected in the outer endosperm (from 30 DAP onwards, Fig. 5), which matched the detection of neoseries-type oligofructans in mature barley grains (Fig. 1). Transcription factors (TFs) showing significant correlation of transcript expression pattern in developing grain with 1-FFT, 1-SST, and 6-SFT (positive), 6G-FFT (negative), and VI-1 (positive) are listed in the inserted text box. Inulin-neoseries represents linear fructans with β(2,1)- and β(2,6)-linked fructosyl units at the glucose (1F, 6G-di-β-d-fructofuranosylsucrose is shown; m=1, n=1), graminan-type represents branched fructans with β(2,1)- and β(2,6)-linked fructosyl units (bifurcose is shown; m=1, n=1), inulin-type illustrates linear fructans with β(2,1)-linked fructosyl units (1-kestose is shown; n=1), and levan-type shows linear fructans with β(2,6)-linked fructosyl units (6-kestose is shown; n=1). The arrows indicate the direction of further polymerization. Abbreviations: 1-FFT, fructan:fructan 1-fructosyltransferase; 1-SST, sucrose:sucrose 1-fructosyl-transferase; 6-SFT, sucrose:fructan 6-fructosyltransferase; 6G-FFT, fructan:fructan 6G-fructosyltransferase; VI-1, vacuolar invertase 1 (unknown role in developing grain).

It is generally agreed that initiation of fructan biosynthesis is triggered by an organ-specific sucrose threshold (Lu et al., 2002; Jin et al., 2017). Also, several molecular components in sucrose-mediated induction of plant fructan biosynthesis, such as protein phosphatases and kinases (Noël et al., 2001), second messenger Ca²⁺ (Martínez-Noël et al., 2006), small GTPases and phosphatidylinositol 3-kinase (Ritsema et al., 2009), as well as the plant hormones abscisic acid and auxin (Valluru, 2015), were shown to be required for activation of fructosyltransferase genes. In addition, the positive metabolic connection between fructan and (1,3;1,4)-β-glucan biosynthesis and the negative association of both pathways with starch biosynthesis has been shown previously (Lim et al., 2019). However, despite the presence of several kinases and β-glucanases within our detected genomic region, no hint at the regulatory involvement was found for any other of the aforementioned molecular components.

An opposing sugar sensing system was recently identified in barley, whereby a single gene on chromosome 2H encodes two functionally distinct TF variants [SUSIBA (sugar signaling in barley) 1 and 2], which respond differently to sucrose concentrations (Jin et al., 2017). However, no distinction was made between different tissues and fructan types, and it remains unclear if this system coordinates fructan and starch biosynthesis in general.

In wheat, TaMYB13, a R2R3-MYB TF, was described as a transcriptional activator of fructan biosynthesis (Xue et al., 2011; Kooiker et al., 2013). In the promoter of the barley genes 1-FFT, 1-SST, 6-SFT, and VI, binding motifs for TaMYB13 were identified, suggesting that the co-expression of these genes may be driven by a TaMYB13 homologue (Huynh et al., 2012). However, the HTH myb-type domain-containing protein (HORVU7Hr1G001830.3) identified here did not show similarity to TaMYB13 and we could not observe a clear homologue in barley. While three Myb-type TFs were also described to activate promoters of genes involved in fructan biosynthesis and degradation in chicory (Wei et al., 2017a, b), involvement of other TF family genes has not yet been reported.

Conclusions

A new genomic region and several causal SNPs involved in the regulation of barley grain fructan content were identified in this study. Further, our data hint towards different regulatory levels regarding the fructan type, the chain length, and the target tissue. The identified genomic region includes a physical cluster of functionally related fructan biosynthetic genes and several potential regulatory genes. While the clustering of fructan biosynthetic genes may hint at the co-evolution of these gene families, a conserved gene co-expression suggesting an equal contribution to grain fructan biosynthesis was not observed. Instead the spatiotemporal dynamics for fructan biosynthetic genes point towards versatile roles for the different fructan types. Phylogenetic relationships between fructosyltransferases and invertases within Poaceae suggest that 6-SFT may have evolved from a Poaceae ancestor genome after the major clade of vacuolar invertases diverged, followed then by 1-FFT and 1-SST (Huynh et al., 2012). The analysis also showed the presence of a unique barley clade of four vacuolar invertase genes, among them the newly annotated 6G-FFT, between the 6-SFT and the 1-FFT and 1-SST clades, suggesting that extra duplication might have occurred in barley. Accordingly, in developing grain, we observed similar co-expression with a set of TFs for 1-FFT, 1-SST, and 6-SFT, which was different from the associations found for 6G-FFT. The proposed dynamics of fructan biosynthesis in barley grain and potential regulators are presented in Fig. 6. Assuming a specific spatiotemporal control of grain fructan biosynthesis, breeding or genetic engineering for high fructan content related to grain-specific traits (e.g. nutritional quality or germination) will require careful approaches targeting certain tissues and developmental stages, as recently suggested for engineering mixed linkage (1,3;1,4)-β-glucan biosynthesis in the endosperm (Lim et al., 2019).

Supplementary data

The following supplementary data are available at JXB online.

Fig. S1. Representative chromatogram of mature barley grain soluble carbohydrates.

Fig. S2. Metabolite distribution among the lines.

Fig. S3. Oligosaccharide profile prototypes as obtained by neural gas clustering.

Fig. S4. Manhattan plots and box plots for all significant results from the ratio-GWAS.

Fig. S5. Box plots for significant effects of SNPs identified from exome capture data.

Table S1. List of germplasms with their corresponding oligosaccharide profile group.

Table S2. Oligosaccharide annotation information.

Table S3. Oligosaccharide peak area entry mean values.

Table S4. Correlation and significance values for associations among the 27 metabolites.

Table S5. List of all detected gene models in the significant QTL interval.

Table S6. Results from exome capture data evaluation.

Table S7. Heatmap of correlations for the expression of fructan genes with other gene models from the QTL interval for all developmental phases and tissues.

Table S8. Expression data and Pearson correlation values for associations among the gene models from the GWA scan for all developmental phases and tissues.

Table S9. Absolute quantification of soluble sugars and fructans

Abbreviations

Abbreviations
 
• DAP

  days after pollination

 
• DM

  dry matter

 
• DP

  degree of polymerization

 
• ELSD

  evaporative light scattering detection

 
• FDR

  false discovery rate

 
• 1-FFT

  fructan:fructan 1-fructosyltransferase

 
• FOS

  fructooligosaccharides

 
• FPKM

  fragments per kilobase, per million mapped reads

 
• 6G-FFT

  6(G)-fructosyltransferase

 
• GWA

  genome-wide association

 
• GWAS

  genome-wide association study

 
• hai

  hours after imbibition

 
• HPAEC-PAD

  high pH anion exchange chromatography with pulsed amperometric detection

 
• KP

  kestopentaose

 
• KT

  kestotetraose

 
• LD

  linkage disequilibrium

 
• LOD

  logarithm of odds

 
• MAF

  minimum allele frequency

 
• NG

  neural gas

 
• PEG

  polyethylene glycol

 
• QTL

  quantitative trait locus

 
• RFO

  raffinose family oligosaccharide

 
• 6-SFT

  sucrose:fructan 6-fructosyltransferase

 
• SNP

  single nucleotide polymorphisms

 
• SPE

  solid phase extraction

 
• 1-SST

  sucrose:sucrose 1-fructosyltransferase

 
• TFA

  trifluoroacetic acid

 
• TPM

  transcripts per million

 
• VI-1

  vacuolar invertase1

Acknowledgements

This work was supported by grants from the German Research Foundation (DFG, MA 4814/3-1 and 3-2) to AM and from the Australian Research Council (ARC) Centre of Excellence in Plant Energy Biology (http://www.plantenergy.uwa.edu.au/) to RAB. KH and RW acknowledge the support of the Rural & Environment Science & Analytical Services Division of the Scottish Government We greatly acknowledge the excellent technical assistance of Shi F. Khor and Jelle Lahnstein. We wish to thank Natalie S. Betts and Helen M. Collins for their help with the geminated grain dataset. The authors acknowledge the use of the facilities, and scientific and technical assistance of the Australian Plant Phenomics Facility, which is supported by the Australian Government’s National Collaborative Research Infrastructure Strategy (NCRIS).

Author contributions

AM and RAB designed and developed the concept of the study. RW provided the material and genomic data of the barley panel and was involved in the design of the GWAS. AM and BB conducted the growth experiments and harvested the mature grain material at The Plant Accelerator. AM performed the oligosaccharide profiling analysis and evaluated the HPAEC-PAD data. AM and KW conducted the experiments related to the identification of fructan structures (isolation and MS identification). US performed the clustering of the data and the distribution analysis. KH performed the GWAS and ratio-GWAS, as well as the evaluation of the exome capture data. MRT, MKA, and LGW conducted the transcriptomic analyses of pistil tissues and developing grain tissues. MS and RW conducted the transcriptomic analysis of developing anther tissues. AM evaluated the transcript expression datasets provided, and conducted the co-expression analyses. AM, KH, and RAB wrote and provided a draft manuscript which has been revised and accepted by all authors.

Conflict of interest

All authors state no conflict of interest concerning this paper.

Data availability

All data supporting the findings of this study are available within the paper and within its supplementary data published online.