1. Introduction

Starch is the most abundant photosynthetic product of green plants, serving as the main calorie source for humans, livestock, and other animals (Ren et al., 2021). Additionally, starch, including components of amylose and amylopectin, is also diversely used for various purposes in both food and non-food industries (Bahaji et al., 2014; Zeeman et al., 2010). Amylose is essentially a linear-chain linked by α-1, 4-glycosidic bonds composed of thousands of glucose residues, while amylopectin is a larger polymer composed of α-1, 4-linked glucose chains with α-1, 6-glycosidic bonds (Buléon et al., 1998; Li and Gilbert, 2018; Zeeman et al., 2010). In green plants, starch is generally biosynthesized through photosynthesis that integrates a series of functional enzymes, mainly including ADP-glucose pyrophosphorylase (AGPase), Brittle-1 (BT1), starch synthase (SS), starch branching enzyme (SBE), starch-debranching enzyme (DBE), and starch phosphorylase (SP) (Huang et al., 2021; Jeon et al., 2010; Zeeman et al., 2010). Starch accumulation, as well as the content ratio of amylose to amylopectin in storage organs directly affect the final performances of yield and quality of cereals (Ma et al., 2023).

Among the functional enzymes involved in the starch biosynthesis pathway, AGPase catalyzes ADPG synthesis, providing a direct substrate for extension of the sugar chain (Bahaji et al., 2014). BT1 serves as a transmembrane protein that transports ADPG into the amyloplasts for the extension of sugar chains (Shannon et al., 1998). SSs are responsible for the extension of liner chains and the extension of branched chains of different lengths (Delrue et al., 1992; Jeon et al., 2010; Keeling and Myers, 2010). SBEs catalyze the formation of α-1, 6-glycosidic bonds to form polysaccharide chain branches, thereby forming amylopectin (Blauth et al., 2001; Keeling and Myers, 2010). DBEs mainly participate in amylopectin synthesis via hydrolyzing the wrong branch chains, involving starch granule formation and starch degradation (Burton et al., 2002; Lü et al., 2008; Zeeman et al., 2007). SPs function by forming complexes with pathway related enzymes to affect starch biosynthesis (Nakamura et al., 2012; Tickle et al., 2009).

Starch biosynthesis in green plants can orderly proceed under the combination of all functional enzymes, and this process can also be affected directly by mutations of related coding genes (Huang et al., 2021; Jeon et al., 2010; Keeling and Myers, 2010). Meanwhile, the regulation of starch biosynthesis related enzymes or coding genes is equally important, including methylation, allosteric, phosphorylation, phytohormones, transcriptional regulation and environmental factors (Huang et al., 2021; Kötting et al., 2010). Low expression of starch biosynthesis related genes (SBRGs) is highly methylated in the rice and maize endosperm (Hu et al., 2021; Zemach et al., 2010). Allosteric and phosphorylation regulation mainly functions on the proteins, including enzyme activities, protein stability and interactions between or among proteins (Ahmed et al., 2015; Grimaud et al., 2008; Tetlow et al., 2004). For example, 3-PGA activates AGPase while Pi inhibits it (Boehlein et al., 2009), and trehalose-6-phosphate (T6P) regulates AGPase redox activity (Kolbe et al., 2005). Meanwhile, in wheat and maize, phosphorylation of key enzymes involved in starch biosynthesis can regulate enzyme activity and complex formation, including TaSSII, TaSBEII, ZmGBSSI, ZmSBEIIb, ZmPHO1, large and small subunits of AGPase in maize, and ZmBT1 (Grimaud et al., 2008; Makhmoudova et al., 2014; Tetlow et al., 2004, 2008; Walley et al., 2013).

Sucrose and ABA also play vital roles in starch biosynthesis in the endosperm (Huang et al., 2021). Sucrose can directly induce transcription of key genes involved in starch biosynthesis (Baguma et al., 2008; Chen et al., 2011). Cis-elements induced by ABA are also found in promoter regions of SBRGs, inducing their expressions (Chen et al., 2011; Hu et al., 2012). The synergistic treatment of ABA and sucrose can also induce transcription of key genes involved in maize starch biosynthesis (Huang et al., 2016). Transcription factors (TFs) that regulate transcription of SBRGs have been identified in rice, maize, wheat, and other crops, and are mainly distributed in the TF families, such as NAC, AP2, DOF, WRKY, bZIP, MADS, bHLH, NF–Y, and MYB (Huang et al., 2021). While TFs mediating starch biosynthesis regulationare still limited in sorghum grains.

Sorghum (Sorghum bicolor L.) is a typical C4 crop that belongs to the Poaceae, and is the 5th most widely cultivated cereal crop around the world. Sorghum is also an important staple food crop in many semiarid areas of the world, and a major feed crop for livestock (Kang et al., 2023; Rakshit et al., 2014). Sorghum grain is comprised of multiple compounds, among which starch and protein are the majors, accounting for 70–75% and 8–15%, respectively (Kang et al., 2023; Pontieri et al., 2022; Zhu, 2014). Therefore, starch accumulation in sorghum grains determines the final performances of sorghum yield and quality. Similar to starch of other crops, sorghum starch is mainly composed of amylose and amylopectin, and also aggregates functional enzymes such as AGPase, SSs, SBEs, and DBEs for biosynthesis (Kang et al., 2023; Xiao et al., 2022a). However, research on the regulatory mechanism of starch biosynthesis in sorghum grains lags far behind rice, maize, wheat, and other crops.

