Transcriptome Analysis Reveals Potential Mechanism in Storage Protein Trafficking within Developing Grains of Common Wheat

Gluten proteins are the major storage protein fraction in the mature wheat grain. They are restricted to the starchy endosperm, which defines the viscoelastic properties of wheat dough. The synthesis of these storage proteins is controlled by the endoplasmic reticulum (ER) and is directed into the vacuole via the Golgi apparatus. In the present study, transcriptome analysis was used to explore the potential mechanism within critical stages of grain development of wheat cultivar “Shaannong 33” and its sister line used as the control (CK). Samples were collected at 10 DPA (days after anthesis), 14 DPA, 20 DPA, and 30 DPA for transcriptomic analysis. The comparative transcriptome analysis identified that a total of 18,875 genes were differentially expressed genes (DEGs) between grains of four groups “T10 vs. CK10, T14 vs. CK14, T20 vs. CK20, and T30 vs. CK30”, including 2824 up-regulated and 5423 down-regulated genes in T30 vs. CK30. Further, the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment highlighted the maximum number of genes regulating protein processing in the endoplasmic reticulum (ER) during grain enlargement stages (10–20 DPA). In addition, KEGG database analysis reported 1362 and 788 DEGs involved in translation, ribosomal structure, biogenesis, flavonoid biosynthesis pathway and intracellular trafficking, secretion, and vesicular transport through protein processing within ER pathway (ko04141). Notably, consistent with the higher expression of intercellular storage protein trafficking genes at the initial 10 DPA, there was relatively low expression at later stages. Expression levels of nine randomly selected genes were verified by qRT-PCR, which were consistent with the transcriptome data. These data suggested that the initial stages of “cell division” played a significant role in protein quality control within the ER, thus maintaining the protein quality characteristics at grain maturity. Furthermore, our data suggested that the protein synthesis, folding, and trafficking pathways directed by a different number of genes during the grain enlargement stage contributed to the observed high-quality characteristics of gluten protein in Shaannong 33 (Triticum aestivum L.).


Functional Annotation between Different Grain Development Stages
First, KEGG enrichment analysis was conducted by using DEGs at all four stages and libraries of all collected samples (T10, T14, T20, and T30 and CK10, CK14, CK20, and CK30), as shown in Figure 5 and Supplementary File S1 sheet 2. The maximum functional annotation of 315 DEGs in pathway id (map04141) is linked with protein processing within the ER; that in the secondary category is described as folding, sorting, and degradation of GSPs; 237 DEGs were observed in starch and sucrose metabolism in pathway id (map00500); 213 in pathway id map04075 were associated with plant hormone signal transduction. Further, the KEGG enrichment analysis was conducted by using all four groups (T10 vs. CK10, T14 vs. CK14, T20 vs. CK20, and T30 vs. CK30), provided in Figure 6 in detail, Supplementary File S1 (sheets 1-4). The rich factor with Padjust (0-0.5) showed phenylpropanoid biosynthesis and carotenoid biosynthesis, of which the pathway id was "map00940". There were 8, 3, 8, and 8 numbers of genes in association with pathway map00906 in all four groups of (T30 vs. CK30, T20 vs. CK20, T14 vs. CK14, and T10 vs. CK10). For further improvement of our knowledge about the total number of genes involved in grain development, function classification of clusters of orthologous groups (COGs) was conducted for all four groups. There was a distribution of 20 COG categories based on COG functional classification, shown in Figure 7. The cell part involved of 1200 genes in group "T30 vs. CK30", in comparison with >200 DEGs involved in protein processing of "T30 vs. CK30" of cellular components. Furthermore, a functional description is given in Supplementary File S1, sheet 3. A total of 1362 DEGs were reported in translation, ribosomal structure, and biogenesis in "T14 vs. CK14", in comparison with other groups. Consistent with functional analysis, we performed GO enrichment analysis of all four groups (T10 vs. CK10, T14 vs. CK14, T20 vs. CK20, and T30 vs. CK30), as shown in Figure 8 and Supplementary File S3, sheets 1-4. There were 20 categories of functional classification, comprised of negative regulation of the protein metabolic process, negative regulation of the cellular protein metabolism process, and so on, with further detail given in Figure 8 and Supplementary File S3.

