Compensatory Modulation of Seed Storage Protein Synthesis and Alteration of Starch Accumulation by Selective Editing of 13 kDa Prolamin Genes by CRISPR-Cas9 in Rice

Rice prolamins are categorized into three groups by molecular size (10, 13, or 16 kDa), while the 13 kDa prolamins are assigned to four subgroups (Pro13a-I, Pro13a-II, Pro13b-I, and Pro13b-II) based on cysteine residue content. Since lowering prolamin content in rice is essential to minimize indigestion and allergy risks, we generated four knockout lines using CRISPR-Cas9, which selectively reduced the expression of a specific subgroup of the 13 kDa prolamins. These four mutant rice lines also showed the compensatory expression of glutelins and non-targeted prolamins and were accompanied by low grain weight, altered starch content, and atypically-shaped starch granules and protein bodies. Transcriptome analysis identified 746 differentially expressed genes associated with 13 kDa prolamins during development. Correlation analysis revealed negative associations between genes in Pro13a-I and those in Pro13a-II and Pro13b-I/II subgroups. Furthermore, alterations in the transcription levels of 9 ER stress and 17 transcription factor genes were also observed in mutant rice lines with suppressed expression of 13 kDa prolamin. Our results provide profound insight into the functional role of 13 kDa rice prolamins in the regulatory mechanisms underlying rice seed development, suggesting their promising potential application to improve nutritional and immunological value.

SSPs are synthesized at the rough endoplasmic reticulum (ER), translocated to the ER lumen, and then transferred to other intracellular compartments of the plant endomem-brane system [9,10].Prolamins accumulate in the ER-derived protein bodies (PB) I (PB-I), whereas globulins and glutelins accumulate in PB-II structures before subsequently moving to protein storage vacuoles (PSVs) [10,11].The PB-I particle has a layered structure, with the 10 kDa prolamin in the core (innermost) layer, the Pro13b-I subgroup in the layer surrounding the core, 16 kDa and Pro13a-I/II subgroups in the middle layer, and Pro13b-II subgroup prolamins in the outer-most layer [8].
Rice prolamins are deficient in lysine and methionine, two of the nine amino acids considered to be essential in the human diet [2,16,17].In addition, prolamins are considered indigestible and nutritionally inferior to glutelins [18,19].Several studies demonstrated that rice SSPs modulate glutathione metabolism, mitigate oxidative damage to lipids and proteins, and provoke antioxidative responses against hypercholesterolemia, although glutelins appeared to be superior to prolamins as antioxidants [19][20][21].Moreover, prolamins found in grains such as wheat, barley, and rye are known to cause allergic reactions such as celiac disease and IgE-mediated food allergies [22][23][24].It is interesting that amino acid sequence analysis has revealed a close proximity between prolamins and protein allergens (RA16/17) in rice seeds [25,26].Furthermore, previous studies have shown that the production of recombinant proteins is more efficient in mutant rice lacking 13 kDa prolamin (13 kDa Pro-less rice) than in wild-type (WT) control rice [27,28].Consequently, there is significant interest in rice lines engineered with low prolamin content, based on the idea that such mutant rice lines will facilitate the production and distribution of recombinant and/or therapeutic proteins and that the grain of such rice lines will be nutritionally superior for humans.Furthermore, engineered mutant rice lines that lack specific prolamin genes are excellent experimental tools for understanding how prolamins and other SSPs influence rice grain development, composition, and properties.Such knowledge will improve our understanding of rice biology and will inform future efforts to improve rice as an agriculturally valuable plant species.
Many published studies demonstrate that the suppression of 13 kDa prolamin expression using RNA-silencing [4,27] and RNA-interference [2,29] was accompanied by compensatory changes in the expression of other SSPs as well as ER stress and the expression of stress-related proteins.ER stress was associated with increased expression of protein folding-related genes, such as binding protein (BiP), protein disulfide isomerase (PDI), and calnexin (CNX) [30][31][32][33].BiP, a key ER chaperone that belongs to the HSP70 family, interacts with and promotes the proper folding of prolamins [34,35].PDI facilitates disulfide bond formation and rearrangement in glutelins and prolamins to promote maturation and the proper folding of SSPs during seed maturation.CNX binds to unfolded glycoproteins, preventing misfolded proteins from being transported from the ER to the Golgi apparatus [25,32,36].As mentioned above, transgenic rice lacking 13 kDa prolamins has increased levels of glutelin, globulin, BiP, and PDI [4].Similarly, in another study, RNAi-treated transgenic rice with low levels of 13 kDa prolamin showed compensatory accumulation of 10 kDa prolamin and glutelins, up-regulation of chaperone proteins such as BiP and PDI, and abnormal PB-I bodies [2].In addition, GPGb-RNAi knock-down (targeting glutelins, prolamins, and globulins) generated a rice line with low levels of glutelin A (GluA), cysteine-rich 13 kDa prolamin, and α-globulin, with higher accumulation of mRNA and protein for glutelin B (GluB) and chaperones (BiP1 and PDIL1-1) [29].A few studies targeted multiple glutelin genes by CRISPR-Cas9 and produced mutant rice lines with significantly reduced glutelin content accompanied by an alteration of starch structure [37,38].We previously suggested that rice glutelin genes may interact with genes that regulate the synthesis of starch and other seed storage proteins and modulate their expressions via post-transcriptional and translational mechanisms [38].However, the editing of the prolamin gene by CRISPR-Cas9 has never been reported.In this study, CRISPR-Cas9 was employed to target the Pro13a-I and Pro13b-I/II subgroups of 13 kDa prolamin genes and generated four mutant lines in rice, and their seeds were characterized in order to elucidate the role of specific 13 kDa prolamin genes in seed development and how deletion of 13 kDa prolamin affects seed characteristics, starch composition, and the organization of protein storage compartments during endosperm maturation.Our results indicated compensatory changes in the rice proteome and transcriptome and revealed alterations affecting ER stress response, transcription factor, and starch synthesis in mutant rice with a selectively edited 13 kDa prolamin gene.