We measured the contents of starch and sucrose in sorghum grains at different developmental stages after pollination, and conducted omics-based analysis of RNA-seq and Assay for Transposase Accessible Chromatin with high-throughput sequencing (ATAC-seq) to dissect starch biosynthesis regulation in sorghum grains. Furthermore, we identified a TF of SbNAC68 that belongs to the NAC family served as the vital protein for transcriptional regulation for starch biosynthesis in sorghum grains. The summary results of the present study provide informative references for the in-depth discovering of the regulatory mechanisms for starch biosynthesis in sorghum grains.

2. Materials and methods 2.1. Plant materials and growth conditions

Sorghum variety BTx623 was grown in the lab field of college farm, Southwest University (Chongqing, China). We labelled the plants and collected sorghum grain samples at the initiation of pollination (0 DAP), and then 3, 6, 9, 12, 15, 20, 25, and 30 days after pollination (DAP). All collected samples were immediately frozen in liquid nitrogen and then stored at −80 ℃. Three biological replicates were prepared each sample.

2.2. Determination of starch and sucrose contents

Total starch content in sorghum grain samples was determined by acidolysis-DNS. Amylose and amylopectin contents in sorghum grain samples were determined through the dual wavelength colorimetric method. Anthrone-sulfuric acid colorimetry was used to determine sorghum grain sucrose contents. SpectraMax 190 (Molecular Devices, USA) was used for absorbance detection. The formula [(C × V_t × D)/(V × m)] × 0.9 was used for calculating total starch content (mg/g). In this formula, C (mg) was the reducing sugar content of the sample in the colorimetric tube obtained according to the standard curve, V_t (mL) was the total volume of sample extraction, V (mL) was the amount of liquid taken, D was the dilution ratio, m (g) was the sample weight, while 0.9 was the conversion factor of glucose to starch. The amylose/amylopectin content (mg/g) was calculated by the formula of C × (V_t/V) × D/m, where C (mg) referred to the starch content corresponding to the sample in the colorimetric tube according to the standard curve. The sucrose content (mg/g) was calculated by the formula of [(C × V_t × D)/(V × m)] × 0.001, where C (mg) referred to sucrose content corresponding to the sample in the colorimetric tube according to the standard curve, while 0.001 was the conversion ratio of μg to mg.

2.3. Determination of enzyme activities

The activities of ADP-glucose pyrophosphorylase (AGP), granule-bound starch synthase (GBSS), and soluble starch synthase (SSS) were determined according to the manufacturer's protocol provided by Suzhou Comin Biotechnology Co., Ltd. (Suzhou, China). 0.1 g of fresh sorghum grains of 5, 10, 15, 20, 25, and 30 DAP were ground evenly and mixed with 1 mL of extraction solution at low temperature. The mixtures were centrifuged at 10, 000 g and 4 ℃ for 10 min, and supernatants were placed on ice for testing. Infinite M200Pro (TECAN, USA) was used for absorbance detection with the wavelength set to 340 nm. Enzyme activity was calculated according to the formula in the instructions.

2.4. RNA isolation and library preparation

Total RNA was extracted through TRIzol reagent (Invitrogen, CA, USA) according to the manufacturer's protocol. A NanoDrop 2000 spectrophotometer (Thermo Scientific, USA) and an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA) were used to detect the quantification and integrity of extracted RNA. Libraries were constructed using a VAHTS Universal v.6 RNA-seq Library Prep Kit according to the manufacturer's instructions.

2.5. RNA-seq

RNA-seq and related analysis were conducted by OE Biotech Co., Ltd. (Shanghai, China). Libraries were sequenced on an Ilumina Novaseq 6000 platform and 150 bp paired-end reads were generated. About 50 M raw reads for each sample were generated. SOAPnuke (v.2.1.0) was used to filter the raw reads to obtain clean reads (Chen et al., 2018). Then about 47 M clean reads for each sample were retained for subsequent analyses. Clean reads were mapped to the Sorghum_bicolor_NCBIv3 using HISAT2 (v.2.2.2.1, Kim et al., 2015). Fragments per kilobase of transcript per million mapped reads (FPKM) of each gene was calculated and each gene read count was obtained by HTSeq-count (Anders et al., 2015; Roberts A et al., 2011). Principal component analysis (PCA) was performed using R (v.3.2.0) to evaluate the biological duplication of samples. To understand the differences among gene expression modules in sorghum grain, we performed Short Time series Expression Miner (STEM) analysis for the transcriptome data (Ernst and Bar-Joseph, 2006). All the genes not expressed at any of nine development stages from 0 to 30 DAP were excluded from the STEM analysis.

2.6. Differentially expressed gene analysis