Functional Annotation between Different Grain Development Stages
First, KEGG enrichment analysis was conducted by using DEGs at all four stages and libraries of all collected samples (T10, T14, T20, and T30 and CK10, CK14, CK20, and CK30), as shown in Figure 5 and Supplementary File S1 sheet 2. The maximum functional annotation of 315 DEGs in pathway id (map04141) is linked with protein processing within the ER; that in the secondary category is described as folding, sorting, and degradation of GSPs; 237 DEGs were observed in starch and sucrose metabolism in pathway id (map00500); 213 in pathway id map04075 were associated with plant hormone signal transduction. Further, the KEGG enrichment analysis was conducted by using all four groups (T10 vs. CK10, T14 vs. CK14, T20 vs. CK20, and T30 vs. CK30), provided in Figure 6 in detail, Supplementary File S1 (sheets 1-4). The rich factor with Padjust (0-0.5) showed phenylpropanoid biosynthesis and carotenoid biosynthesis, of which the pathway id was "map00940". There were 8, 3, 8, and 8 numbers of genes in association with pathway map00906 in all four groups of (T30 vs. CK30, T20 vs. CK20, T14 vs. CK14, and T10 vs. CK10). For further improvement of our knowledge about the total number of genes involved in grain development, function classi- in comparison with other groups. Consistent with functional analysis, we performed GO enrichment analysis of all four groups (T10 vs. CK10, T14 vs. CK14, T20 vs. CK20, and T30 vs. CK30), as shown in Figure 8 and Supplementary File S3, sheets 1-4. There were 20 categories of functional classification, comprised of negative regulation of the protein metabolic process, negative regulation of the cellular protein metabolism process, and so on, with further detail given in Figure 8 and Supplementary File S3. Figure 5. KEGG pathway enrichment analysis for 127 candidate genes. The x-axis shows the gene ratio; the y-axis corresponds to KEGG pathways. The dot color represents the corrected P value of <0.05, and the dot size represents the number of genes enriched in the reference pathway. The protein processing in endoplasmic reticulum, starch and sucrose metabolism, amino sugar and nucleotide sugar metabolism, glycolysis/gluconeogenesis, necroptosis, and pyruvate metabolism was significantly enriched, detail given in Supplementary File S1, sheet 1. Figure 5. KEGG pathway enrichment analysis for 127 candidate genes. The x-axis shows the gene ratio; the y-axis corresponds to KEGG pathways. The dot color represents the corrected P value of <0.05, and the dot size represents the number of genes enriched in the reference pathway. The protein processing in endoplasmic reticulum, starch and sucrose metabolism, amino sugar and nucleotide sugar metabolism, glycolysis/gluconeogenesis, necroptosis, and pyruvate metabolism was significantly enriched, detail given in Supplementary File S1, sheet 1.

Pathway Analysis of Protein Processing in Endoplasmic Reticulum (ER)
The synthesis, folding, and deposition of the gluten proteins take place within the endomembrane system of a developing grain. They fold in the lumen or membrane of the ER from where they are sorted toward their site of action [26,58,59]. ER protein quality control in the biological processes has long been a subject of intensive study. In the present study, KEGG pathway analysis identified 7,9,19, and 17 DEGs for the cellular amino acid metabolic process, protein metabolic process, Golgi vesicle transport, and protein serine/threonine kinase activity, respectively, in association with protein processing in the ER pathway (ko04141) (Figure 9), and for the protein export pathway (map03060) ( Figure  S1). The protein quality is controlled within the ER. The protein synthesis, folding, and trafficking are marked by the composition of N-glycan. Here, we observed that the folding sensor UGGT acts as an unusual molecular chaperone and covalently modifies (G1M9) the folding intermediate to correctly folded M9 ERMan1. The putative genes take part in recognition of folding and misfolded proteins within the ER; in addition, genes are listed in Table 3.