Design of sgRNA for Targeting 13 kDa Prolamin Genes Using CRISPR-Cas9 Technology
There are four genes encoding 10 kDa prolamins, 28 genes encoding 13 kDa prolamins, and two genes encoding 16 kDa prolamins in the rice genome.The 13 kDa prolamin genes are divided into four subgroups (Pro13a-I, Pro13a-II, Pro13b-I, and Pro13b-II) based on cysteine content and sequence similarity (Figure 1A).The nucleotide sequence similarity between Pro13a.1 and Pro13a.2 of the Pro13a-I subgroup is 92.03%, while these two genes share 75.33 to 78.17% similarity with the four genes in the Pro13a-II subgroup (Pro13a.3~6),which are highly similar to each other (97.93 to 99.68%).Lastly, the Pro13b-II subgroup includes 18 prolamin genes, which demonstrate 87.59 to 100% sequence similarity to each other.The high similarity between 13 kDa prolamin genes in each subgroup indicates that multiple prolamin genes can be edited using one sgRNA.We designed two sgRNAs, as shown in Figure S1 and Table S1.sgRNA-pro13a targets two genes in the Pro13a-I subgroup, and sgRNA-pro13b targets 17 genes in the Pro13b-I/II subgroups.An in vitro assay was used to demonstrate that the sgRNAs efficiently and specifically induce cleavage of their target genes (Figure S2).Subsequently, sgRNA-pro13a and sgRNA-pro13b were cloned into a binary vector for the CRISPR-Cas9 system (Figure 1B and Figure S3).Each pCAMBIA-Cas9-sgRNA binary vector was then introduced into the Korean rice cv.Ilmi (Oryza sativa subsp.Japonica) via Agrobacterium-mediated transformation to generate 13 kDa prolamin-knockout rice.

Generating and Characterizing 13 kDa Prolamin-Knockout Rice
PCR analysis was performed using SpCas9 and HygR primers specific to the pCAMBIA-Cas9 vector to determine the presence of the binary vector transformed in the T 0 lines.Results showed that 8 and 10 T 0 transgenic rice plants were transformed, with editing genes inserted in the target loci of the Pro13a-I and Pro13b-I/II subgroups, respectively (Figure S4B).Next-generation sequencing was carried out to identify mutations in target genes in two T 0 transgenic plants (1a and 2a) in the Pro13a-I subgroup and two T 0 transgenic plants (4b and 8b) in the Pro13b-I/II subgroups (Figure S4A and Table S2).The results showed that genes in the Pro13a-I subgroup were mutated in the 1a and 2a T 0 plants at frequencies of 51.3 and 28.5%, respectively, and genes in the Pro13b-I/II subgroups were mutated in the 4b and 8b T 0 plants at frequencies of 4.8 and 12.6%, respectively.

Generating and Characterizing 13 kDa Prolamin-Knockout Rice
PCR analysis was performed using SpCas9 and HygR primers specific to the pCAM-BIA-Cas9 vector to determine the presence of the binary vector transformed in the T0 lines.Results showed that 8 and 10 T0 transgenic rice plants were transformed, with editing genes inserted in the target loci of the Pro13a-I and Pro13b-I/II subgroups, respectively (Figure S4B).Next-generation sequencing was carried out to identify mutations in target genes in two T0 transgenic plants (1a and 2a) in the Pro13a-I subgroup and two T0 transgenic plants (4b and 8b) in the Pro13b-I/II subgroups (Figure S4A and Table S2).The results showed that genes in the Pro13a-I subgroup were mutated in the 1a and 2a T0 plants at frequencies of 51.3 and 28.5%, respectively, and genes in the Pro13b-I/II subgroups were mutated in the 4b and 8b T0 plants at frequencies of 4.8 and 12.6%, respectively.

Target Gene Expression at Transcriptional and Translational Levels
The expression of prolamin-encoding transcripts was compared in the imma seeds of 13 kDa prolamin-knockout T2 and WT control plants 2 weeks after flowerin WAF).The results (Figure 2B) show that transcripts of mutant Pro13a.1 and Pro13a.2ge were remarkably decreased in 1a-8-1 plants; transcripts of mutant Pro13a.1 gene w strongly reduced in 2a-2-1 plants; transcripts of mutant Pro13b.1 and Pro13b.13genes w decreased in 4b-9-7 plants; and transcripts of mutant Pro13b.3were decreased in 8b plants.These results are consistent with expectations based on the targeted mutations troduced into these four mutant plant lines.At the translational level, SSP content analyzed in the mature seeds of T2 plants by SDS-PAGE and western blot.SDS-PA analysis did not show dramatic differences between SSP content/expression in the f mutant plant lines and WT control plants.However, western blot analysis showed prolamins in Pro13a-I/II subgroups were strongly and weakly decreased in mature se of 1a-8-1 and 2a-2-1 plants, respectively, due to the number of mutated genes.Simila Total SSPs (5 µg) were separated on a gradient SDS-PAGE gel (10-17.5%).The resolved SSPs were transferred onto a polyvinylidene fluoride (PVDF) membrane, and the membrane was incubated with anti-Pro13a-I or anti-Pro13b-I antibodies, as indicated.S/M: Size marker.

