22KD Zein Content Coordinates Transcriptional Activity during Starch Synthesis in Maize Endosperm

: Starch, the main form of stored energy in plants, plays an important role in maize ( Zea mays L.) kernel development. The Shrunken-2 ( Sh2 ) gene encodes the large subunit of the rate-limiting starch biosynthetic enzyme ADP-glucose pyrophosphorylase (AGPase). The sh2 mutant exhibits impaired AGPase activity, resulting in the partial or complete loss of starch synthesis. Here, we investigated the transcriptional regulatory framework of sh2 through transcriptome and co-expression network analysis using an F 2 population derived from the maize reference line B73 and sweet corn inbred line HZ508. We identiﬁed 5175 di ﬀ erentially expressed genes (DEGs), including 2878 upregulated and 2297 downregulated genes in sh2 mutant lines. DEGs are associated with various biological processes including nutrient reservoir activity, transferase activity, catalytic activity, water deprivation and glycogen metabolism. At the genetic level, 2465 DEGs, including 357 transcription factors, were involved in transcription. In addition, the maize ﬂoury and opaque mutant genes ﬂ1 , ndk2, o7 and o2 , which regulate the biosynthesis of 22KD zein, were co-expressed with the di ﬀ erential expressed transcription factor genes, thus suggesting that zein content might be a key regulator coordinating the expression of genes determining starch accumulation in maize endosperm.


Introduction
Starch regulation draws the attention of the scientific community, as this polymer is essential to meet energy needs in early plant development, human diet, feedstock and bioenergy. It accounts for up to 70% of the seed storage components in cereal endosperm and is thus a major contributor to yield as well as the main source of energy for seed germination and seedling growth. Starch production is a multistep mechanism starting with the conversion of Glc-1-P and ATP into ADP-Glc by ADP-glucose pyrophosphorylase (AGPase). This step is particularly important, as the inhibition of AGPase activity results in the partial or complete termination of starch biosynthesis, which greatly reduces yields [1,2]. AGPase is made of two subunits known as the large subunit (LSU) and small subunit (SSU). Genes encoding for AGPase subunits depend on the organ and tissue where they are localized. The fate of LSU and SSU are determined by Agpslzm and Agpllzm in the leaf, shrunken2 (Sh2) and brittle2 (Bt2) in the endosperm while Agplsmzm and Agplemzm are specifically expressed in the embryo [3]. Knocking down of Sh2 and Bt2 genes yields the starch-deficient mutant with a pericarp shrivel upon maturation [4]. As the mutant sh2 and bt2 genes behave in a complementary manner, the complete loss of function, which is lethal, requires a double homozygous mutant [4]. Heterologous expression of the gene encoding Escherichia coli glgC AGPase increased tuber yield by 35% in potato (Solanum tuberosum L.) [5], grain yield by 11% in rice (Oryza sativa L.) [6] and seed weight by 13% to 25% in maize [7]. Likewise, the inhibition of AGPase activity using site-specific mutagenesis system [8] and overexpression of Sh2 [9] also lead to an increase of maize seed weight by 11%-18%. Heterologous expression of a modified form of maize Sh2 (Sh2r6hs) increased yield by 38% in wheat (Triticum aestivum L.) [10], 23% in rice [11] and 64% in maize [12]. Although AGPase is crucial for starch biosynthesis and seed development in maize, the genetic and molecular mechanisms regulating this process need to be further clarified.
In this study, we aimed to get further insights into the genetic regulation of starch biosynthesis during maize kernel development.
We performed transcriptomic analysis using mutant sh2 kernels and wild-type kernels collected at 15 days after pollination (DAP), as well as co-expression network analysis using RNA-seq data from some well-known maize kernel mutants. Finally, we confirmed the reliability of our RNA-seq data via qRT-PCR using various differentially expressed genes (DEGs) related to starch or zein biosynthesis. Our analyses uncovered huge changes in the transcriptome of sh2 and point to four zein-related genes as major coordinators of transcriptional activity occurring during endosperm development in maize kernels.