Pathway Analysis of Protein Processing in Endoplasmic Reticulum (ER)
The synthesis, folding, and deposition of the gluten proteins take place within the endomembrane system of a developing grain. They fold in the lumen or membrane of the ER from where they are sorted toward their site of action [26,58,59]. ER protein quality control in the biological processes has long been a subject of intensive study. In the present study, KEGG pathway analysis identified 7,9,19, and 17 DEGs for the cellular amino acid metabolic process, protein metabolic process, Golgi vesicle transport, and protein serine/threonine kinase activity, respectively, in association with protein processing in the ER pathway (ko04141) (Figure 9), and for the protein export pathway (map03060) ( Figure S1). The protein quality is controlled within the ER. The protein synthesis, folding, and trafficking are marked by the composition of N-glycan. Here, we observed that the folding sensor UGGT acts as an unusual molecular chaperone and covalently modifies (G1M9) the folding intermediate to correctly folded M9 ERMan1. The putative genes take part in recognition of folding and misfolded proteins within the ER; in addition, genes are listed in Table 3.

Analysis of Different Groups of Differently Expressed Genes
For comprehensive analysis, to understand the comparative difference between different groups in association with DEGs, which controlled the protein procession within the ER, all four constructed groups "T10 vs. CK10, T14 vs. CK14, T20 vs. CK20, and T30 vs. CK30" were analyzed. The similarity index of DEGs between all used groups is provided in Figure 10. There were 1838 genes reported in "T10 vs. CK10" in comparison with "T14 vs. CK14" were 4137. The minimum number of DEGs were seen in "T10 vs. CK10" in comparison with the three other groups. This suggested the primary involvement of the "cell division" phase of 10 DPA in protein quality control of the ER. Further, to identify the transcript level of DEGs in comparison with different developing stages of both used wheat samples, the transcript level of each gene in all four groups was compared and filtered with |log10 (TPM)|> = 1 and FDR < 0.001; for further detail, see Supplementary File S1, sheet 4.  Green points indicate down-regulated genes and red points indicate up-regulated genes in each contrast, detail given in Supplementary file S1, sheet 3 (Table S1).

Putative Genes Related to Protein Quality Control
To further investigate the role of ER-associated protein folding, sorting, and degradation, at different stages of grain development in wheat, the candidate genes related to cell wall organization or biogenesis, protein transporter activity, the structural constituent of ribosomes, protein folding, unfolded protein, intracellular protein transport, GTPase activity, ATPase activity, vesicle-mediated transport, N-glycan processing (ER), intracellular protein transport (COPI vesicle coat), intra-Golgi vesicle-mediated transport, mismatch repair, retrograde vesicle-mediated transport, Golgi to ER, Full -Glutenin, protein serine/threonine kinase activity, beta-amylase activity, low-molecular-weight glutenin (LMW) subunit, glycogen (starch) synthase activity, and cysteine-type endopeptidase activity were identified according to the GO annotation and local TBLAST search. Only the genes with more than a two-fold change in FPKM within the developing grain of all selected stages are enlisted, with details given in (Supplementary File S4, sheets 1-5).
A total of 219 and 175 genes were identified to be linked with protein folding, unfolding, a structural constituent of ribosomes (ER), and intracellular protein trafficking such as GTPase activity and ATPase activity, respectively, given in Supplementary File CK14, (C) Group T20 vs. CK20, and (D) Group T30 vs. CK30. Points outside the gray area indicate >100-fold-differences in expression. Green points indicate down-regulated genes and red points indicate up-regulated genes in each contrast, detail given in Supplementary file S1, sheet 3 (Table S1).