Target Gene Expression at Transcriptional and Translational Levels
The expression of prolamin-encoding transcripts was compared in the immature seeds of 13 kDa prolamin-knockout T 2 and WT control plants 2 weeks after flowering (2 WAF).The results (Figure 2B) show that transcripts of mutant Pro13a.1 and Pro13a.2genes were remarkably decreased in 1a-8-1 plants; transcripts of mutant Pro13a.1 gene were strongly reduced in 2a-2-1 plants; transcripts of mutant Pro13b.1 and Pro13b.13genes were decreased in 4b-9-7 plants; and transcripts of mutant Pro13b.3were decreased in 8b-3-9 plants.These results are consistent with expectations based on the targeted mutations introduced into these four mutant plant lines.At the translational level, SSP content was analyzed in the mature seeds of T 2 plants by SDS-PAGE and western blot.SDS-PAGE analysis did not show dramatic differences between SSP content/expression in the four mutant plant lines and WT control plants.However, western blot analysis showed that prolamins in Pro13a-I/II subgroups were strongly and weakly decreased in mature seeds of 1a-8-1 and 2a-2-1 plants, respectively, due to the number of mutated genes.Similarly, expression of Pro13b-I/II subgroups were strongly and weakly decreased in mature seeds of 4b-9-7 and 8b-3-9 plants, respectively (Figure 2C).

Protein Body Formation and SSP Composition in 13 kDa Prolamin-Knockout Lines
Transmission electron microscopy (TEM) was used to investigate PB-I formation and subcellular structures in the immature seeds collected (two WAF) from WT and 13 kDa prolamin-knockout T 2 plants.In WT, PB-I structures form by the sequential accumulation of concentric layers composed of Pro10, Pro13b-I, Pro13a-I/II, Pro16, and 13b-II prolamins outward from the core, thus exhibiting a distinct pattern of electron-dense spherical layers (Figure 5A).In contrast, in pro13a.1/13a.2-knockout(1a-8-1) and pro13a.1-knockout(2a-2-1) plants, PB-I particles were smaller, and the layered structure was entirely lacking.On the other hand, in pro13b.1/13b.13-knockout(4b-9-7) and pro13b.3-knockout(8b-3-9) plants, PB-I particles stained more darkly in the outer layers and core than WT PB-I particles.This indicates that the 13 kDa prolamins are critical for forming layered PB-I particles in the seed endosperm of WT rice.

Protein Body Formation and SSP Composition in 13 kDa Prolamin-Knockout Lines
Transmission electron microscopy (TEM) was used to investigate PB-I formation and subcellular structures in the immature seeds collected (two WAF) from WT and 13 kDa prolamin-knockout T2 plants.In WT, PB-I structures form by the sequential accumulation of concentric layers composed of Pro10, Pro13b-I, Pro13a-I/II, Pro16, and 13b-II prolamins outward from the core, thus exhibiting a distinct pattern of electron-dense spherical layers (Figure 5A).In contrast, in pro13a.1/13a.2-knockout(1a-8-1) and pro13a.1-knockout(2a-2-1) plants, PB-I particles were smaller, and the layered structure was entirely lacking.On the other hand, in pro13b.1/13b.13-knockout(4b-9-7) and pro13b.3-knockout(8b-3-9) plants, PB-I particles stained more darkly in the outer layers and core than WT PB-I particles.This indicates that the 13 kDa prolamins are critical for forming layered PB-I particles in the seed endosperm of WT rice.Furthermore, the total content of SSP per grain was higher in 1a-8-1 plant seeds than in WT, while the total SSP content per grain was similar in WT and the other mutant plant Furthermore, the total content of SSP per grain was higher in 1a-8-1 plant seeds than in WT, while the total SSP content per grain was similar in WT and the other mutant plant lines (Figure 5B).In contrast, prolamins (all sizes) represented 15.8, 15.3, 14.8, and 17.5% of all SSPs in 1a-8-1, 2a-2-1, 4b-9-7, and 8b-3-9 mutant plants, respectively, while prolamins represented 20.0% of WT SSPs (Figure 5C).Conversely, the percent content of glutelins increased significantly from 60.2% in WT to 73.9, 70.0, 69.8, and 65.8% in 1a-8-1, 2a-2-1, 4b-9-7, and 8b-3-9 mutant plant lines, respectively.These results indicate that suppression of expression of certain SSPs was accompanied by compensatory increases in the expression of other SSPs; more specifically, down-regulation of genes in the Pro13a-I subgroup was accompanied by up-regulation of genes in the Pro13a-II, Pro13b-I/II, GluB, and GluD SSP subgroups.Similarly, decreased transcription of Pro13b-I subgroup prolamins was accompanied by increased transcription of Pro13a-II, Pro10, and GluB SSPs.Nevertheless, down-regulation of genes in Pro13b-I/II subgroups was not accompanied by compensatory changes in expression of SSP genes in other groups (Figure 5D).

Correlation Network Involving DEGs and 13 kDa Prolamins
To investigate patterns of transcriptional change among prolamin genes in the immature seeds of WT, 1a-8-1, 2a-2-1, 4b-9-7, and 8b-3-9 plants, Pearson correlation analysis was performed among 927 DEGs and prolamin genes.This analysis identified 746 transcripts with a correlation coefficient ≥|0.7| (Table S6).We then constructed a correlation network based on GO enrichment categories such as RNA processing (BIN27.1),RNA transcription and RNA regulation of transcription (BIN27.2and 27.3), protein synthesis ribosome (BIN29.2.1 and 29.2.6), protein synthesis transfer RNA (BIN29.2.7), protein synthesis elongation (BIN29.2.4), protein folding (BIN29.6),protein degradation (BIN29.5),protein posttranslational modification (BIN29.4),protein targeting (BIN29.3),stress (BIN20), and transport (BIN34) (Figure 7).Furthermore, the results revealed that Pro13a.1 and Pro13a.2genes in the Pro13a-I subgroup are negatively correlated with genes in the Pro13a-II and Pro13b-I/II subgroups.Genes in the network formed six clusters with distinct patterns of correlation with prolamin genes, as follows: Cluster 1 shows positive correlation with genes in the Pro13a-II and Pro13b-I/II subgroups; cluster 2 demonstrates an inverse association with cluster 1, displaying a negative correlation with genes in the Pro13a-II and Pro13b-I/II subgroups; cluster 3 exhibits a positive correlation with the Pro13a.1 gene but a negative correlation with genes in the Pro13a-II and Pro13b-I/II subgroups, whereas its opposite is cluster 4; similarly, cluster 5 shows a positive correlation with Pro13a.2 gene and a negative correlation with genes in the Pro13a-II and Pro13b-I/II subgroups, with cluster 6 being its opposite (details in Table S7).The network reveals which genes are correlated with each prolamin gene.Thus, the results suggest that transcriptional regulatory mechanisms may exist and may explain how and why altered expression of some prolamins is accompanied by compensatory changes in expression of other prolamins and other network genes.