Plant Materials and Sample Collection
Maize kernels harboring the starch deficiency (shrunken) phenotype in the HZ508 background (HZ508 sh2 ) were received from the agricultural academy of Enshi Tujia and Miao autonomous prefecture. This material was crossed with the maize inbred B73 to succinctly produce the F 1 and F 2 population used to carry this experiment. The F 1 plants were self-crossed; 15 days after pollination (DAP), samples were collected on developing F 2 ears. On each ear, 10 homozygous mutant kernels (sh2sh2) and homozygous wild-type kernels (Sh2Sh2) were snap-frozen in liquid nitrogen and stored at −80 • C until RNA extraction. The identity of each kernel was confirmed by genotype analysis, and each set of 10 kernels constitute a sample that was replicated three times.

RNA Extraction and Library Construction
Total RNA was extracted from pooled wild-type and sh2 mutant kernels obtained from the same F 2 ears derived from the cross B73 × HZ508 sh2 at 15 DAP (10 kernels per sample) using TRIzol reagent, followed by cleanup and DNase I treatment using a Qiagen RNeasy Mini Kit, according to the manufacturer's protocol. We used the Illumina NEBNext Ultra RNA Library Prep Kit ( cat.no E7530L New England Biolabs, MA, USA) to generate RNA-seq libraries following the manufacturer's recommendations. These libraries were sequenced on the Illumina NovaSeq platform by Novogene.