DESeq2 (v.1.22.2) were used to perform differential expression analysis (Love et al., 2014). Q value < 0.05 and foldchange > 2 or foldchange < 0.5 were set as the threshold for differentially expressed genes (DEGs) screening. Foldchange (FC) here was calculated as the FPKM ratio of grain sample pairs at the corresponding development stages, i.e., FPKM-0DAP/FPKM-3DAP. At a certain development stage, for example 0DAP, if the FCs of some gene(s), the FPKM-0DAP versus FPKMs at other eight development stages, were all larger than 2 or all less than 0.5, then such gene(s) was/were determined as highly or low expressed gene(s), respectively. Hierarchical cluster analysis of DEGs was performed through R (v.3.2.0) to demonstrate the expression pattern of genes within different groups and samples. We used Gene Ontology (GO) (http://www.geneontology.org/) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway (http://www.genome.jp/kegg/) enrichment analysis of DEGs to screen the significantly enriched terms through R (v.3.2.0) (Consortium, 2019; Kanehisa et al., 2007). Column, chord, and bubble diagrams of the significantly enriched terms were also drawn by R (v.3.2.0). Gene Set Enrichment Analysis (GSEA) was performed through GSEA software (Subramanian et al., 2005).

2.7. ATAC-seq and enrichment analysis

Sorghum grains at 12 DAP were collected and ground into powder in liquid nitrogen, and incubated with lysis buffer at 4 ℃ for 10 min. The ATAC-seq libraries were constructed according to Bajic et al. (2018). To perform tagmentation, 50, 000 nuclides were added to the transposition reaction solution. Tn5 transposed DNA was purified by AMPure DNA magnetic beads and PCR amplification, and the final qualified library was sequenced on the Illumina Nova-seq platform (San Diego, CA, USA) with PE150 mode. Bowtie2 (v.2.3.5) was used to align Clean Data with the reference genome sequence, and Samtools (v.1.9) was used to summary the alignment results. When MACS3 was used to detect the candidate peaks, nomodel was selected, and the sift and extsize were set as 75 and 150, respectively. The annotations of detected peaks, visualizing of IGV, and detection of differential Peaks were correspondingly performed through the modules of ChIPseeker, Gviz, and Diffbind that implemented in R (v.3.2.0) with default settings. For example, the default setting for differential Peak detecting were Log2FoldChange > 0.58 or FDR < 0.05.

GO and KEGG were used for gene annotation and enrichment of the peak associated genes. MEME Suite was used for motif scanning and enrichment analysis of peak sequences in the JASPAR database (http://jaspar.genereg.net/). In addition, Venn analysis was performed between the gene set that selected from RNA-seq data profile with FPKM ≥ 1 (Xiao et al., 2022a) and that set with detected Peaks within 5′URT and promoter regions (≤ 3 Kb) from ATAC-seq data profile through the online Venn diagrams (https://bioinformatics.psb.ugent.be/webtools/Venn/).

2.8. qRT-PCR

Cultured protoplasts were used for total RNA extraction through TRIzol regent (Invitrogen, Carlsbad, CA, USA). First-strand cDNA was synthesized from 1.5 μg of total RNA using a PrimeScript^TM RT reagent Kit with gDNA Eraser (TaKaRa, Kusatsu, Japan) according to the instructions. The qRT-PCR assays were carried out according to Xiao et al. (2022a) with β-actin as the internal control, and performed via the Bio-Rad CFX96 real-time system. Relative transcription levels were calculated via the 2^−ΔΔCTmethod, possessed three biological replications, and p-values were calculated by t-test in R (v.3.2.0). Primers information for qRT-PCR assay was listed in Table S1.

2.9. Electrophoretic mobility shift assay (EMSA)

Based on the transcriptome analysis of different sorghum tissues (Xiao et al., 2022a), as well as the transcriptome data analysis of different sorghum grain developmental stages and co-expression analysis, SbNAC68 was identified as a candidate TF for transcriptional regulation of starch and sucrose related genes. The binding motifs of NAC family members were discovered from the ATAC-seq analysis results. Simultaneously the presence of ATAC-seq sites and NAC protein binding sites were identified in the promoter region (3000 bp upstream of ATG) of key genes involved in starch and sucrose metabolism. Promoter fragments with both ATAC-seq and NAC protein binding sites were further considered for EMSA analysis.

Prokaryotic expression vector pET32a was used for the prokaryotic expression of SbNAC68-His fusion protein, Escherichia coli strain of BL21 was selected as host cells for prokaryotic expression, and isopropyl β-D-1-thiogalactopyranoside (IPTG) was added to induce fusion protein expression at a final concentration of 0.5 mM. To purify the recombinant proteins, a Ni-Agarose His label Kit (CWBIO, Beijing, China) was used based on the manufacturer's instructions. A double-strand DNA fragment containing 5′-ACGCAA-3′ within the promoter regions of sorghum starch and sucrose metabolism related genes was used as the probe. All probes were synthesized with biotin labeling at 5′ end by Sangon Biotech Co. Ltd. (Shanghai, China). The binding reaction, competitive reaction, and EMSA detection were performed according to the instructions of a Chemiluminescent EMSA Kit (Beyotime Biotechnology, Shanghai, China). All sequence information of primers for vector construction was listed in Table S1.

2.10. Dual-luciferase assay in maize leaf protoplasts

The pGreenII0800-LUC double-reporter vector and pBI221 were used to detect the relationship between promoter activity and regulatory factors. The promoters of SbSUS4 (1920 bp), SbAGPS1 (1956 bp), SbSBEI (1869 bp), and SbSBEIIb (1971 bp) were subcloned into the pGreenII0800-LUC vector to drive Luciferase (Luc) expression. SbNAC68 was amplified and integrated into the plant expression vector pBI221, driven by the Ubi promoter. The pUbi-SbNAC68: pGreenII0800-Pro-Luc (1:2) was set as the experimental group, and pUbi-Gus: pGreenII0800-Pro-Luc (1:2) as the control group. All constructs were transformed into the maize protoplast prepared according to Xiao et al. (2021). LUC and Renilla (REN) luciferase activities were measured with a dual-luciferase assay kit (Yeasen Biotechnology, Shanghai, China) and analyzed with a GloMax_2020 (Thermo Fisher Scientific, Waltham, MA, United States). LUC/REN ratio was calculated to measure the relationship between experimental and control groups. Three independent experiments were performed, and each consisted of three replicates. The corresponding significant differences were detected by t-test of R (v.3.2.0). Primer information for vector construction was listed in Table S1.

2.11. Preparation of sorghum leaf protoplasts and transient overexpression

Enzymatic hydrolysis was used to prepare sorghum leaves protoplasts, and the enzymolysis liquid consisted of 0.5% Macerozyme R-10 (Yakult, Japan), 1.5% Cellulase R-10 (Yakult, Japan), 0.5 M mannitol, 10 mM CaCl₂, 10 mM 2-Morpholinoethanesulphonic acid monohydrate (MES monohydrate) and 0.1% Bovine serum albumin (BSA). The sorghum leaf samples that were cut into 1 mm wide pieces were immersed in the enzymolysis liquid and cracked at 25 ℃ at 50 rpm for 4 h. W5 buffer (125 mM CaCl₂, 154 mM NaCl, 5 mM KCl, and 2 mM MES) was used for sorghum leaf protoplast washing, and MMG buffer (0.4 M mannitol, 15 mM MgCl₂, and 4 mM CaCl₂) for sorghum leaf protoplast suspension. The transformation solution was composed of 40% PEG4000, 0.2 M mannitol, and 100 mM CaCl₂. The amount of transformed plasmid was 10 μg, and transformed protoplasts were cultured in the dark at 25 ℃ for 24 h.

3. Results 3.1. Determination of starch content and enzyme activity in sorghum grain

To investigate starch accumulation in sorghum grains, we characterized the dynamic changes in starch and sucrose content, and key enzyme activity at different developmental stages. KI-I₂ staining results indicated that starch gradually accumulated during sorghum grain development (Fig. 1A), which was also confirmed by starch content determination (Fig. 1B). Sucrose content initially increased between 5 (2.65 mg/g) and 10 DAP (13.65 mg/g), followed by a decrease that stabilized at around 8 mg/g (Fig. 1C). GBSS, SSS, and AGPase activities exhibited similar patterns during sorghum grain development, gradually increasing from 0 to 20 DAP (AGPase) and 0 to 25 DAP (GBSS and SSS), and then decreasing trends after 20 DAP (AGPase) and 25 DAP (GBSS and SSS) (Fig. 1D).

3.2. Dynamic RNA-seq of sorghum grains

To investigate transcriptome dynamics among different development stages of sorghum grains, we selected the whole grains at nine development stages of BTx623, i.e., 0, 3, 6, 9, 12, 15, 20, 25, and 30 DAP, for RNA-seq analysis (Fig. 2A).

A total of 188.59 G clean bases with Q30 > 95% were generated (Table S2), and 86.78% (Sb_20DAP_03) to 96.69% (Sb_3DAP_03) reads could be uniquely mapped to the Sorghum_bicolor_NCBIv3 reference genome (Table S2). Gene IDs with annotation to reference genome and levels of gene expressions were normalized as FPKM with the uniquely mapped reads (Table S3). Strategies of PCA, correlation analysis, and hierarchical clustering (HCL) were conducted to reveal the transcriptomic dynamics among different sorghum grain samples (Fig. S1). The PCA results implied significant differences across developmental stages (Fig. S1A). RNA-seq correlation analysis indicated obvious differences between pairs of grain samples. Correlation coefficients (CCs) among 3, 6, and 9 DAP grain sample pairs were > 0.7, those among 12, 15, 20, 25, and 30 DAP were > 0.8, while those among 3 to 9 DPAs and 12 to 30 DAP were < 0.3 (Fig. S1B). Results of hierarchical clustering based on the RNA-seq datasets were generally consistent with those from the correlation analysis. In summary, the overall development stages of sorghum grains could be divided into four periods (Fig. S1B), i.e., flowering (0 DAP), early grain development (3–6 DAP), critical filling (9–15 DAP), and latter grain development (20–30 DAP) (Fig. S1C).

3.3. Differential gene analysis of sorghum grains

Eight clusters were obtained based on the expression patterns from the results of RNA-seq among different sorghum tissues via the k-means clustering algorithm (Xiao et al., 2022a) (Fig. S2). Clusters 1 and 2 covered more genes, with transcripts detected in all tissues, while genes in clusters 4, 5, and 8 exhibited relatively high expression levels during sorghum grain development (Fig. S2).

Nine profiles with significant changes in gene expression trends were identified among 50 analyzed modules (Fig. 2B, Fig. S3; Table S4), covering 17, 362 genes and 1057 TFs (17362/1057). Among these profiles, Profile_10 covered the most genes and TFs (6837/417) that mainly enriched in photosynthesis, carbon fixation in photosynthetic organisms, and other pathways via KEGG (Fig. 2B, Fig. S3; Table S4). Notably, the KEGG results showed that genes and TFs in Profile_29 (1278/65) were mainly enriched in starch and sucrose metabolism and other pathways, and those of Profile_39 (2339/182) mainly enriched in MAPK and phytohormone signal pathways (Fig. 2B, Fig. S3; Table S4). Additionally, genes and TFs of Profile_39 were mainly expressed in grains after 9 DAP, some members of Profile_11, Profile_29, Profile_21, and Profile_40 mainly at middle to late grain development stages (20 and 30 DAP), while other genes and TFs of these profiles, and all genes and TFs of Profile_46/10/1/47 were mainly expressed at early development stages (0–15 DAP) in sorghum grains (Fig. 2B, Fig. S3; Table S4). Moreover, all identified 1057 TFs of nine profiles belonged to 51 families and exhibited differential expression patterns among different development stages of sorghum grains (Fig. 2C). Profile_10 contained the most TFs (417), among which 164 belonged to the bHLH family (Fig. 2C).

Meanwhile, a total of 32, 021 genes, including 1805 TFs were detected in sorghum grains (Tables S3 and S5), and 5789 with 334 TFs highly expressed in all 9 grain samples (Figs. 2D and S4A). The highly expressed genes including TFs among 0, 3, 6, 9, 12, 15, 20, 25, and 30 DAP grain samples were counted as 3926 (highly expressed genes)/221 (TFs), 339/24, 87/13, 16/1, 1/0, 1/0, 57/4, 30/1, and 987/70 (Fig. 2D). These results indicated that more genes and TFs were involved in the early (0 and 3 DAP) and late (30 DAP) development stages of sorghum grains (Fig. S4B). Co-expression analysis was conducted between TFs and DEGs and showed that within 3 and 6 DAP grain samples, 99% of TFs exhibited PCCs > 0.8 with at least 70% of highly expressed genes, while this proportion was significantly lower than those in grain samples at other development stages (Fig. S4C).

3.4. Transcription levels of SBRGs in sorghum grains

SBRGs expression patterns in sorghum grains were analyzed through the FPKM of RNA-seq (Fig. 3A; Table S6). There were six genes encoding the subunits of AGPase in sorghum, among which both SbAGPS1 and SbAGPS2 encoded the small subunits, while SbAGPLS1, SbAGPLS2, SbAGPLS3, and SbAGPLS4 encoded the large subunits. All these genes exhibited high expression levels during sorghum grain development (Fig. 3A and B; Table S6). Sorghum Brittle 1 (SbBT1), encoding ADPG transporter protein, was highly expressed from 9 to 30 DAP in sorghum grains (Fig. 3A and B, Table S6). There were 18 other genes that function on the extension and branching of sugar chains in sorghum grains (Xiao et al., 2022a). Among these genes, SbGBSSI, SbSSI, SSIIa, SbISA1, and SbSBEIIb were highly expressed in grains from 9 to 30 DAP; SbSSI, SbSBEI, and SbSBEIIa were high expressed during grain development; while the remaining 10 genes exhibited relatively lower expression levels (Fig. 3A and B; Table S6). Meanwhile, Nine genes of SbSPS1, SbSPS4, SbSPP1, SbSUS1, SbSUS4, SbCIN4, SbCIN5, SbCWIN5, and SbCWIN7 involved in the sucrose metabolism pathway possessed high expression levels in sorghum grains (Fig. 3A and B; Table S6). The differential transcriptional levels of SBRGs also indicated the existence of transcriptional regulation in starch biosynthesis in sorghum grains.

3.5. ATAC-seq analysis of 12DAP sorghum grains

Based on ATAC-seq of three biological replicates of 12-DAP sorghum grains through Illumina PE150, 59.07 Gb data was obtained, and a total of 84, 884 peaks were identified (Fig. 4A and Table S7). The distribution trend of ATAC-seq signal in the gene body region of three independent samples were summarized, and these peaks were associated with 28, 587 genes, among which 85.34% (24, 397) possessed peaks within 3 Kb regions of promoters and 5′UTRs (Fig. 4B, Table S8). Apart from 27.17% of detected peaks that fell into distal intergenic regions, the majority were located in gene related regions, i.e., 39.98% in 3 Kb regions adjacent to promoters, 12.33% in UTRs, 14.51% in both exons and intros, and 6.01% in downstream gene regions (Fig. 4C, Fig. S5). Though these peak-associated genes primarily enriched in the biosynthesis of acids and carbon metabolism, the pathways of plant hormone signal transduction, MAPK signaling, and starch and sucrose metabolism were still enriched via KEGG analysis (Fig. 4D). In addition, a total of 267 peak motifs were identified to associate with various TFs, among which 188 motifs were identified in all three replicates (Figs. 4E and S6; Table S9).

The conserved motifs identified via ATAC-seq analysis corresponded to binding sites of multiple TF families, including ERF, ARF, bZIP, DOF, MYB, TCP, and NAC. Notably, distinct motifs within the same TF family exhibited high sequence conservation. For example, ERF-type motifs were predominantly characterized by the core sequence of GCCG, WRKY-type by GTCAA, bHLH-type by CACGTG, DOF-type by AAGA, TCP-type by GG∗∗CC, and NAC-type by ACGCAA (Fig. 4F and Table S9). These findings indicate that TFs within the same family maintain conserved binding site recognition.

3.6. Combined analysis of ATAC-seq and RNA-seq revealed the binding sites of SBRGs

Based on RNA-seq of 12-DAP grains, a total of 19, 961 expressed genes (FMPK ≥ 1) were identified, most of which (16, 097, 80.64%) exhibited peaks within the promoter regions or 5′UTRs, similar high proportion of 65.98% (16, 097/24, 397) as those identified by ATAC-seq (Fig. 5A). GO and KEGG analysis results showed that starch biosynthesis was noted as an important pathway (Fig. 5B and C). KEGG enrichments supported that the top 6 enriched pathways of these 16, 097 genes fell into Ribosome, Plant hormone signal transduction, Spliceosome, Protein processing in endoplasmic reticulum, Endocytosis, and Starch and sucrose metabolism (Fig. 5C). Notably, 129 genes were enriched in starch and sucrose metabolism, which exhibited as one of the most enriched pathways among the remaining (Fig. S7). Several genes, including SbSUS1, SbAGPS1, SbAGPLS2, SbGBSSI, SbSSI, SbSSIIa, SbSBEI, SbSBEIIa, SbSBEIIb, SbISA1 and SbPHOL, were mainly highly expressed in the grain (Fig. 5D). These genes encode key enzymes involved in starch biosynthesis in sorghum grain.

Based on ATAC-seq, a total of 532 peaks associated 142 genes were involved in starch and sucrose metabolism, among which 182 association sites were located in the promoter regions, 152 peaks in intergenic, and 198 peaks were located in other regions including UTR, Extron, Intron, and Downstream (Fig. 5E and Table S10).

Meanwhile, 131 peaks were identified to associate with SBRGs (Fig. S7 and Table S11). Among these peaks, 44 were located in gene promoters (39) and 5′UTR (5), 61 in intergenic regions, and the remaining 26 peaks were predominantly within exon (13), intron (2), 3′UTR (7), and downstream 3000 bp sequences (4) of genes (Fig. S8 and Table S11). The 44 peak sites identified in the promoter regions and 5′UTR regions were directly distributed across 26 SBRGs, including SbAGPS1, SbAGPS2, SbAGPLS1, SbAGPLS2, SbBt1, SbGBSSI, SbSSI, SbSSIIa, SbSBEI, SbSBEIIb, SbISA1 (Fig. S9 and Table S11).

3.7. Key factors involved in transcriptional regulation of SBRGs

To identify candidate factors involved in transcriptional regulation of starch biosynthesis, we performed a coordinated analysis of multi-tissue RNA-seq data (Xiao et al., 2022a) and ATAC-Seq. The integrated approach identified 20, 129 genes with expression level of FPKM > 1 (Fig. 6A), and 575 grain-specific highly expressed genes (FPKM > 1) were characterized, including 61 TFs (Table S12). GO and KEGG enrichment confirmed significant representation of starch and sucrose metabolism pathways (Fig. 6C and D). The spatiotemporal expression patterns across tissues confirmed the enriched transcript of 61 TFs accumulation in developing grains (Fig. 6E). Co-expression analysis showed that multiple TFs, such as SbNAC68, SbNAC107, SbNAC38, co-expressed with starch biosynthesis and sucrose metabolism related genes (Fig. 6F).

Meanwhile, ATAC-seq profiling revealed 1760 peaks associated with these 575 genes. Among these peaks, 779 (44.26%) and 50 (2.84%) correspondingly distributed within the corresponding promoter regions (< 3 Kb upstream of TSS) and 5′UTR, 320 within overlapping structural sequences of genes (introns or/and exons), 485 in distal intergenic regions, and 126 within downstream of genes (Fig. S10A and Table S13). Notably, 61 grain-enriched TFs possessed 203 peaks, including 87 promoter-associated peaks, 48 intragenic peaks, and 68 peaks in distal intergenic/downstream 3000 bp regions (Fig. S10B and Table S14).

3.8. SbNAC68 directly participates in transcriptional regulation of SBRGs

Based on the multi-tissue RNA-seq analysis, several highly expressed TFs were identified in sorghum grain, including members of AP2, NAC, DOF, and MYB families (Xiao et al., 2022a). Combined with the RNA-seq data of sorghum grains at different developmental stages, eight NAC TFs exhibited much higher FPKM and specific expression patterns in the grains (Fig. S11), and four TFs were detected at highly transcript levels in sorghum grains after 9 DAP (Fig. S11B). The qRT-PCR further confirmed that SbNAC68 possessed abundant transcripts in the grains (Fig. 7A), as well as co-expressed with diverse SBRGs, including SbSUS2, SbSBEIIa, SbISA1, SbSPS1, SbSSI, SbAGPS2, SbSUS1, SbAGPLS1, SbGBSSI, and SbBt1 (Fig. 6F).

The ATAC-seq results revealed that NAC protein binding motif types existed the transposase-accessible chromatin, NAC028, NAC058, and NAC080 (Fig. 4F). Interestingly, these three NAC-type motifs were generally observed within 3000 bp upstream ATG of SbSPS1, SbSUS1, SbSUS2, SbAGPS2, SbAGPLS1, SbBt1, SbGBSSI, SbSSI, SbSBEIIa, and SbISA1 (Fig. 7B, Table S15). EMSA results showed that SbNAC68 directly bound to the conserved motif of 5'-ACGCAA-3' that identified through ATAC-seq (Fig. 7C), which was confirmed by the combination of SbNAC68 with promoters of SbGBSSI, SbSSIIa, SbAGPLS1, SbSBEI, SbAGPS1, and SbSS1 (Fig. 7C). Meanwhile, SbNAC68 promoted the activities of SbGBSSI and SbBt1 (Xiao et al., 2022a), could also promote the activities of SbSBEIIb and SbAGPS1, and repressed that of SbSBEI (Fig. 7D). Furthermore, transient overexpression of SbNAC68 in sorghum leaf protoplasts could significantly promote the transcription of some SBRGs, including SbSSV, SbAGPLS1, SbSSI, SbAGPLS4, SbGBSSIIa, and SbGBSSI, and significantly inhibit the transcription of some other SBRGs, i.e., SbSBEIIb, SbAGPS2, SbISA1, SbSSIIc, and SbSSIIIb (Fig. 7E). Meanwhile, transient overexpression of SbNAC68 in maize leaf protoplasts also affected the transcription levels of maize SBRGs, such as significantly increased the transcription levels of ZmGBSSI and ZmSSIIIa (Fig. S12), and promoted the transcription levels of ZmSSI and ZmSS2b. These results suggested the vital role of SbNAC68 in transcriptional regulation of starch biosynthesis in sorghum grains.

4. Discussion 4.1. Conservation of cereal starch biosynthesis

Starch is the vital carbohydrate for living organisms, serving as an energy storage carrier in plants and an energy supply source for animals. Cereal crops, such as maize, rice, and wheat, are the most utilized crop types, providing about 50% of human energy sources (Huang et al., 2021). For cereal crop utilization, grain related traits, including yield and quality, were the earliest concerns of humans. Starch is one of the main components of cereal crop grains, and its contents and composition are important for crop yield and quality (Huang et al., 2021; James et al., 2003; Jeon et al., 2010).

Starch is composed of two d-glucose (D-Glc) polymers in the form of amylose and amylopectin, and stored in tissues as insoluble semi-crystalline starch granules (Seung and Smith, 2019; Zeeman et al., 2010). In cereal crops, starch, including temporary and storage starch, are biosynthesized in photosynthetic and non-photosynthetic cell plastids (Jeon et al., 2010; Keeling and Myers, 2010; Zeeman et al., 2010). And grain development also exhibits similar trend among different cereals. For example, in both rice and wheat, 0 DAP represented the flowering, while grains will quickly initiate the development of embryo and endosperm in a few days after double fertilization (Yu et al., 2016; Liu et al., 2025). Starch begins to be biosynthesized and rapidly accumulate upon ~DAP15, while the grains will enter latter development stage and turn to maturation stage after DAP20 (Liu et al., 2025). Similar development trends were also observed via the omics-based profiling among different development stages of sorghum grains (Fig. S1C). The development similarities among cereal grains might imply the conversion of starch biosynthesis in the endosperm cells among the cereals.

Starch biosynthesis is highly conserved among cereal crops and depends on various enzymes and proteins, including AGPase, BT1, GBSS, SSS, SBE, DBE, and SP (Huang et al., 2021; Jeon et al., 2010). All these enzymes and functional proteins coding genes related to starch biosynthesis have been intensively documented in rice (Hirose and Terao, 2004; Ohdan et al., 2005), maize (Chen et al., 2011; Xiao et al., 2021), wheat (Stamova et al., 2009), barley (Radchuk et al., 2009), and sorghum (Kang et al., 2023; Xiao et al., 2022a). These functional enzymes possess conserved domains or sequence features that play crucial roles in substrate recognition or catalysis in driving functions. For example, the C-terminus of SSs contains conserved domains of Glyco-transfer_1 and Glyco-transfer_5, CBM_25 is the conserved domain of SSIII, and CBM_48 is a shared domain at the N-terminal of SBE and DBE (Jeon et al., 2010; Keeling and Myers, 2010; Xiao et al., 2022a). Meanwhile, changes in functional enzymes could affect starch biosynthesis, serving as an important potential way to improve crop yield (Bahaji et al., 2014; Ma et al., 2023). Certainly, starch biosynthesis improvement is also important to improve sorghum quality and yield (Huang et al., 2021; Kang et al., 2023). The expression patterns of genes related to sorghum starch biosynthesis (Fig. 3B and Fig. S13) are consistent with those of genes related to maize starch biosynthesis (Xiao et al., 2021). SbNAC68 can enhance the transcription of sorghum SBRGs (Fig. 7E), and transient overexpression in maize leaf protoplasts can also regulate the transcription of maize SBRGs (Fig. S12). Meanwhile, the same type of transcription factors, such as NAC TFs, have been found to regulate starch biosynthesis in the grains of crops like maize (Chen et al., 2023; Xiao et al., 2021), rice (Wang et al., 2020; Jin et al., 2023), and wheat (Liu et al., 2020; Wang et al., 2023). This also highlights the evolutionary conservation of starch biosynthesis and its regulatory.

Additionally, we also observed notable fluctuations of uniquely mapped percentages to sorghum reference genome among gene profiles of 0–30 DAP generated from RNA-seq. For example, even though high uniquely mapped percentages of grain samples at 20, 25, and 30 DAP (86.94%–89.38%) were recorded, these proportions were somewhat lower than those at 0, 3, 9 DAP, the early development stages of grains in sorghum (Table S2). This fluctuation might be caused by both development and maturation of grains. As grains develop and mature, more storage products, such as starch and protein, are accumulated, and the accumulation will in turn affect the purity of extracted RNA from the corresponding tissues. On the other hand, when the development status altered from early to later stages, gene activities in grains tend to reduce rapidly in wheat (Guan et al., 2022). Similar biological process might also occur in sorghum grains, which will also affect the total active RNAs genome-widely and cause the fluctuations in percentage between early and late development stages.

4.2. Transcriptional regulation plays important role in sorghum starch biosynthesis

Transcription levels of SBRGs are regulated by a series of TFs (Huang et al., 2021), such as TFs from NAC families, i.e., ZmNAC126/128/130 (Chen et al., 2023; Xiao et al., 2021; Zhang et al., 2019), OsNAC20/24/26 (Wang et al., 2020; Jin et al., 2023), TaNAC019/A19 (Liu et al., 2020; Wang et al., 2023); members of bZIP, i.e., OsbZIP58/RISBZ1/OsOpaque3 (Cao et al., 2022; Wang et al., 2013), ZmbZIP91/ZmOpaque2 (Chen et al., 2016; Zhang et al., 2016), TubZIP28 (Song et al., 2020), DOF family members of OsRPBF (Kawakatsu et al., 2009) and ZmDOF3/36 (Qi et al., 2016; Wu et al., 2019), AP2 family members of OsEBP-89/OsRSR1 (Fu and Xue, 2010; Zhu et al., 2003) and ZmEREB156 (Huang et al., 2016); and other TF family members of ZmMYB14 (Xiao et al., 2017), OsbHLH144/OsNF_YB1/OsNF_YC12/OsSGL (Bello et al., 2019; Liu et al., 2022), and TabHLH95/TaNF_YB1 (Liu et al., 2023). These TFs could positively regulate starch biosynthesis or negatively exert their effects, or they could simultaneously regulate the transcription of multiple key genes involved in starch biosynthesis or specifically target a single gene. However, the commonality of those TFs was that their gene transcripts could be detected in grains, or even expressed specifically in grains. We also found some TFs shared high expression levels in grains and co-expressed with SBRGs in sorghum grains (Fig. 6E and F). Furthermore, our previous work focused on SbDOF21 and SbNACs showed that these TFs could significantly increase the promoter activities of SBRGs in sorghum (Xiao et al., 2022a, 2022b). We further verified that SbNAC68 could activate SBRGs promoters in sorghum grains, and regulate the transcription of some genes in sorghum leaf protoplasts (Fig. 7E). Additionally, transient overexpression of SbNAC68 in maize leaf protoplasts could also significantly promote the transcription of ZmGBSSI (Fig. S12).

TFs are proteins involved in gene transcriptional regulation, and cis-elements are crucial sites for TF binding. In rice, maize, and wheat, TFs could bind to specific motifs to regulate gene expressions. For example, the motifs of 5′-ACGCAA-3′, 5′-CACG-3′, and 5′-TTGACAA-3′ were identified as the binding sites of NAC TFs (Jin et al., 2023; Wang et al., 2020; Xiao et al., 2021; Zhang et al., 2019). 5′-ACTCAT-3′ served as the binding site of ZmbZIP91 (Chen et al., 2016), 5′-TGACG-3′ was the binding site of OPAQUE3 (Cao et al., 2022), and 5′-TGTAAAG-3′ (P-Box) was the binding site of DOF proteins (Qi et al., 2016; Xiao et al., 2022b). In the present study, vast sites within the promoter regions of SBRGs in sorghum grains were identified as specific binding sites to different TF families, such as DOF, ERF, MYB, and NAC (Table S18). Additionally, three NAC binding types of NAC028, NAC058, and NAC080 identified through ATAC-seq could be detected among numerous gene promoters involved in sorghum grain starch biosynthesis (Fig. 7B and Table S15), with conserved motif sequence of 5′-ACGCAA-3′ (Fig. 4F). The EMSA results also confirmed that NAC family TF of SbNAC68 could directly bind to the 5′-ACGCAA-3′ motif, which was generally detected within the SBRGs promoter regions in sorghum (Fig. 7C). TFs usually co-expressed with SBRGs in sorghum grains, and ATAC-seq analysis further confirmed the existence of specific binding sites between TFs and related genes (Tables S13 and S14). Overall, transcriptional regulation plays a vital role in starch biosynthesis, and SbNAC68 is an important regulator for starch biosynthesis in sorghum grain.

5. Conclusions

We characterized sorghum starch biosynthesis, discovered transcriptional regulatory networks at different levels through a combination of multiple omics, and constructed a preliminary transcriptional regulatory model for sorghum grain starch biosynthesis. Starch accumulation in sorghum grains gradually increases, with peak sucrose contents at 10 DAP. Differences in gene expression occur during grain development, and co-expression occurs between TFs and SBRGs. A case study of SbNAC68 confirmed that the motifs contained in the peak of key enzymes of sorghum starch biosynthesis identified by ATAC-seq are also important binding sites for TFs. The summary profiling of multiple omics in the present study provides molecular insights into the transcriptional mechanisms regulating of sorghum starch biosynthesis.

Acknowledgements

This work was supported by the National Natural Science Foundation of China (32372076, 32001607), the Fundamental Research Funds for the Central Universities (SWU-XDJH202315). The authors would like to express their gratitude to EditSprings (https://www.editsprings.cn) for the expert linguistic services provided.

CRediT authorship contribution statement

Tingting Dai: Investigation. Hong Duan: Methodology, Investigation. Jing Li: Visualization, Investigation, Data curation, Formal analysis, Writing – original draft, Conceptualization. Chaoxian Liu: Resources, Methodology. Jiuguang Wang: Methodology. Xiangling Gong: Visualization, Investigation, Data curation, Formal analysis, Writing – original draft. Qianlin Xiao: Visualization, Conceptualization, Supervision, Project administration, Writing – review & editing, Resources, Formal analysis, Writing – original draft, Funding acquisition. Zhizhai Liu: Supervision, Investigation, Data curation, Resources, Software, Visualization, Project administration, Writing – review & editing. Min Yin: Investigation. Wei Wei: Investigation. Hameed Gul: Methodology. Yi Zheng: Investigation

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Appendix A. Supplementary data

Supplementary data to this article can be found online at https://doi.org/10.1016/j.pld.2025.11.006.