Putative Genes Related to Protein Quality Control
To further investigate the role of ER-associated protein folding, sorting, and degradation, at different stages of grain development in wheat, the candidate genes related to cell wall organization or biogenesis, protein transporter activity, the structural constituent of ribosomes, protein folding, unfolded protein, intracellular protein transport, GTPase activity, ATPase activity, vesicle-mediated transport, N-glycan processing (ER), intracellular protein transport (COPI vesicle coat), intra-Golgi vesicle-mediated transport, mismatch repair, retrograde vesicle-mediated transport, Golgi to ER, Full -Glutenin, protein serine/threonine kinase activity, beta-amylase activity, low-molecular-weight glutenin (LMW) subunit, glycogen (starch) synthase activity, and cysteine-type endopeptidase activity were identified according to the GO annotation and local TBLAST search. Only the genes with more than a two-fold change in FPKM within the developing grain of all selected stages are enlisted, with details given in (Supplementary File S4, sheets 1-5).
A total of 219 and 175 genes were identified to be linked with protein folding, unfolding, a structural constituent of ribosomes (ER), and intracellular protein trafficking such as GTPase activity and ATPase activity, respectively, given in Supplementary File S4, sheet 1. The genes regarding the defense mechanism are mainly involved in response to the biotic stimulus and disease resistance protein RPP13 (Supplementary File S4, sheet 3), and intracellular-trafficking-, secretion-, and vesicular-transport-related DEGs (Supplementary File S1, Sheet 3 (Table S1) & S4, sheet 2). The 84 genes identified here encode for grain quality characteristics such as grain hardness (protein serine/threonine kinase activity), starch synthesis, and viscoelasticity (cysteine-type endopeptidase activity) ( Supplementary File S4 and sheet 4). Further, all the expressed genes in all collected samples of Shaannong 33 and CK are shown in Figure 11, with detail given in Supplementary File S5. There were a total of 696 DEGs in all collected samples, from which TraesCS1A02G011300, TraesCS1B02G017000, and TraesCS1B02G276200 showed a two-fold transcript level in T30 samples. TraesCS4A02G316000 and TraesCS5D02G551600 reported an optimum two-fold transcript level in T10. In comparison, TraesCS4A02G316000 reported an optimum transcript level, with further detail shown in Figure 11, Supplementary File S5. The expression profile of all genes involved in protein trafficking and protein folding and unfolding were well consistent with FPKM values from different developing stages ( Figure 11). The expression profiles of all genes involved in protein trafficking and protein folding and unfolding were well consistent with FPKM values from different development stages ( Figure 11). TraesCS1A02G133100 was associated with GO: 0006457 (20 DPA) and T30 (30 DPA). This suggested intracellular protein quality control within the ER. GTPase activity has been approved as an important regulator in gluten quality characteristics, through controlling ER-to-Golgi vesicle trafficking [59][60][61][62][63]. TraesCS1A02G137213 linked with GO: 0043547 (positive regulation of GTPase activity) showed a higher expression at 10 DPA (cell division phase), which decreased at 30 DPA. Similarly, TraesCS1A02G131500 showed a higher GO: 0007264 (mall-GTPase-mediated signal transduction) expression at 10 DPA, which then decreased with exposure to increasing DPA. TraesCS1A02G060700 suggested a linkage with translation, ribosomal structure, and biogenesis expression decrease until 30 DPA. TraesCS1A02G092900 gene showed that GO: 0003735 (structural constituent of ribosome) expression increases with increasing development, a putative regulator during protein processing within the ER quality controlling process. TraesCS1A02G126800 showed that the GO: 0080163 (regulation of protein serine/threonine phosphatase activity) expression level was medium when exposed to grain development stages. Defense-mechanisms-related putative genes such as TraesCS1A02G081300 showed optimum expression at 30 DPA. A total of nine randomly selected genes were validated through RT-qPCR analysis ( Figure 12).