Identification of DEGs
The sequence reads were trimmed using Trimmomatic (version: 0.33) and mapped to the B73 reference genome (RefGen_v4.34) using HISAT2 (version: 2.1.0). StringTie (version: 1.3.3b) was employed to reconstruct the transcripts and to estimate gene expression levels [13]. HTSeq (version: 0.6.1) was used to count the number of reads per gene. When using the DEseq2 package, we applied the p-value adjusted for multiple testing (padj) as a threshold. Significant DEGs are those with a padj of differential expression above the threshold (p < 0.05, absolute value of log 2 (fold change) > 0.5). Gene ontology (GO) enrichment was performed for the significant DEGs using the agriGO Singular Enrichment Analysis tool with default parameters (http://bioinfo.cau.edu.cn/agriGO/) and REVIGO [14].

Confirmation of DEG
qRT-PCR was performed on 33 genes using SYBR Select Master Mix (Aidlab) on a CFX96 real-time PCR system (Bio-Rad, C1000 Touch). Relative expression levels were calculated using the comparative Ct method [15] and the maize actin gene (Zm00001d010159) as the internal control. All primers used are listed in Table S1.

Co-Expression Network Analysis of Well-Known Maize Kernel Mutant Transcriptomes
To construct a co-expression network, RNA-seq data from several well-known maize kernel mutants were obtained from the sequence read archive SRA database. Raw data filtering, mapping to the reference genome and gene expression level estimation were performed as described for RNA-seq data analysis. The weighted correlation network analysis (WGCNA) package in R was used to construct the gene co-expression network using the TPM matrix from StringTie. The analysis was conducted using an unsigned type of topological overlap matrix (TOM). For any gene pair (i and j), we applied the exponential adjacency function in the WGCNA algorithm to measure their relation index, the exponential weighted β square of the correlation coefficient (aij = power (Sij, β) = |Sij|β). Exponential weighted β is the power of the correlation coefficient. We selected β = 12 after the analysis (fit value R 2 to approximately 0.9). After determining the adjacency function parameter β, the correlation matrix S = [Sij] was switched into the adjacency matrix A = [aij] and converted into the topological overlap matrix Ω = [ωij]. ki or kj indicate the sum of one node's adjacency coefficients. Genes for which the sum in each sample was below 80 (<80) and the standard deviation (SD) was below 100 (<100) or above 6000 (>6000) were filtered out, yielding 14,528 genes for WGCNA co-expression network analysis.
The node is a gene (i or j). The hierarchical clustering tree was built using the dissimilarity coefficient, and the different branches represent the gene modules. We used the function of cutreeDynamic in WGCNA package to identify gene modules. The soft thresholding power was 12, the relatively minimum module size was 30 and the absolute value of correlation was greater than 0.25.

Isolation of Shrunken Mutant and DEGs
Using total RNA from the endosperm of sh2 mutant kernels and wild-type kernels at 15 DAP, we generated three mRNA libraries. After high-throughput sequencing and filtering of reads, we obtained 24.8-27.3 million clean reads (Table S2). We mapped these clean reads to the B73 reference genome (RefGen_v4.34). More than 91% of the clean reads were uniquely mapped to the reference genome for each sample (Table S3), and more than 94% of the mapped reads aligned to exons (Table S4). Comparing the structure of Sh2 gene in pooled wild-type and F 2 samples, we observe that the structure of Sh2 gene was altered in this mutant ( Figure S1). The wild-type Sh2 line has 15 transcribed exons while sh2 mutant line has only one exon. That difference also served as a basis for separating between wild-type and mutated shrunken kernels.
Identification of differentially expressed genes started by calculating the relative expression level of each gene using the unique mapped reads. Then, we normalized relative expression levels to the number of fragments per kilobase of the transcript sequence per million base pairs sequenced (FPKM) [16]. Principal component analysis based on these FPKM values showed that the three biological replicates of each sample were clustered into a single group ( Figure S2).

GO Enrichment Analysis of DEGs
Out of the 2297 downregulated genes in sh2, 74.44% could be functionally annotated using AgriGO. Most of the significantly enriched gene ontology (GO) terms belong to the molecular function category, including genes involved in nutrient reservoir activity, transferase activity, ATP binding, and catalytic activity. The biological process-related DEGs are involved in glutathione catabolic processes, organonitrogen compound catabolism, and sulfur compound catabolic processes. There were no significantly enriched GO terms in the cellular component category (Figure 1a, Table S5).
Out of the 2297 downregulated genes in sh2, 74.44% could be functionally annotated using AgriGO. Most of the significantly enriched gene ontology (GO) terms belong to the molecular function category, including genes involved in nutrient reservoir activity, transferase activity, ATP binding, and catalytic activity. The biological process-related DEGs are involved in glutathione catabolic processes, organonitrogen compound catabolism, and sulfur compound catabolic processes. There were no significantly enriched GO terms in the cellular component category ( Figure  1a, Table S5).
Out of the 2878 upregulated genes in sh2 mutant lines, 74.08% could be functionally annotated. DEGs in the biological process category include genes involved in the response to water deprivation, glycogen metabolism, nucleosome assembly and single-organism processes. Cellular-componentrelated genes are enriched in various GO categories, including nucleosome, cell junction, plasmodesmata and cell periphery. Molecular-function-related genes are involved in protein heterodimerization activity, oxidoreductase activity, cofactor binding, coenzyme binding, monooxygenase activity, acid phosphatase activity and DNA binding ( Figure 1b, Table S6).

Starch Biosynthesis and Zein-Related Genes Are Differentially Regulated in sh2 Kernels
We identified seven AGPase-encoding DEGs between mutant and wild-type kernels. Out of these DEGs, five were upregulated and two were downregulated genes along with Sh2 in mutant lines. We identified nine additional DEGs involved in starch biosynthesis, including one phosphoglucomutase (PGM) gene, three glycogen synthase genes, three glucan-branching enzyme (GBE) genes, and one starch debranching enzyme gene. Out of these, only the starch debranching enzyme gene was downregulated in sh2 mutant lines (Figure 2a). In addition, 27 zein genes were downregulated and five were upregulated, which may lead to the reduced zein levels in the mutant kernels ( Figure 2b). Many genes influencing zein accumulation and maize kernel development were differentially expressed between mutant and wild-type kernels, such as the classic semi-dominant Out of the 2878 upregulated genes in sh2 mutant lines, 74.08% could be functionally annotated. DEGs in the biological process category include genes involved in the response to water deprivation, glycogen metabolism, nucleosome assembly and single-organism processes. Cellular-component-related genes are enriched in various GO categories, including nucleosome, cell junction, plasmodesmata and cell periphery. Molecular-function-related genes are involved in protein heterodimerization activity, oxidoreductase activity, cofactor binding, coenzyme binding, monooxygenase activity, acid phosphatase activity and DNA binding ( Figure 1b, Table S6).

Starch Biosynthesis and Zein-Related Genes Are Differentially Regulated in sh2 Kernels
We identified seven AGPase-encoding DEGs between mutant and wild-type kernels. Out of these DEGs, five were upregulated and two were downregulated genes along with Sh2 in mutant lines. We identified nine additional DEGs involved in starch biosynthesis, including one phosphoglucomutase (PGM) gene, three glycogen synthase genes, three glucan-branching enzyme (GBE) genes, and one starch debranching enzyme gene. Out of these, only the starch debranching enzyme gene was downregulated in sh2 mutant lines (Figure 2a). In addition, 27 zein genes were downregulated and five were upregulated, which may lead to the reduced zein levels in the mutant kernels ( Figure 2b). Many genes influencing zein accumulation and maize kernel development were differentially expressed between mutant and wild-type kernels, such as the classic semi-dominant maize endosperm mutant genes fl1, fl2, fl3 and fl4 [18][19][20][21]; the classic recessive maize endosperm mutant genes o1 and o7 [21,22]; and the maize endosperm mutant genes pbf1, bZIP17, and bZIP22 [23,24] and PPDK1 (Figure 2c), all of which could affect storage protein gene expression [25].
To confirm the reliability of the RNA-seq data, we analyzed the expression levels of individual DEGs using qRT-PCR analysis in mutant and wild-type kernels and found a high association between the two datasets ( Figure 2d). The expression levels of many genes related to starch and zein biosynthesis were altered in the mutant, suggesting that the levels of starch and zein were altered due to the dysfunction of AGPase activity. All expression levels were normalized to that of the wild type. The data are from three biological replicates per sample and are presented as mean values ± SD. (Student's t-test, "·" represents p-value < 0.1, "*" represents p-value < 0.05, "**" represents p-value < 0.01).

Carbohydrate Catabolism and Electron Transport Chain Related Genes Are Differentially Regulated in sh2 Kernels.
As the weight of mature sh2 kernels was substantially reduced, we focused on genes in the carbohydrate catabolism pathway. We identified 103 DEGs in various carbohydrate catabolism pathways. For every enzyme in the Embden-Meyerhof-Parnas (EMP) pathway, at least one gene was differentially expressed, and most of these DEGs were upregulated in sh2 (Figure 3). Some of these enzymes are related to energy supply and maize kernel development, such as phosphofructokinase [26], and there is a competitive relationship between the EMP pathway and starch biosynthesis [27].
For the pentose phosphate pathway (PPP), 11 genes encoding five enzymes were differentially expressed, 10 of which were upregulated in sh2 mutant lines ( Figure 3). For every enzyme in the sucrose degradation pathway, at least one gene was differentially expressed, and most were upregulated in sh2 (Figure 3). Sus1 and Sus2 play important roles in sucrose metabolism and plant growth, but the functions of these genes are different, perhaps explaining why they had opposite All expression levels were normalized to that of the wild type. The data are from three biological replicates per sample and are presented as mean values ± SD. (Student's t-test, "•" represents p-value < 0.1, "*" represents p-value < 0.05, "**" represents p-value < 0.01).
To confirm the reliability of the RNA-seq data, we analyzed the expression levels of individual DEGs using qRT-PCR analysis in mutant and wild-type kernels and found a high association between the two datasets ( Figure 2d). The expression levels of many genes related to starch and zein biosynthesis were altered in the mutant, suggesting that the levels of starch and zein were altered due to the dysfunction of AGPase activity.

Carbohydrate Catabolism and Electron Transport Chain Related Genes Are Differentially Regulated in sh2 Kernels
As the weight of mature sh2 kernels was substantially reduced, we focused on genes in the carbohydrate catabolism pathway. We identified 103 DEGs in various carbohydrate catabolism pathways. For every enzyme in the Embden-Meyerhof-Parnas (EMP) pathway, at least one gene was differentially expressed, and most of these DEGs were upregulated in sh2 (Figure 3). Some of these enzymes are related to energy supply and maize kernel development, such as phosphofructokinase [26], and there is a competitive relationship between the EMP pathway and starch biosynthesis [27]. For the pentose phosphate pathway (PPP), 11 genes encoding five enzymes were differentially expressed, 10 of which were upregulated in sh2 mutant lines (Figure 3). For every enzyme in the sucrose degradation pathway, at least one gene was differentially expressed, and most were upregulated in sh2 (Figure 3). Sus1 and Sus2 play important roles in sucrose metabolism and plant growth, but the functions of these genes are different, perhaps explaining why they had opposite expression patterns in mutant and wild-type kernels [28]. No starch degradation genes were upregulated in mutant kernels (Figure 3).
Agronomy 2020, 10, x FOR PEER REVIEW 6 of 14 expression patterns in mutant and wild-type kernels [28]. No starch degradation genes were upregulated in mutant kernels (Figure 3). We identified 12 and 27 DEGs involved in the tricarboxylic acid cycle (TCA) and electron transport chain (ETC) (Figure 3), respectively, which help supply energy for plant growth and function in stress resistance, but only 17.9% (7 of 39) of these were downregulated.  [29] were used for pathway analysis described above from left to right: Pentose phosphate pathway (PPP), sucrose degradation, starch degradation, Electron Transport Chain (ETC), Tricarboxylic Acid Cycle (TCA) and the Embden-Meyerhof-Parnas pathway (EMP).

Hormone-Related DEGs
Hormones play important roles in regulating plant growth and kernel development. Ethylene negatively regulates those genes encoding enzymes related to starch biosynthesis [30]. We identified nine DEGs related to the ethylene metabolism pathway (Figure 4), and 34 ERFs (ethylene-responsive transcription factors) were differentially expressed between mutant and wild-type kernels, as revealed by RNA-seq (Table S7).
Abscisic acid (ABA) functions in abiotic stress resistance, and its levels of accumulation are highly correlated with starch accumulation [31]. Four genes related to ABA biosynthesis pathways and three genes related to ABA catabolic pathways were differentially expressed between mutant and wild-type kernels (Figure 4), including nced2, which codes for 9-cis-epoxycarotenoid dioxygenase (NCED) and catalyzes the conversion of 9-cis-epoxycarotenoids to xanthoxin, a key regulatory step in ABA biosynthesis [32]. Three genes (abh1, abh2 and abh5) encoding for ABA 8'hydroxylase, an enzyme involved in the ABA catabolic pathway, were also differentially expressed [33].
We identified 12 and 27 DEGs involved in the tricarboxylic acid cycle (TCA) and electron transport chain (ETC) (Figure 3), respectively, which help supply energy for plant growth and function in stress resistance, but only 17.9% (7 of 39) of these were downregulated.

Hormone-Related DEGs
Hormones play important roles in regulating plant growth and kernel development. Ethylene negatively regulates those genes encoding enzymes related to starch biosynthesis [30]. We identified nine DEGs related to the ethylene metabolism pathway (Figure 4), and 34 ERFs (ethylene-responsive transcription factors) were differentially expressed between mutant and wild-type kernels, as revealed by RNA-seq (Table S7).
Abscisic acid (ABA) functions in abiotic stress resistance, and its levels of accumulation are highly correlated with starch accumulation [31]. Four genes related to ABA biosynthesis pathways and three genes related to ABA catabolic pathways were differentially expressed between mutant and wild-type kernels (Figure 4), including nced2, which codes for 9-cis-epoxycarotenoid dioxygenase (NCED) and catalyzes the conversion of 9-cis-epoxycarotenoids to xanthoxin, a key regulatory step in ABA biosynthesis [32]. Three genes (abh1, abh2 and abh5) encoding for ABA 8'-hydroxylase, an enzyme involved in the ABA catabolic pathway, were also differentially expressed [33].
Agronomy 2020, 10, x FOR PEER REVIEW 7 of 14 increased CK (cytokinin) activity leads to higher grain number and improved grain yield [38]. We identified five DEGs related to CK degradation and biosynthesis pathways (Figure 4), including cko1, the orthologue of which regulates grain shape in rice [39].

Differentially Expressed TF Genes between the Wild Type and sh2
To investigate potential regulators involved in starch biosynthesis and kernel development, we focused on differentially expressed TF genes based on annotations in PlantTFDB v4.0 [17]. We identified 160 (6.97%) downregulated genes and 197 (6.85%) upregulated genes that were annotated as TF genes (Table S7).
According to PlantTFDB these TF genes are classified into 45 TF families ( Figure 5). Transcription factors belonging to the ethylene-responsive family (ERF), NAC factors, MYB, GRAS, G2-like, bZIP and bHLH are the most abundant, with DEG amount ranging from 34 to 14 (Table S7).

Differentially Expressed TF Genes between the Wild Type and sh2
To investigate potential regulators involved in starch biosynthesis and kernel development, we focused on differentially expressed TF genes based on annotations in PlantTFDB v4.0 [17]. We identified 160 (6.97%) downregulated genes and 197 (6.85%) upregulated genes that were annotated as TF genes (Table S7).
According to PlantTFDB these TF genes are classified into 45 TF families ( Figure 5). Transcription factors belonging to the ethylene-responsive family (ERF), NAC factors, MYB, GRAS, G2-like, bZIP and bHLH are the most abundant, with DEG amount ranging from 34 to 14 (Table S7).

Transcriptional Activity is a Key Regulator of Starch Synthesis at 15 DAP
We found numerous DEGs which related to many biological processes in our RNA-seq data but did not reveal key regulators of maize kernel development and starch biosynthesis that were not altered in the sh2 mutant. To resolve this problem, we performed weighted correlation network analysis (WGCNA) [41] to build a co-expression network using published RNA-seq data from various maize kernel mutants, including o2, pbf, mads47, o11, fl3, fl4, dek2, dek10, dek35, dek37, sh2 and the o2-pbf double mutant (Table S8). We used HISAT2 and StringTie to obtain the transcripts per kilobase of exon model per million mapped read (TPM) [13] for every gene in each sample. Genes expressed in each sample were filtered out with moderate SD in different samples, leading to the identification of 20 modules containing these genes (Figure 6a). Finally, 6623 nodes and 19 modules were obtained by applying link weights above 0.1 with each module containing >3 and <2465 nodes (Table 1).
To investigate potential regulators involved in starch biosynthesis and kernel development, we focused on differentially expressed TF genes based on annotations in PlantTFDB v4.0 [17]. We identified 160 (6.97%) downregulated genes and 197 (6.85%) upregulated genes that were annotated as TF genes (Table S7).
According to PlantTFDB these TF genes are classified into 45 TF families ( Figure 5). Transcription factors belonging to the ethylene-responsive family (ERF), NAC factors, MYB, GRAS, G2-like, bZIP and bHLH are the most abundant, with DEG amount ranging from 34 to 14 (Table S7).  The graph shows the number of upregulated genes (log 2 (fold change) > 1 or log 2 (fold change) < −1) in each family identified in the mutant (red) and wild-type (blue) and the number of upregulated genes (log 2 (fold change) > 0.5 or log 2 (fold change) < −0.5) in each family identified in the mutant (black) and wild-type (gray).

Co-Expression Network Analysis
To further understand the transcriptional mechanism regulating starch biosynthesis in early endosperm development (15 DAP) in maize, we selected four maize endosperm mutants belonging to the "orange" module ( Figure 6). RNA-seq data from these mutants (fl1, o7, nkd2 and o2) were downloaded from the SRA database. The imprinted floury1 mutant gene fl1 presented greater coexpression with genes related to transcriptional activity (Figure 7). It is followed by nkd2 and o7, while  The "orange" module had a high module membership (MM) with sh2, with a correlation coefficient of 0.97 and a p-value of 8 × 10 −8 (Figure 6b, Figure S4). We subjected the genes in that module to GO enrichment analysis and identified only one significant GO term: TF activity (sequence-specific DNA binding; GO:0003700, p-value = 3.5 × 10 −9 ). This implies that many genes related to TF activity are included in this module. Out of the 2465 genes found in this module, 249 are TF genes, which include 70 upregulated and 15 downregulated DEGs in sh2, as determined by annotation using PlantTFDB ( Figure S4). Most of these genes encode members of the ERF, bZIP, bHLH, GRAS, Trihelix, MYB, GATA, C2H2 and NAC TF families, and nine TF families include more than 10 members. The bZIP TF family, with 24 members, includes o2, a zein mutant known to network regulate maize endosperm development by directly regulating zein and other TF genes [24,42]. We identified 22 bHLH TF genes, including o11, encoding a key regulator of endosperm development and nutrient metabolism [43]. Twenty MYB and MYB-related TF genes are included, such as MYB14, which regulates the starch biosynthesis pathway in maize by regulating ZmBT1 promoter activity [44]. Other genes in the "orange" module, such as Zm00001d021537 and Zm00001d038288, are targets of o11 in the MYB TF family [43]. Two targets of o11, Zm00001d023294 and Zm00001d049860, encode NAC TFs [43]. Although transcriptional regulation is essential in starch synthesis, its global network is yet to be identified.

Co-Expression Network Analysis
To further understand the transcriptional mechanism regulating starch biosynthesis in early endosperm development (15 DAP) in maize, we selected four maize endosperm mutants belonging to the "orange" module ( Figure 6). RNA-seq data from these mutants (fl1, o7, nkd2 and o2) were downloaded from the SRA database. The imprinted floury1 mutant gene fl1 presented greater co-expression with genes related to transcriptional activity ( Figure 7). It is followed by nkd2 and o7, while o2 is only associated to five genes.

Discussion
Malfunction of AGPase resulting from loss of Sh2 function causes disruption of starch biosynthesis pathway, resulting in an increase of soluble sugar in kernels that are not fully formed and causing mature kernels to have reduced weight [45,46]. Insight into the transcriptional regulatory framework of sh2 mutant in immature maize kernels can improve our understanding of genetic regulation of starch accumulation during endosperm development. To be enzymatically active, AGPase requires the complementary activity of its two subunits encoded by Bt2 and Sh2 genes in maize endosperm [4,47,48]. The effect of mutational changes in these genes leads to various phenotypes, which range from similarity to wild-type to lethality depending on the combination of the Bt2 and Sh2 alleles [49]. We identified a sh2 allele in which only one exon was transcribed, while the wild type had 15 exons transcribed in the B73 background ( Figure S1). RNA-seq analysis performed on this mutant at 15 DAP revealed 5175 DEGs identified from sh2 mutant line kernels and wild-type kernels. DEGs related to endosperm include 11 genes for kernel development, 16 genes for starch biosynthesis and 32 zein-encoding genes. DEGs related to starch, including bt2, were upregulated ( Figure 2). This alteration of the shrunken2 gene also affects the expression of proteinrelated genes. Among the 32 zein-encoding genes identified in our experimental conditions, 25 genes were downregulated (Figure 2b). These data support that starch and zein synthesis are strongly associated biological processes which likely depend on transcript content of AGPase-related genes [48,50,51].
Starch synthesis occurs at 10-35 DAP, when the carbohydrate, storage protein and amino acid metabolism enzymes have accumulated [52]. Transcriptional regulation of expression of proteins was

Discussion
Malfunction of AGPase resulting from loss of Sh2 function causes disruption of starch biosynthesis pathway, resulting in an increase of soluble sugar in kernels that are not fully formed and causing mature kernels to have reduced weight [45,46]. Insight into the transcriptional regulatory framework of sh2 mutant in immature maize kernels can improve our understanding of genetic regulation of starch accumulation during endosperm development. To be enzymatically active, AGPase requires the complementary activity of its two subunits encoded by Bt2 and Sh2 genes in maize endosperm [4,47,48]. The effect of mutational changes in these genes leads to various phenotypes, which range from similarity to wild-type to lethality depending on the combination of the Bt2 and Sh2 alleles [49]. We identified a sh2 allele in which only one exon was transcribed, while the wild type had 15 exons transcribed in the B73 background ( Figure S1). RNA-seq analysis performed on this mutant at 15 DAP revealed 5175 DEGs identified from sh2 mutant line kernels and wild-type kernels. DEGs related to endosperm include 11 genes for kernel development, 16 genes for starch biosynthesis and 32 zein-encoding genes. DEGs related to starch, including bt2, were upregulated ( Figure 2). This alteration of the shrunken2 gene also affects the expression of protein-related genes. Among the 32 zein-encoding genes identified in our experimental conditions, 25 genes were downregulated (Figure 2b). These data support that starch and zein synthesis are strongly associated biological processes which likely depend on transcript content of AGPase-related genes [48,50,51].
Starch synthesis occurs at 10-35 DAP, when the carbohydrate, storage protein and amino acid metabolism enzymes have accumulated [52]. Transcriptional regulation of expression of proteins was initially suggested to be the regulatory mechanism determining this coordinated expression [52]. Our co-expression network analysis shows that the module with the highest module membership had 2465 genes involved in transcriptional activity. A total of 357 TF genes, or 6.89% of the genome-wide DEGs, were identified, thus supporting the importance of transcriptional control in the regulation of starch biosynthesis [43,50,53,54]. In addition, four known mutants fl1, o7, nkd2 and o2 are co-expressed (Figure 7). The mutant of floury gene 1 (fl1) shows the greatest amount of co-expressed genes. Previous work reports that fl1 encodes a zein protein body membrane accumulating 22-kD α-zeins in the γ-zein-rich periphery and center of the core rather than their normal discrete location in a ring at outer edge of the core [18]. This mutation results in the formation of starchy endosperm, while the wild-type allele Fl1 causes the vitreous endosperm. This is consistent with a putative role as a coordinator of transcriptional activity in starch synthesis at 15 DAP. The maize opaque 7 mutant gene (o7) appear to be the second TF gene highly co-expressed with maize endosperm genes at this developmental stage. In addition to decreasing starch α-zein contents, o7 influences amino acid biosynthesis by affecting a-ketoglutaric acid and oxaloacetic acid both in maize and rice [22,55]. Although less important than fl1, the naked endosperm 2 (nkd2) and the opaque 7 (o7) mutant TF genes are also co-expressed with a wide number of genes in the "orange" module. ndk2 TF is a member of the indeterminant domain (IDD) TF family, which is implicated in many biological functions, including the regulation of carbohydrate metabolism, gravitropism, seed germination, lateral organ morphogenesis, gravitropism, hormone signaling or DNA binding factor in Arabidopsis [56]. In maize, characterization of ndk1 ndk2 double mutant into B73 background revealed that ndk1 and ndk2 influence 6.7% of transcriptome in starchy endosperm (at 15 DAP) [57], and NDK2 activates transcription from the 22KD zein promoter. The fourth connector is o2, a gene that was previously identified in the regulation of starch and protein synthesis in maize endosperm at 16 DAP [23]. However, the genetic network regulated by o2 is reduced to five genes only in our experimental conditions. The difference is likely due to the pathway controlled by these two genes. While sh2 controls starch synthesis by regulating the activity of AGPase, the transcription factor o2 induces the downregulation of genes involved in the pentose phosphate pathway. Altogether, these data confirm that coordinated expression of starch is controlled by protein content. Our co-expression data highlight the specific role of 22KD zein content in sh2-deficient mutant.

Conclusions
Control of starch production and storage in the endosperm determines seed weight and kernel properties. In maize endosperm, this mechanism is regulated by several enzymes, among which AGPase plays a determinant role. Transcriptome profiling supports that the modification of sh2 transcript content influences not only AGPase activity, but also the expression of genes controlling starch and protein synthesis during endosperm development. This regulation is carried out through a transcriptional network capped by four TF, namely the floury mutant gene 1 (fl1), the opaque 7 gene, the naked endosperm 2 (nkd2) gene and opaque 2 TF gene. Previous functional characterization of these genes demonstrated that they activate the transcription of genes involved in the storage of 22KD zein, thus supporting that accumulation of zein is a determinant for starch accumulation.