Correlation Network Involving DEGs and 13 kDa Prolamins
To investigate patterns of transcriptional change among prolamin genes in the im mature seeds of WT, 1a-8-1, 2a-2-1, 4b-9-7, and 8b-3-9 plants, Pearson correlation analys was performed among 927 DEGs and prolamin genes.This analysis identified 746 tran scripts with a correlation coefficient ≥|0.7| (Table S6).We then constructed a correlatio network based on GO enrichment categories such as RNA processing (BIN27.1),RNA transcription and RNA regulation of transcription (BIN27.2and 27.3), protein synthes ribosome (BIN29.2.1 and 29.2.6), protein synthesis transfer RNA (BIN29.2.7), protein syn thesis elongation (BIN29.2.4), protein folding (BIN29.6),protein degradation (BIN29.5 protein posttranslational modification (BIN29.4),protein targeting (BIN29.3),stres (BIN20), and transport (BIN34) (Figure 7).Furthermore, the results revealed that Pro13a.and Pro13a.2genes in the Pro13a-I subgroup are negatively correlated with genes in th Pro13a-II and Pro13b-I/II subgroups.Genes in the network formed six clusters with dis tinct patterns of correlation with prolamin genes, as follows: Cluster 1 shows positive co relation with genes in the Pro13a-II and Pro13b-I/II subgroups; cluster 2 demonstrates a inverse association with cluster 1, displaying a negative correlation with genes in th Pro13a-II and Pro13b-I/II subgroups; cluster 3 exhibits a positive correlation with th Pro13a.1 gene but a negative correlation with genes in the Pro13a-II and Pro13b-I/II sub groups, whereas its opposite is cluster 4; similarly, cluster 5 shows a positive correlatio with Pro13a.2 gene and a negative correlation with genes in the Pro13a-II and Pro13b-I/ subgroups, with cluster 6 being its opposite (details in Table S7).The network revea which genes are correlated with each prolamin gene.Thus, the results suggest that tran scriptional regulatory mechanisms may exist and may explain how and why altered ex pression of some prolamins is accompanied by compensatory changes in expression o other prolamins and other network genes.
Starch synthesis is initiated by converting glucose-1-phosphate to ADP-glucose by ADPglucose pyrophosphorylase (AGPase).ADP-glucose is used for amylose synthesis by granulebound starch synthases (GBSSs) and for amylopectin synthesis via a series of enzymatic reactions such as soluble starch synthases (SSs), branching enzymes (BEs), and debranching enzymes (DBEs) [39].In GPGb-RNAi and Glu-knockout lines, it has been presented that reduced starch content and altered starch composition are associated with a significant decrease in some starch synthesis genes at the transcriptional level [29,38].In this study, the expression of SSI and BEIIa genes is significantly lower in the immature seeds of the four 13 kDa prolamin-knockout mutant lines, which could explain the reduced total starch content.On the other hand, the up-regulation of the GBSSI gene in the immature seeds of the 4b-9-7 (pro13b.1/13b.13-knockout)line and its down-regulation in the immature seeds of other mutant lines are consistent with higher amylose content in starch granules in 4b-9-7 plants and lower amylose content in the other mutant plants.Moreover, rod/filamentous starch granules have been reported previously in grains with high amylose content [40,41], which is consistent with the idea that rod/filamentous starch granules observed in 4b-9-7 plants are caused by starch with increased amylose content.

Abnormal PB-I Morphology in 13 kDa Prolamin-Knockout Plants
In WT rice, PB-I particles form by sequential deposition of concentric circular layers of specific SSPs, where 10 kDa prolamin (Pro10) accumulates in the core layer, followed sequentially by cys-poor 13 kDa prolamin (Pro13b-I), cys-rich 13 kDa prolamin (Pro13a-I/II), 16 kDa prolamin (Pro16), and an outermost layer of cys-poor 13 kDa prolamin (Pro13b-II) [8,42].It can be asserted that the specific interactions among rice prolamins play a crucial role in the formation of PB-I structures.Cys-rich prolamins (Pro10, Pro13a, and Pro16) can form intermolecular disulfide bonds in the layer surrounding the central core, while cys-poor 13 kDa prolamins (Pro13b-I/II) may interact with cys-rich prolamins through hydrophobic interactions [42,43].It has been reported that the composition of prolamins in PB-I has a profound impact on its layered structure [44,45].Previous studies reported aberrant PB-I morphology in rice seeds from 13 kDa-prolamin RNAi-targeted transgenic plants.Specifically, seeds with reduced expression of RM1 (former name of Pro13a.2 in Pro13a-I subgroup), RM2 (former name of Pro13b.2 in the Pro13b-I subgroup), RM4 (former name of Pro13b.11/12/13 in the Pro13b-II subgroup), and RM9 (former name of Pro13a.4 in the Pro13a-II subgroup) displayed PB-Is with cracks and jagged periphery structures [4].In addition, small PB-Is with no lamella structure were observed when expression of the Pro13a-I subgroup was suppressed [2].Furthermore, the ERAD protein degradation system plays a critical role during seed development, when the biosynthesis of SSPs is intensely up-regulated.ERAD polyubiquitinates and then helps eliminate misfolded, damaged, and unfolded proteins by a process involving the ubiquitin ligase complex composed of OsHrd1 and OsOS-9 [35,37,46,47].The functional loss of OsHrd3 and the ubiquitin ligase complex is thought to play a role in forming deformed PB-Is [48], and reduced expression of OsDER1, a protein-interacting partner of OsHrd3, also led to the formation of cracked PB-Is.These results indicate that the quality of SSPs likely has a significant impact on PB-I structure.In this study, knockout of some of the genes in the Pro13a-I subgroup (pro13a.1/13a.2 in 1a-8-1 line and pro13a.1 in 2a-2-1) and the Pro13b-I/II subgroups (pro13b.1/13b.13 in 4b-9-7 and pro13b.3 in 8b-3-9) also led to notable structural changes in PB-I morphology.Compared to WT, mutant lines with reduced expression of Cys-rich Pro13a-I subgroup 13 kDa prolamins had small PB-I particles lacking lamellar spherical structure.In contrast, mutant lines with low expression of Cys-poor Pro13b-I/II subgroups had PB-I structures with a dark appearance, due to the changes in electron density (Figure 5A).

Functional Analysis of Genes Correlated with 13 kDa Prolamin Genes
In this study, subsets of 13 kDa prolamin genes were targeted/edited using CRISPR-Cas9, but compensatory changes in the expression of non-targeted genes were observed (Figure S7 and Figure 5).Specifically, in the immature seeds of the pro13a.1/13a.2-knockoutmutant line (1a-8-1), the transcription of genes in the Pro13a-II and Pro13b-I/II subgroups were up-regulated; in the pro13a.1-knockoutmutant line (2a-2-1), the transcription of genes in the Pro13a-II subgroup was up-regulated but those in the Pro13b-I/II subgroups were not changed.Moreover, in the immature seeds in which transcription of the genes in the Pro13b-I/II subgroups were suppressed (4b-9-7), a weak up-regulation of genes in the Pro13a-I subgroup was observed but not changed in the Pro13a-II subgroup.In the pro13b.3knockoutmutant line (8b-3-9), down-regulation of the genes in the Pro13b-I subgroup was accompanied by up-regulation of genes in the Pro13a-II, Pro13b-II, and Pro10 groups.Moreover, to gain insight into molecular events that occur in mutants suppressing the expression of specific prolamin subgroups in more detail, correlation network analysis identified six clusters of DEGs, each of which displayed unique patterns of expression in the four mutant plant lines.The DEGs played roles in RNA processing, RNA transcription, transport, stress response, protein synthesis, protein targeting, protein posttranslational modification, protein folding, and protein degradation.
During rice seed development, large amounts of SSPs are translated and matured in ER, in which the accumulation of unfolded proteins activates unfolded protein response (UPR), leading to enhanced expression of genes encoding chaperones such as BiPs and PDILs to help proper protein folding, proteins involved in protein secretion to target organelles, and proteins implicated in ER-associated degradation (ERAD) of misfolded proteins [49].Recent articles report that the suppression of SSP using RNAi and CRISPR-Cas9 systems can impact the expression of ER chaperones and co-chaperones such as BiPs, protein folding machinery (CNX, PDIs, and HSPs), and protein-homeostasis-related sHSPs, resulting in ER-stress and altered SSP composition [2,4,29,38].Moreover, the overexpression of ER stress marker genes (OsbZIP50, OsBiP1, OsBiP2, and OsBiP3) suppresses starch synthesis-related genes such as GPA3, FSE1, FLO7, and OsNF-YB1, resulting in grain chalkiness and lower amylose and SSP contents [16].The characterization of a rice opaque endosperm mutant, opaque3 (o3), with over-accumulation of 57 kDa proglutelins and lower SSP and starch contents indicated that O3 encodes OsbZIP60 that controls ER-homeostasis and grain quality by regulating genes involved in starch synthesis (GBSSI, AGPL2, SBEI, and ISA2), SSP (OsGluA2, Prol14, and Glb1), and chaperone (OsBiP1 and PDIL1-1) [50].
In this study, we showed that nine ERAD-related genes encoding heat shock proteins (three HSP70 and one HSP90), small heat shock proteins (HSP26, HSP20, HSP17.7, and HSP22a), and one co-chaperone (OsDjC35) were grouped in clusters 1 and 3 (Figure 8, Table S7).Cluster 1 includes genes that are positively correlated with prolamin genes in the Pro13a-II and Pro13b-I/II subgroups, including three OsHSP70 genes that are upregulated in pro13a.1-knockout(2a-2-1) and pro13b.3-knockout(8b-3-9) mutant lines.Cluster 3 contains genes that are negatively correlated with prolamin genes in the Pro13a-II and Pro13b-I/II subgroups but positively correlated with the Pro13a.1 gene in the Pro13a-I subgroup; these include OsHSP90, OsHSP26, OsHSP20, OsHSP17.7,OsHSP22a, and Os-DjC35.The genes were only up-regulated in the pro13b.1/13b.13-knockout(4b-9-7) mutant line.HSP70 functions as a molecular chaperone, facilitating the proper folding of nascent polypeptides, refolding misfolded proteins, and contributing to substrate degradation via the ubiquitin-proteasome system [51,52].HSP90, along with its co-chaperones, is involved in ERAD, a sophisticated cellular process responsible for identifying and degrading misfolded proteins through the ubiquitin-proteasome system [53][54][55].The small heat shock proteins (sHSPs) transfer denatured proteins to the HSP 70/90 chaperone system to perform continuous ATP-dependent refolding of these proteins, preventing protein misfolding and aggregation [56].DjC35, as a member of HSP40 co-chaperones (J-domain proteins), enhances HSP70 function by promoting ATPase activity, stabilizing substrate interactions, and preventing the aggregation of unfolded proteins [57][58][59].Up-regulation of the genes involved in protein processing pathways in the ER probably indicates that prolamin gene-knockout affects ER chaperones and SSP composition.These results suggest that the expression of ERAD-related genes and SSPs varied in the four mutant plant lines, causing a variable impact on grain quality.
Recently, it has been reported that TFs spontaneously regulate starch and SSP synthesis during seed development in rice [12], including OsbZIP58 [14,60], Nhd1 [61], OsNAC20 [15], OsNAC26 [15], FLO2 [62], and OsMADS1 [63].OsbZIP58 activates the expression of starch synthesis-related genes (OsAGPL3, GBSSI, SSIIa, SBE1, BEIIb, and ISA2) by directly binding their promoters [60] and interacts with rice prolamin box binding factor, activating the expression of SSP genes (GluA1~3, GluB1, GluD1, and 10/13/16 kDa prolamins) [14].Upregulation of N-mediated-heading-date (Nhd1), a MYB TF, resulted in reducing amylose and enhancing amylopectin through the down-regulation of GBSSI and up-regulation of SSII-3/SSIIa and FLO5, and suppressing the expression of GluA2, one of the SSP genes [61].Wang's group [15] reported that OsNAC20 and 26 up-regulate the expression of SSP genes (GluA1, GluB4/5, globulin, and 16 kDa prolamin) and starch synthesis-related genes (SSI, Pul, AGPS2b, AGPL2, and SBEI) together.The function-loss mutation of the FLOURY ENDOSPERM2 (FLO2) gene caused a decrease in grain size and starch content by reducing the expression of starch and SSP-related genes, while its overexpression resulted in an increase in amylose content and grain weight/size by enhancing the expression of starch (AGPL1, AGPL4, and GBSSI) and SSP genes (GluA1 and globulin) [62].Moreover, Oat-like rice (Olr), which was generated by a mutated OsMADS1 allele, showed a loose arrangement of starch granules, and reduced starch content but increased SSP content [63].The results indicate that the regulatory mechanisms of SSP and starch synthesis are coupled by TFs, and investigating novel TFs and understanding their regulatory networks are crucial factors in improving seed quality and yield.
In this study, we provide the involvement of 17 TFs belonging to the correlation network affected by 13 kDa prolamin gene expression, most of which are responsive to abiotic stress.These TFs belong to clusters 1, 2, 3, and 4 in the correlation network of genes.
In conclusion, based on the high homology among genes in a prolamin subgroup, we attempted to knockout two and seventeen genes in prolamin 13a-I and 13b-I/II subgroups using each sgRNA, which was designed at a specific region conserved among genes of each prolamin subgroup.The editing of up to two of the seventeen genes in the prolamin 13b-I/II subgroup suggests the editing limitations of one sgRNA and the need to use multiple sgRNAs.The sgRNA design in conserved regions represents a very useful strategy for editing genes with redundant functions and high homology.Furthermore, we investigated the relationship between 13 kDa prolamin suppression via the CRISPR-Cas9 tool and morphological changes in starch granules and protein bodies using RNA-Seq analysis, suggesting the involvement of ER stress-responsive TFs and (co-)chaperones.Although understanding the regulatory mechanisms is very complex and exploring their upstream genes has been limited, this study will provide important information for molecular breeding to improve grain quality.
The efficiency of the designed sgRNAs was tested in vitro using the Guide-it sgRNA In Vitro Transcription and Screening System.The designed sgRNAs were produced by in vitro transcription reactions according to the instructions of the Guide-it sgRNA In Vitro Transcription Kit (Cat.#632635, Takara Bio, Shiga, Japan).DNA templates were synthesized by PCR with specific primers (Table S8).Then, the synthesized DNA template, sgRNA, and recombinant Cas9 nuclease were combined as indicated in the Guide-it sgRNA screening kit (Cat.#632639, Takara) and incubated at 37 • C for 1 h.Reaction products were analyzed by agarose gel electrophoresis (Figure S2).

Cloning CRISPR-Cas9 Vector and Agrobacterium-Mediated Generation of 13 kDa Prolamin-Knockout Rice Plants
To insert the designed sgRNAs into the pRGE31 vector, two single-stranded 70 nt DNA oligonucleotides (ssDNA oligo) containing 20 nt sgRNA sequences (Table S8) were designed, synthesized, and purified with a final concentration of 0.2µM.The pRGE31 backbone vector was digested with Bsa I (Cat.#R3773, NEB, Ipswich, MA, USA) for 16 h at 37 • C. A reaction mixture consisting of 5 µL of ssDNA oligo (0.2 µM), 30 ng of restriction enzyme-linearized vector, and 10 µL of NEBuilder HiFi DNA Assembly Master Mix (Cat.#E2621, NEB) with a final volume of 20 µL was incubated for 1 h at 50 • C for ligation and then transformed into E. coli DH5-α cells.The recombinant clones were extracted and sequenced with the undigested pRGE31 as the negative control and the insertion donor as the positive control, as previously described [38].Furthermore, to construct binary vectors for Agrobacterium-mediated rice transformation, the sgRNA cassettes and pCAMBIA-Cas9 binary vector were double-digested with Hind III (Cat.#R3104, NEB,) and Bgl II (Cat.#R0144S, NEB) for 5 h at 37 • C, separated by gel electrophoresis, extracted, and then ligated with the ratio of vector and insert of 1:10 using T4 DNA ligase (Cat.#B0202A, NEB) at 16 • C for 13 h (Figure S3).Finally, the ligated binary vector was confirmed by colony PCR and DNA sequencing, using the undigested pCAMBIA-Cas9 as the negative control and the insertion cassettes as the positive control, with primers listed in Table S8.
Finally, the constructed binary vectors (Figure 1B) were transformed into the Agrobacterium tumefaciens strain (EHA105) and co-cultured with the calli induced from mature seeds of rice (Oryza sativa L. Japonica cv.'Ilmi').Then, the selection and regeneration of transformed calli were carried out as previously described [77].

Screening Transgenic Rice Plants
In order to detect mutations in transgenic rice plants, genomic DNA were extracted from the leaves (50 mg) of thirty-day-old rice plants as previously described [78].SpCas9 and HygR genes were amplified by PCR to screen for the presence/absence of the binary vector transformed using gene-specific primers (Table S9), and reaction products were analyzed by agarose gel electrophoresis (Figure S4).Furthermore, the targeted genes were amplified by PCR with specific primers for deep-sequencing (listed in Table S9), sequenced by the Mini-seq service (Bio-Core center, KAIST, Korea), and then analyzed using RGEN tools (http://www.rgenome.net/(accessed on 9 January 2024)).Mutation frequency was calculated based on the percentage of sequences with mutations (insertion or deletion of nucleotides) over the total number of sequences sequenced.We calculated that mutation frequency from deep sequencing data and compared the target regions in a sample of T 0 , T 1 , and T 2 plants.

RNA Isolation and qRT-PCR
Total RNA was isolated from three immature rice seeds at 2 WAF according to a previously published protocol [79].The cDNA was synthesized using 1 µg RNA and the QuantiTect Reverse Transcription Kit (Cat.#205311,Qiagen, Hilden, Germany) according to the manufacturer's instructions.qRT-PCR was performed as described in the QuantiTect SYBR Green PCR Kit (Cat.#204343,Qiagen, Hilden, Germany) in 20 µL reactions containing 5 ng cDNA, 0.5 µM forward primer, 0.5 µM reverse primer, and 2X Master Mix (10 µL).The Qiagen Rotor-gene Q cycler was used with the following settings: an initial denaturation at 95 • C for 10 min, 40 PCR cycles at 95 • C for 30 s, 60 • C for 30 s, and 72 • C for 30 s, a final extension at 72 • C for 10 min, and a melting curve analysis from 72 to 95 • C. The primers for qRT-PCR are listed in Table S8.Expression levels were normalized using the 2 −∆∆CT method with Ubiquitin (Os02g0161900) as the reference gene and WT as the reference sample [80] and represented as the log 2 value of the average fold change calculated from three biological replicates.

SSP Fractionation
Albumin/globulin, prolamin, and glutelin proteins were fractionated based on their solubility using different solvents as previously described [38,82].Protein content in each fraction was determined by bicinchoninic acid (BCA) assay (Pierce BCA Protein Assay Kit, Cat.#BCA1, Sigma-Aldrich) according to the manufacturer's instructions.

Measurement of Rice Agronomic Traits and Starch Content
Rice agronomic traits were quantified, including the average weight of 100 grains and the average length, width, and thickness of one grain.Ten de-husked mature seeds were ground in a mortar and pestle using liquid nitrogen.Starch and amylose content in 25 mg grain powder were measured using an Amylose/Amylopectin Assay Kit (Cat.#K-AMYL,Megazyme, Wicklow, Ireland) according to the manufacturer's instructions.

Microscopic Analysis
To observe the morphology of starch granules in rice endosperm, dry, mature seeds were transversely cut with a razor blade, mounted on SEM stubs, and coated with platinum particles.The mounted specimens were examined by a ZEISS Gemini 500 Scanning Electron Microscope (SEM) (Carl Zeiss AG, Jena, Germany) at 15 kV, as previous described [38].To observe PB-I and PB-II particles, immature rice seeds at 2 WAF were fixed in 50 mM cacodylate buffer (pH 7.2) containing 2% glutaraldehyde and 2% paraformaldehyde for 4 h at room temperature, washed three times in 50 mM cacodylate buffer (pH 7.2) for 30 min, fixed in 50 mM cacodylate buffer (pH 7.2) containing 1% OsO4 for 1 h, washed three times in 50 mM cacodylate buffer (pH 7.2) for 30 min, dehydrated in a gradient ethanol series (30,50,70,90,95, and 100%) at 30 min intervals, and then sequentially immersed in a mixture of LR White resin and 100% ethyl alcohol in ratios of 1:2, 1:1, 2:1, and 1:0 at 60 • C and 6 h intervals.Samples with LR White resin in a gelatin capsule were cured at 60 • C for 24 h, cut to 80-100 µm thickness using an ultra-microtome, mounted on carbon-coated nickel grids, double-stained with 4% uranyl acetate and 0.4% lead citrate, and observed with a JEM-2100F transmission electron microscope (TEM) (JEOL Ltd., Tokyo, Japan).

Statistical Analysis
A student's t-test (two samples assuming equal variances) was conducted using Excel (Microsoft office 365), which is used to compare means between two independent groups.For each comparison, p-values were calculated to determine statistical significance.Significant differences are indicated with asterisks (*), where * p < 0.1, ** p < 0.01, and *** p < 0.001.This analysis ensures that the reported differences in seed characteristics, starch composition, protein composition, and expression genes are statistically validated.

Figure 4 .
Figure 4. Starch granule appearance and starch content in WT and 13 kDa prolamin-knockout seeds.(A) Representative images of seed transverse sections showing the structure of starch granules.The red arrow indicates rod/filamentous granules.(B) Average weight of starch per grain.(C) Relative composition of amylose and amylopectin.(D) Transcripts of genes involved in starch metabolism were analyzed by qRT-PCR.Normalized expression of target genes was calculated using the 2 −ΔΔCT method and is shown as a log2 value.Error bars represent ± SD of three replicates.P-values were calculated using the Student's t-test (* p < 0.1, ** p < 0.01, and *** p < 0.001).

Figure 4 .
Figure 4. Starch granule appearance and starch content in WT and 13 kDa prolamin-knockout seeds.(A) Representative images of seed transverse sections showing the structure of starch granules.The red arrow indicates rod/filamentous granules.(B) Average weight of starch per grain.(C) Relative composition of amylose and amylopectin.(D) Transcripts of genes involved in starch metabolism were analyzed by qRT-PCR.Normalized expression of target genes was calculated using the 2 −∆∆CT method and is shown as a log 2 value.Error bars represent ± SD of three replicates.P-values were calculated using the Student's t-test (* p < 0.1, ** p < 0.01, and *** p < 0.001).

Figure 5 .
Figure 5. Structure of PB-I and composition of SSPs in Ilmi WT and 13 kDa prolamin-knockout lines.(A) TEM image of PB-I in developing endosperm at 2 WAF.(B) Average weight of SSP per grain.(C) Relative composition of SSPs (Albumin/globulin, Prolamin, and Glutelin).(D) qRT-PCR analysis of SSP genes in 13 kDa prolamin-knockout lines.Transcripts encoding different SSPs were analyzed in the immature seeds of WT and mutant lines.Normalized expression of the target genes was calculated using the 2 −ΔΔCT method and is represented as a log2 value.Error bars denote ± SD of three replicates.p-values were calculated using the Student's t-test (* p < 0.1, ** p < 0.01, and *** p < 0.001).

Figure 5 .
Figure 5. Structure of PB-I and composition of SSPs in Ilmi WT and 13 kDa prolamin-knockout lines.(A) TEM image of PB-I in developing endosperm at 2 WAF.(B) Average weight of SSP per grain.(C) Relative composition of SSPs (Albumin/globulin, Prolamin, and Glutelin).(D) qRT-PCR analysis of SSP genes in 13 kDa prolamin-knockout lines.Transcripts encoding different SSPs were analyzed in the immature seeds of WT and mutant lines.Normalized expression of the target genes was calculated using the 2 −∆∆CT method and is represented as a log 2 value.Error bars denote ± SD of three replicates.p-values were calculated using the Student's t-test (* p < 0.1, ** p < 0.01, and *** p < 0.001).

Figure 7 .
Figure 7. Correlation analysis among DEGs and 13 kDa prolamin genes.A Pearson correlation test was performed for DEGs and 13 kDa prolamin genes.The inclusion criteria was a correlation coefficient (r) ≥|0.7|.Transcripts were identified in nine functional groups: RNA processing, RNA transcription, transport, stress, protein synthesis, protein targeting, protein posttranslational modification, protein degradation, and protein folding.The data were visualized using Cytoscape software (ver.3.10.1).Circles indicate 13 kDa prolamin genes, and squares indicate correlated genes.The blue line indicates a negative correlation, and the red line indicates a positive correlation.

Figure 7 .
Figure 7. Correlation analysis among DEGs and 13 kDa prolamin genes.A Pearson correlation test was performed for DEGs and 13 kDa prolamin genes.The inclusion criteria was a correlation coefficient (r) ≥|0.7|.Transcripts were identified in nine functional groups: RNA processing, RNA transcription, transport, stress, protein synthesis, protein targeting, protein posttranslational modification, protein degradation, and protein folding.The data were visualized using Cytoscape software (ver.3.10.1).Circles indicate 13 kDa prolamin genes, and squares indicate correlated genes.The blue line indicates a negative correlation, and the red line indicates a positive correlation.

Figure 8 .
Figure 8. Expression of ER pathway genes in 13 kDa prolamin-knockout lines.(A) Relative expression of DEGs involved in ER-related pathways in 13 kDa prolamin-knockout lines.Normalized expression of each gene is represented by the Z-score.(B) qRT-PCR analysis of chaperone genes (BIP, PDIs, CNX, HSP90, HSP20, and HSP70) in 13 kDa prolamin-knockout lines.Normalized expression of the target genes was quantified using the 2 −ΔΔCT method, and the relative expression of each gene is represented as a log2 value.Error bars denote ± SD of three replicates.p-values were calculated using the Student's t-test (* p < 0.1, ** p < 0.01, and *** p < 0.001).

Figure 8 .
Figure 8. Expression of ER pathway genes in 13 kDa prolamin-knockout lines.(A) Relative expression of DEGs involved in ER-related pathways in 13 kDa prolamin-knockout lines.Normalized expression of each gene is represented by the Z-score.(B) qRT-PCR analysis of chaperone genes (BIP, PDIs, CNX, HSP90, HSP20, and HSP70) in 13 kDa prolamin-knockout lines.Normalized expression of the target genes was quantified using the 2 −∆∆CT method, and the relative expression of each gene is represented as a log 2 value.Error bars denote ± SD of three replicates.p-values were calculated using the Student's t-test (* p < 0.1, ** p < 0.01, and *** p < 0.001).

Figure 9 .
Figure 9. Expression of TF genes associated with 13 kDa prolamins.(A) The heatmap displays the expression profiles of the TF genes.This map shows altered expression of genes in 13 kDa prolaminknockout lines relative to WT. Normalized expression of each gene is represented by the Z-score.(B) Confirmation of the expression patterns of a few putative TF genes by qRT-PCR.Transcripts of TF genes were analyzed in the immature seeds.Normalized expression of the target genes was

Figure 9 .
Figure 9. Expression of TF genes associated with 13 kDa prolamins.(A) The heatmap displays the expression profiles of the TF genes.This map shows altered expression of genes in 13 kDa prolaminknockout lines relative to WT. Normalized expression of each gene is represented by the Z-score.(B) Confirmation of the expression patterns of a few putative TF genes by qRT-PCR.Transcripts of TF genes were analyzed in the immature seeds.Normalized expression of the target genes was calculated using the 2 −∆∆CT method, and is represented as a log 2 value.Error bars denote ± SD of three replicates.p-values were calculated by the Student's t-test (* p < 0.1 and ** p < 0.01).