Discussion
In the developing endosperm of bread wheat (Tritium aestivum L.), seed storage proteins are produced on the rough Endoplasmic Reticulum (rER) and transported to protein bodies, specialized vacuoles for the GSP (Gluten) synthesis. The GSP is the main source of gluten protein synthesis within endosperm, governing its end-use value [27,28,[30][31][32][33][34][35][36]64]. However, the underlying mechanism of protein quality during ER-to-Golgi trafficking, specifically folding and unfolding, is still a challenging factor in molecular terms [65]. The present study conducted a compressive analysis by using a deep transcriptomic survey at four critical stages (10 DPA, 14 DPA, 20 DPA, and 30 DPA), to observe systematic changes and the potential mechanism associated within the developing grain of highquality characteristics of wheat cultivar Shaannong33 (Triticum aestivum L.), in association with its low-quality sister line (CK). The peak mitotic division within starchy endosperm occurs after 10 DPA in maize, after 12 in barley, and in wheat, it remains until 16 DPA [66][67][68][69][70]. A sum total of 125,729 DEGs were retrieved at all selected stages and 94,972 genes were classified in all four groups (T10 vs. CK10, T14 vs. CK14, T20 vs. CK20, and T30 vs. CK30) for further study. They were mainly involved in translation, ribosomal structure and biogenesis, intracellular trafficking, secretion, vesicular transport, defense mechanism, serine synthesis, cell wall modeling, carbohydrate metabolism, starch and sucrose metabolism, and many unknown functions (Table 3 and

Discussion
In the developing endosperm of bread wheat (Tritium aestivum L.), seed storage proteins are produced on the rough Endoplasmic Reticulum (rER) and transported to protein bodies, specialized vacuoles for the GSP (Gluten) synthesis. The GSP is the main source of gluten protein synthesis within endosperm, governing its end-use value [27,28,[30][31][32][33][34][35][36]64]. However, the underlying mechanism of protein quality during ER-to-Golgi trafficking, specifically folding and unfolding, is still a challenging factor in molecular terms [65]. The present study conducted a compressive analysis by using a deep transcriptomic survey at four critical stages (10 DPA, 14 DPA, 20 DPA, and 30 DPA), to observe systematic changes and the potential mechanism associated within the developing grain of high-quality characteristics of wheat cultivar Shaannong33 (Triticum aestivum L.), in association with its low-quality sister line (CK). The peak mitotic division within starchy endosperm occurs after 10 DPA in maize, after 12 in barley, and in wheat, it remains until 16 DPA [66][67][68][69][70]. A sum total of 125,729 DEGs were retrieved at all selected stages and 94,972 genes were classified in all four groups (T10 vs. CK10, T14 vs. CK14, T20 vs. CK20, and T30 vs. CK30) for further study. They were mainly involved in translation, ribosomal structure and biogenesis, intracellular trafficking, secretion, vesicular transport, defense mechanism, serine synthesis, cell wall modeling, carbohydrate metabolism, starch and sucrose metabolism, and many unknown functions (Table 3 and Figures 5-8), which were also widely studied in other plant species. A number of studies had confirmed that the rule of the effects of GSP trafficking between the ER to Golgi body is linked with grain processing quality characteristics, especially during grain filling [27,28,65,[70][71][72]. Here, we identified 315 DEGs in pathway id (map04141) that are linked with protein processing within the ER in the secondary category described as folding, sorting, and degradation, linked with GSPs degradation (Figure 9 and Table 3). The maximum translation, ribosomal structure, and biogenesis took place within T10 vs. CK10, T14 vs. CK14, and T20 vs. CK20 suggested optimum mitotic cell division, as previously reported by [68][69][70][71][72][73]. Further, we found that protein folding, sorting, and degradation-associated genes were enriched within 10 DPA-20 DPA. This might indicate the systematic role in quality control within ER to Golgi bodies, resulting in the enhancement of GSP quality improvement with storage protein activator (SPA). The SPA as a transcriptional regulator plays a significant role in GSP trafficking, ultimately defining dough viscoelasticity and grain hardness [57,73,74]. The 12 cysteine residues have a prominent effect on HMW-GS subunit "5 + 10" (more possible disulfide bonds in its chemical structure), leading to higher strength compared to the HMW-GS subunit (2 + 12), which possesses 11 cysteines [74,75]. Therefore, we suggested the significant number of genes in association with serine/threonine kinase activity, enlisted in Supplementary File S4. TraesCS1A02G148428 showed a 63.71-fold expression level at 10 DPA, which was recognized in GO: 0006535 (cysteine biosynthetic process from serine), suggesting the putative role in increasing tensile resistance of the used wheat cultivar "Shaanong 33". It further signifies the role of GSP trafficking to the ER. Furthermore, the specific stage of grain development could be helpful for further investigation and could open a new venue toward an ER protein quality control mechanism.
The ER is the largest organelle in the cell and is a major site of protein synthesis and transport, protein folding, calcium storage, and carbohydrate metabolism [58,59,76,77]. The maximum activity of starchy endosperm development took place until 16 DPA [31,56,78,79]. The transcript driving the progression of endosperm development was differently expressed at different stages of harvest [9,76,79], as shown in Figure 11. The 60-80% GSP consists of gluten, developed within starchy endosperm. It determines the dough cohesiveness and viscoelasticity [9,74,80,81], which is a highly recommended wheat for bread and Chinese yellow alkaline noodles making. The extensive role of GSP highly endorses it for further extensive studies. In addition, the dynamic role of quantitative genetic variation in GSP composition has been reported for wheat [3,9,39,40,82]. In this study, we detected 315 DEGs encoding GSPs translocation in association with map0414, shown in Figure 9, Supplementary File S4 (sheets 1-4). They trigger protein processing pathways to maintain a homeostasis of the total amount of protein per grain [75,82,83]. The Rab GTPase is best characterized in Arabidopsis thaliana for its key involvement in specifically the trafficking mechanism. It works as a molecular switch to drive the transport of vesicles between membranous compartments. The motor proteins within the membrane regulate vesicle and compartment motility [28][29][30]33,61,62,84]. Hence, different Rab sub-clades have been classified in protein vesicle bodies trafficking to the Golgi apparatus and vacuole [85]. Consistent with the previous successful work, the present study also reported 177 DEGs in protein vesicle bodies trafficking to Golgi apparatus (Supplementary File S4, sheet 2), from which TraesCS1A02G223800 showed a 651.31-fold expression at 14 DPA, linked with GO: 0043547 functional association of positive regulation of GTPase activity. Further, at 14 DPA, the second major transition stage is where starch and seed storage proteins accumulate within cells to make semi-solid endosperm [56,79,86]. This suggested the role of storage protein activator SPA, in the key functional mechanism to maintain a homeostasis supply of nutrients between the maternal cells and the starchy endosperm, and linked with the serine-to-cysteine biosynthesis process. Rab GTPase has been approved in different quality characteristics of mango, and others [63,87]. Therefore, further investigation is needed to improve our knowledge. However, the current study could lead to a specific stage and candidate gene in association with protein processing within the developing grain.

Plant Materials
Shaannong 33 is a newly released high-quality, strong gluten wheat variety that has been approved in Shaanxi Province and introduced in Henan Province, China [88,89]. A summary and details of Shaannong 33 are shown in Table 4, and the sister line of Shaannong 33, used as the control (CK), reported with low-quality traits, given in Table 1, represents three-year field trails. The plant material was planted in a 2 m row length, the plant-toplant distance was 0.25 m, and each line was comprised of 40 plants, in three biological replicates [88][89][90][91], grown at Northwest Agriculture and Forestry University (Field station) under non-stressed conditions from October 2019 to 2020. The main individual spikes were tagged at the first day of anthesis. The tagged spikes were harvested at four main seed developmental stages, i.e., 10 DPA, 14 DPA, 20 DPA, and 30 DPA, immediately frozen in liquid nitrogen, and stored at −80 • C. The embryos of these samples were used for RNA extraction to isolate the entire RNA expressed in both used planted materials.

Library Preparation and Illumina Hiseq4000 Sequencing
There were 24 RNA-seq transcriptome libraries that were prepared following the TruSeqTM RNA sample preparation Kit from Illumina (San Diego, CA, USA) using 5 µg of total RNA. First, messenger RNA was isolated according to the polyA selection method by oligo (dT) beads and then fragmented by fragmentation buffer. Secondly, doublestranded cDNA was synthesized using a Superscript double-stranded cDNA synthesis kit (Invitrogen, Carlsbad, CA, USA) with random hexamer primers (Illumina). Then, the synthesized cDNA was subjected to end-repair, phosphorylation, and "A" base addition according to Illumina's library construction protocol. Libraries were size-selected for cDNA target fragments of 200-300 bp on 2% low-range ultra-Agarose followed by PCR amplified using Phusion DNA polymerase (NEB) for 15 PCR cycles. After quantification by TBS380, the paired-end RNA-seq sequencing library was sequenced with the Illumina HiSeq 4000 (2 × 150 bp read length) [91,92].

Read Mapping
The raw paired-end reads were trimmed and quality-controlled by 11 October 2021 Se-qPrep (https://github.com/justjohn/SeqPrep) and 11 October 2021 Sickle (https://github.com/ najoshi/sickle) with default parameters. Then, clean reads were separately aligned to reference Triticum_aestivum version: iwgsc_refseqv1.0 source: http://www.wheatgenome.org/News/ Latest-news/All-IWGSC-reference-sequence-resources-now-publicly-available-at-URGI with the orientation mode using TopHat 2.0.12 (Johns Hopkins University, Baltimore, Maryland) (http://tophat.cbcb.umd.edu/.version2.0.0) [93] software. The mapping of a criterion of bowtie was as follows: sequencing reads should be uniquely mated to the genome, allowing up to 2 mismatches, without insertions or deletions. Then, the region of the gene was expanded following depths of sites, and the operon was obtained. In addition, the whole genome was spilt into multiple 15 kbp windows that share 5 kbp. New transcribed regions were defined as more than 2 consecutive windows without overlapped regions of the gene, where there were at least 2 reads mapped per window in the same orientation.

Differential Expression Analysis and Functional Enrichment
A total of 24 RNA samples from Shaannong 33 and CK samples at four spike developmental stages (10,14,20, and 30 DAA) were sent to Shanghai Majorbio Bio-Pharm Technology Co., Ltd. (Shanghai, China) using paired-end sequencing on an Illumina HiSeq PE150 platform Shanghai Meiji Biomedical Technology (Shanghai, China) for RNA-sequencing and transcriptome assembly using previously described protocols [91,92]. Fragments per kilo-base of transcript per million mapped reads (FPKM) values were used to calculate transcript abundance in Shaannong 33 and CK collected samples. 21 October 2021 RSEM (http://dewelab.biostat.wisc.edu/rsem/) [94] was used to quantify gene abundances. A summary of the sequencing statistics for reads and bases obtained in each sample and statistics for the mapping of clean reads from each sample to the reference genome are given in Supplementary File S1. Further, the significance of differentially expressed genes (DEGs) was analyzed using R Statistical package software at 30 October 2021, EdgeR (empirical analysis digital gene expression in R, http://www.bioconductor.org/packages/2.12/bioc/ html/edeR.html), a published method [95]. The biological processes associated with each gene in each sample were evaluated using Gene Ontology (GO) analysis 1 November 2021 (http://www.geneontology.org/). KEGG pathways with corrected p-values < 0.01 were considered to be significantly enriched (http://www.genome.jp/kegg/) [96].

New Isoforms Prediction
The TOpHat-Cufflinks pipeline was used to predict gene isoforms from our RNA-seq data. In 1 January 2022 TopHat (http://tophat.cbcb.umd.edu/, version 2.0.0) [93], the option "min-isofom-fraction" was disabled; instead, "coverage-search", "butterfly-search", and "microexon-search" were used. The expected fragment length was set to 200 bp and the "small-anchor-fraction" was set to 0.08, which requires at least 8 bp on each side of an exon junction for our 100-bp RNA-seq data. Cuffcompare was used to compare and merge the reference annotation and the isoform predictions.

Alternative Splice Events Identification
All the alternative splice events that occurred in our sample were identified by using the recently released program Multivariate analysis of 1 December 2021 Transcript Splicing (MATS, http://rnaseq-mats.sourceforage.net/) [97]. Only the isoforms that were similar to the reference or comprised novel splice junctions were considered, and the splicing differences were detected as exon inclusion, exclusion, alternative 5', 3', and intron retention events.

Gene Expression Level Validation by qRT-PCR
Total RNA was isolated by using a total RNA kit (TIAGEN, Beijing, China) following the manufacturer's instructions, followed by [88,89]. For the complementary DNA (cDNA) synthesis, 2 mg of total RNA was reverse-transcribed using a Prime-Script TM II First-Strand cDNA synthesis Kit (TakaRa, Dalian, China). The quantitative real-time PCR (qRT-PCR) analysis was performed with a Light Cycle 480 II system (Roche, Basel, Switzerland) using an SYNR Premix Ex TaqTM kit (TakaRa, Dalian, China). The UBI-eq gene was used as an internal control. The primer sequences used for the PCR are given in (Supplementary File S5, sheet 2). Three replicates were made for three separate RNA extracts from the three samples.

Conclusions
The accumulation of protein bodies is taken either directly within the lumen of the ER or through the Golgi apparatus into the vacuole [27][28][29][30][31][59][60][61][62][63]70,98]. Here, we proposed a model through combining our results and previous successful works, which could be used to test how trafficking of GSP within the developing grain improves the wheat grain endproduct quality (Figure 13). This proposal begins with the expression of the translation-, ribosomal-structure-and biogenesis-, and protein-folding-and unfolding-related DEGs. The TFs bind with the promoter region of a suite of genes, improving storage proteins and protein synthesis regulator proteins in wheat grains. Thus, the composition changes within the ribosome lead to a higher demand for protein synthesis and more molecular chaperones assisting nascent proteins to fold properly and be stored. The abundances of chaperones and other storage proteins into the protein storage vacuoles with enhanced stability led to a broad increase in the abundance of HMW-GSs, further enhancing storage protein synthesis and deposition, and leading to a significant improvement in the dough quality with changes in gluten composition. Therefore, this study suggested the identification of putative DEGs that could contribute to enhanced protein content, at different stages of DPA in wheat grains. synthesis, 2 mg of total RNA was reverse-transcribed using a Prime-Script TM II First-Strand cDNA synthesis Kit (TakaRa, Dalian, China). The quantitative real-time PCR (qRT-PCR) analysis was performed with a Light Cycle 480 II system (Roche, Basel, Switzerland) using an SYNR Premix Ex TaqTM kit (TakaRa, Dalian, China). The UBI-eq gene was used as an internal control. The primer sequences used for the PCR are given in (Supplementary File S5, sheet 2). Three replicates were made for three separate RNA extracts from the three samples.

Conclusions
The accumulation of protein bodies is taken either directly within the lumen of the ER or through the Golgi apparatus into the vacuole [27][28][29][30][31][59][60][61][62][63]70,98]. Here, we proposed a model through combining our results and previous successful works, which could be used to test how trafficking of GSP within the developing grain improves the wheat grain end-product quality ( Figure 13). This proposal begins with the expression of the translation-, ribosomal-structure-and biogenesis-, and protein-folding-and unfolding-related DEGs. The TFs bind with the promoter region of a suite of genes, improving storage proteins and protein synthesis regulator proteins in wheat grains. Thus, the composition changes within the ribosome lead to a higher demand for protein synthesis and more molecular chaperones assisting nascent proteins to fold properly and be stored. The abundances of chaperones and other storage proteins into the protein storage vacuoles with enhanced stability led to a broad increase in the abundance of HMW-GSs, further enhancing storage protein synthesis and deposition, and leading to a significant improvement in the dough quality with changes in gluten composition. Therefore, this study suggested the identification of putative DEGs that could contribute to enhanced protein content, at different stages of DPA in wheat grains. Figure 13. Potential mechanism of protein synthesis and trafficking from ER to Golgi body and ultimate linkage with processing quality characteristics in wheat grain. Figure 13. Potential mechanism of protein synthesis and trafficking from ER to Golgi body and ultimate linkage with processing quality characteristics in wheat grain.