Using transcriptome sequencing (RNA-Seq) to screen genes involved in β-glucan biosynthesis and accumulation during oat seed development

Plant materials and sampling of oat seeds

The plant material used in this study was the high β-glucan content variety BY and the low β-glucan content variety DY. These varieties were selected from among 130 oat germplasm resources (Table S1). Seeds were planted in the experimental field of Inner Mongolia Agricultural University, Hohhot, Inner Mongolia (40°80′77″N, 111°62′29″E). The two varieties were sowed on April 20 and April 30, 2021, respectively, to achieve a consistent flowering period. A randomized zonal design with three replications was used. The area of the plot is 5 m*5 m = 25 m² with a row spacing of 30 cm. A total of 50 plants of each variety with the same growth were randomly selected at the spike stage for numbering and listing. To ensure the consistency of the seed growth period, the other whorl layers were cut off with scissors and only the first whorl layer was retained. We treated day 7 after flowering as day 1 of the seed-filling period (1 d). Seed samples were collected at five seed-filling periods: 1 d (BY1 and DY1), 6 d (BY2 and DY2), 11 d (BY3 and DY3), 16 d (BY4 and DY4), and 21 d (BY5 and DY5), and harvest mature seeds (26 d). Seeds collected on 1 d were used as the control. The samples were placed into a sterile frozen storage tube, quickly frozen with liquid nitrogen, and returned to the laboratory for storage in an ultra-low temperature refrigerator (−80 °C). Three biological replicates were taken for each sample, and a total of 30 samples were obtained.

Determination of β-glucan content in oat seeds

The β-glucan content of oat seeds was determined using the Megazyme^® mixed-linkage β-glucan assay kit (Biostest, Beijing, China), and the experiment was performed according to the kit instructions with three replicates.

RNA extraction and cDNA library construction

The total RNA of the oat seeds was extracted using RNA Prep Pure Plant Kit (Tiangen, Tianjin, China), using approximately 0.1 g of seeds. RNA integrity was determined by 1% agarose gel electrophoresis, and RNA integrity values (>6.5) were detected by Agilent 5400 Bioanalyzer. RNA concentration and purity were measured using a NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA). Oligo(dT) magnetic beads were used to enrich the mRNA with poly-A tails, and mRNA was chemically cleaved using divalent metal ion solutions. Randomly fragmented mRNA was used as a template to synthesize the first and second cDNA strands. cDNA fragments were enriched with 14 PCR cycles and library concentrations were precisely quantified to ensure that the libraries were available for subsequent sequencing. A total of 30 cDNA libraries were constructed.

RNA-seq data processing

High-throughput double-end sequencing was performed using Illumina NovaSeq 6000 (Illumina, USA). Reads with connectors, undeterminable base information, low quality (Q<20), and lengths less than 20 bp were filtered, and the data were assembled and spliced for subsequent bioinformatics analysis. The obtained clean data were compared to the reference genome (https://wheat.pw.usda.gov/GG3/graingenes-downloads/pepsico-oat-ot3098-v2-files-2021) using HISAT2 (Mortazavi et al., 2008). Fragments per kilobase of transcript per million mapped reads (FPKM) were calculated to evaluate the gene expression level of each sample. Differentially expressed gene (DEGs) analysis was performed between comparative combinations using DESeq2 (1.20.0) software. P-adj ≤0.05 and —log₂FC—≥1 were set as the thresholds of significant differential expression (Mao et al., 2005). Gene Ontology (GO) enrichment analysis, Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis, and gene expression analysis were performed using the R package embedded in Novomagic (https://magic.novogene.com/). The data were plotted graphically.

Weighted gene co-expression network analysis

To identify DEGs associated with β-glucan more comprehensively, weighted gene co-expression network analysis (WGCNA) was applied to RNA-seq data. The expression levels of differentially expressed genes in different samples were clustered by WGCNA, and correlation analysis was performed to find the set of genes directly associated with β-glucan content in oat seeds.

Detecting the gene expression level by quantitative real-time PCR (RT-qPCR) analysis

Total RNA was extracted from oat seeds using RNA prep Pure Plant Kit (Tiangen, Tianjin, China) using approximately 0.1 g of seeds. RNA integrity was determined by 1% agarose gel electrophoresis, and RNA concentration and purity were determined by NanoDrop2000 spectrophotometer. RNA with an A260/280 ratio of 1.8–2.0 and an A260/A230 ratio of 2–2.5 may be used for subsequent experiments (Table S2). The actin (KP257585.1) gene of oats was used as the internal reference gene, and primers were designed for the reference genes and differentially expressed genes according to the RT-qPCR primer design principle using Primer 3 Plus (https://www.primer3plus.com/), which were synthesized by Sangon Biotech (Shanghai, China), and the specificity was verified in pre-experiments (Table S3). The first strand of cDNA was synthesized with TIANGEN^® FasKting cDNA (Tiangen, Tianjin, China) First Strand Synthesis Kit (RNA 2 µl; 5 ×gDNA Buffer 2 µl; RNase-Free ddH₂O 6 µl; 10 ×King RT Buffer 2 µl; FastKing RT Enzyme Mix 1 µl; FQ-RT Primer Mix 2 µl; RNase-Free ddH₂O 5 µl.), incubated at 42 °C for 15 min, then 95 °C for 3 min. The cDNA obtained was subjected to fluorescence quantitative PCR using the SGExcel FastSYBR qPCR Mix (Sangon Biotech, Shanghai, China) KitI (1 ×MonAmpTM SYBR^® Green qPCR Mix 10 µl; 10 µM forward primer 0.4 µl; 10 µM reverse primer 0.4 µl; cDNA 2 µl; RNase-Free Water 7.2 µl). The reaction procedure was: 95 °C pre-denaturation for 30 s, denaturation at 95 °C for 30 s, annealing at 60 °C for 40 s, extension at 72 °C for 30 s. The last three steps were performed for a total of 40 cycles, three replicates for each sample. The lysis curves and amplification curves were collected to verify primer specificity. Finally, the relative expression values were calculated by the 2^−ΔΔCq method. To ensure data availability, all 30 samples were subjected to four concurrent identical technical replications.

Statistical analysis

We analyzed the data using one-way ANOVA for significance using IBM SPSS Statistics 25 (Meyers, Gamst & Guarino, 2013). Origin 2021 (https://www.originlab.com/) was used by us to plot dot lines and bar plots. We used Novomagic (https://magic.novogene.com/) to analyze the transcriptome data and plot all charts.

Dynamics of β-glucan content in oat seeds at different filling stages

In order to study the accumulation pattern of β-glucan in oat seeds during the filling period (1 d–26 d), we determined the β-glucan content in BY and DY seeds at different filling periods (Fig. 1, Table S4). The results showed that the β-glucan content of the two varieties maintained a gradually increasing trend during the filling period. The accumulation rate of β-glucan in both varieties showed the characteristic of “fast-slow-fast”. The period from days 6–16 was the stage with the fastest accumulation rate of β-glucan, and more than 72% and 70% of the β-glucan in BY and DY, respectively, were accumulated in this stage. From 6 d onwards, both varieties showed significant differences (P ≤ 0.001) in β-glucan content. The β-glucan content of BY was 5.414% in mature seeds (26 d), which was significantly higher than that of DY (2.092%).

Transcriptome sequencing and data quality control

Illumina sequencing data quality analysis is shown in Table 1. There were 2,216,337,008 original downstream sequences obtained from 30 samples, with a total valid sequence data volume of 324.09 G and data volume variability between 10.02 G–11.66 G; the data error rate of all 30 samples was about 0.03%; the Q30 content was between 92.73%–93.73; the GC content ranged from 53.14%–57.66%. The sequencing was of high quality and fully satisfied the requirements of subsequent analysis.

To assess the correlation between the three biological replicates, we performed Pearson linear correlation analysis on the sequencing data of three biological replicates of each variety and plotted the correlation heat map (Fig. 2A). The results indicated that the gene expression similarity among the three biological replicates of the 10 samples was high enough for subsequent analysis of differential genes. Principal component analysis was then performed (Fig. 2B). The tight clustering between the three biological replicates of each sample indicated that the samples were reproducible, and the two varieties were separated from each other, indicating that there were significant differences in transcriptional programs between the two varieties at each filling stage.

Screening and functional annotation of DEGs among different species

A two-by-two comparison of BY and DY was performed for each filling period to investigate the transcriptional regulatory mechanisms associated with the synthesis and accumulation of β-glucan in both genotypes. A total of 8,443, 8,866, 7,137, 8,875, and 5,522 DEGs were obtained for the five comparison groups (BY1 vs. DY1, BY2 vs. DY2, BY3 vs. DY3, BY4 vs. DY4 and BY5 vs. DY5, respectively), of which 2,066 were common to the five periods (Fig. 3A, Table S5). Further analysis of the DEGs revealed that the number of up-regulated genes was greater than the number of down-regulated genes on 1 d, 6 d, and 21 d of grouting; this trend was reversed on 16 d, and the number of up-regulated and down-regulated genes was almost identical on 11 d (Fig. 3B). A cross-sectional comparison between the five comparison groups revealed that 1,459 DEGs were consistently up-regulated as the grouting process proceeded and 605 genes were consistently down-regulated with the perfusion process (Table S6).

All of the screened DEGs were enriched for GO functions, and a total of 213 GO-terms were significantly (P-value < 0.05) enriched for the five comparison groups (Fig. S1). These terms involved the glucan metabolic process, cellular polysaccharide metabolic process, cellular glucan metabolic process, cell wall, protein heterodimerization activity, xyloglucan:xyloglucosyl transferase activity, and nutrient reservoir activity GO.

KEGG pathway enrichment was performed for the DEGs screened in each comparison group and a total of 21 significantly (P-value < 0.05) enriched pathways were obtained (Fig. 3C), of which protein processing in the endoplasmic reticulum, pentose phosphate pathway and flavonoid biosynthesis were significantly enriched in at least two perfusion periods.

Screening and functional annotation of DEGs during seed development

To investigate the DEGs associated with the synthesis and accumulation of β-glucan in oat seeds at different filling stages, we used the gene expression levels of DY and BY at the 1 d of filling as controls. We compared the gene expression changes of the two varieties at 6, 11, 16, and 21 days of filling. A total of 77,534 DEGs were obtained from the eight comparison groups (BY2 vs. BY1, BY3 vs. BY1, BY4 vs. BY1, BY5 vs. BY1, DY2 vs. DY1, DY3 vs. DY1, DY4 vs. DY1, DY5 vs. DY1), of which 38,675 were obtained for BY and 38,859 for DY. A total of 1,557 and 1,656 DEGs were co-expressed for BY and DY at the four perfusion stages, respectively (Fig. 4A, Table S7). Further analysis of the expression of all DEGs revealed that the number of DEGs in both varieties tended to increase as the filling process proceeded, and the number of up-regulated genes was greater than that of down-regulated genes in each period (Fig. 4B, Table S8).

All screened DEGs were enriched for GO function, and DEGs of BY and DY were significantly (P-value < 0.05) enriched for 203 and 284 GO-terms, respectively (Fig. S2). The main GO entries included responses to abiotic stimulus, cell wall, and protein heterodimerization activity. Notably, DEGs of DY at 11 d, 16 d, and 21 d were significantly enriched for the GO term xyloglucan:xyloglucosyl transferase activity, whereas BY was significantly enriched in this entry only at 6 d.

KEGG pathway enrichment was performed for DEGs screened in each comparison group, and DEGs from BY and DY were significantly (P-value < 0.05) enriched into 16 and 18 pathways, respectively (Fig. 4C). Of these 34 significantly enriched pathways, 11 are common to both species, including starch and sucrose metabolism, pentose and glucuronate interconversions, protein processing in the endoplasmic reticulum, and galactose metabolism, protein processing in the endoplasmic reticulum, and galactose metabolism.

Candidate genes related to β-glucan biosynthesis and accumulation during seed development

Starch and sucrose metabolism (map00500) has been shown to be closely associated with β-glucan biosynthesis and accumulation, but β-glucan biosynthesis and accumulation is a complex biological process. We analyzed the expression pattern of each gene in the starch and sucrose metabolism pathway (Fig. 5A) to further explore the key genes for β-glucan biosynthesis and accumulation in oat seeds. The results showed that the expression of TRINITY_DN34723, TRINITY_DN15420, and TRINITY_DN34252, which encode glucan endonucleases, decreased in both genotypes of oats as the seeds developed, and the expression in BY was lower than that in DY, which was consistent with the change in β-glucan content of the seeds.

It is well documented that members of two cellulose synthase-like (CslF and CslH) subfamilies are responsible for β-glucan synthesis (Liepman & Cavalier, 2012; Doblin et al., 2009; Schreiber et al., 2014). Therefore, we screened the DEGs of the five filling periods of BY and DY and identified a total of 56 CslF and seven CslH genes. However, only one CslF2, two CslF8, and one CslF9 gene reached a certain expression abundance (Fig. 5B). Among them, the gene TRINITY_DN16295, encoding CslF2, was expressed on average 4.6 times more than DY in all periods of BY, and the expression of this gene showed a continuous increase with the development of the seed filling process.

Transcription factors related to β-glucan biosynthesis and accumulation during seed development

Different transcription factors (TFs) regulate a range of life processes during the growth and development of a plant. To identify the TFs involved in regulating β-glucan biosynthesis and accumulation during oat seed development, we investigated the genes encoding TFs in the RNA-Seq dataset. The low β-glucan genotype oat DY was used as a control to compare the FPKM values of TFs at each developmental period. The results showed that there were a total of 2,256 up-regulated TFs expressed in the five comparison groups, of which 213 TFs were common to the five comparison groups (Fig. S3). These TFs belong to 41 families such as bHLH, WRKY, NAC, MYB-related, bZIP, and M-type-MADS, respectively (Table S9). We found that these 13 TFs belonged to 10 families including B3, FAR1 and NAC. Some of them were up-regulated expression in BY, while they were not expressed in DY. The expression pattern of another part in DY was consistent with that of BY, but the expression amount was much lower than that of BY (Fig. 6A). In addition to this, we found 19 TFs belonging to 13 families respectively whose expression first increased and then decreased (Fig. 6B). Their expression peaked at 11–16 d, which was consistent with the trend in the accumulation rate of β-glucan. Notably, five of these 19 TFs belonged to the Dof family in common.

Weighted gene co-expression network analysis of DEGs

A total of 39,211 DEGs were clustered based on gene co-expression patterns. To make the constructed network more consistent with the scale-free topology, we chose five as the optimal soft threshold (Fig. S4), identified a network with 32 co-expression modules (Fig. 7A), and analyzed the characteristic genes of each co-expression module and their relationship with β-glucan content in oat seeds (Fig. 7B). Five modules associated with β-glucan were clearly identified (p < 0.05), namely the blue module (r = 0.93, p = 9 ×10⁻⁵), magenta module (r = 0.82, p = 0.004), light cyan module (r = 0.74, p = 0.01), white module (r = 0.70, p = 0.03) and dark green module (r = 0.67, p = 0.03). These modules were significantly and positively correlated with β-glucan content in oat seeds and contained 6,559, 1,223, 364, 80, and 139 DEGs, respectively.

RT-qPCR validation of RNA-Seq data

We validated RT-qPCR on seven randomly selected genes using actin as an internal reference and performed a linear fit of RT-qPCR to RNA-Seq data, which showed that R² was 0.81 and 0.48 in BY and DY, respectively. The Pearson correlation coefficients were 0.90 and 0.69, respectively. This indicates that the RNA-Seq results were reliable (Fig. 8, Table S10).