Comparative transcriptomics reveals the difference in early endosperm development between maize with different amylose contents

Plant materials and sequencing

The common maize (Zea mays. L) hybrid Shandan 609 (SD609) and the high-amylose maize hybrid high-amylose 68 (HS68) were grown under field conditions at a density of approximately 67,500 plants ha⁻¹ in the summer of 2013 in Yangling, Shaanxi Province, China. The ears were bagged prior to silk emergence and then subjected to manual self-pollination. Three biological replicates of the ears were harvested at 5, 10, 15 and 20 DAP. To ensure the integrity and specificity of the endosperm tissues and the absence of contamination by any other tissues, all slightly damaged endosperm tissues were discarded, and only complete tissues were frozen immediately in liquid nitrogen and stored at −80 °C prior to RNA extraction.

The total RNA of three biological replicates of the endosperm was extracted separately from each sample using the TRIzol-A+ reagent (Invitrogen) following the manufacturer’s instructions. RNA sequencing libraries were constructed according to the standard Illumina protocol and were sequenced using the Illumina HiSeq™ 2000 platform to generate 2 ×100-nucleotide paired-end reads (Qu et al., 2016).

The sequence data obtained in this study have been deposited in the National Center for Biotechnology Information Sequence Read Archive (http://www.ncbi.nlm.nih.gov/sra) under the accession numbers SRP065059 and SRP149609.

Read mapping and expression analysis

RNA-Seq reads from each sample were aligned to the maize B73 reference genome (ZmB73_RefGen_v2; http://www.maizesequence.org/) using the Tophat2 (v2.09; http://ccb.jhu.edu/software/tophat/index.shtml) programme with Bowtie2 read alignment software (Trapnell, Pachter & Salzberg, 2009; Langmead et al., 2009). The aligned reads obtained after mapping using Tophat2 were subjected to transcript assembly using Cufflinks software (version 2.1.0; http://cole-trapnell-lab.github.io/cufflinks). The expression level was calculated in fragments per kilobase of transcript per million fragments mapped (FPKM) as described previously (Qu et al., 2016). The differentially expressed genes (DEGs) were identified through pairwise comparison using EBSeq software, and the P-values were adjusted using the Benjamini–Hochberg procedure to determine the false discovery rate (FDR) (Benjamini & Hochberg, 1995; Robinson, McCarthy & Smyth, 2010; Leng et al., 2013). Only the genes with FPKM results that met the criteria of FDR <0.01 and fold change  ≥2 or ≤0.5 between the two conditions were deemed to be significantly differentially expressed.

Gene coexpression network analysis

A gene coexpression network analysis was performed using the WGCNA R package (version1.63; https://horvath.genetics.ucla.edu/html/CoexpressionNetwork/Rpackages/WGCNA/) (Langfelder & Horvath, 2008). We selected all DEGs for the network analysis based on their variance across different samples. We subsequently calculated the Pearson correlation coefficients for each gene-gene comparison and then constructed an adjacency matrix of the connection strengths. The power β was optimized to 6 to adjust the scale-free property of the coexpression network and the sparsity of the connections between genes. Modules were defined as sets of genes showing high topological overlap. An adjacency matrix was used to calculate the topographical overlap matrix (TOM), and a TOM-based dissimilarity measure (1–TOM) was used for hierarchical clustering. The cutreeDynamic function was used to define the gene modules corresponding to the branches of the resulting dendrogram. Highly similar clusters were merged in the network using the mergeCloseModules function with a cutHeight value of 0.25. Each module was assigned a unique number. The NetworkAnalyzer plugin available in Cytoscape was used to calculate the relevant network parameters, such as the degree of connection (Assenov et al., 2008), which measures the number of incoming and outgoing edges of a node. A higher degree of connection indicates that the node has a higher number of edges, and thus, this property can be used to identify whether a gene is a major hub of the network and thus potentially a major regulatory gene in a pathway. The top 5% of genes in each module were identified as hub genes. The modules were visualized using Cytoscape (version 3.6.1; http://www.cytoscape.org/download.php), and the layout was set based on an edge-weighted spring-embedded layout (Shannon et al., 2003).

Statistical analyses

Functional enrichment analyses were performed using the KEGG database (https://www.genome.jp/kegg/) and agriGO (http://bioinfo.cau.edu.cn/agriGO/analysis.php; version 2.0) (Du et al., 2010). A gene ontology (GO) term was considered significantly enriched if the adjusted p-value was lower than 0.05 compared with the genome-wide background. The K-mean clustering method, Euclidean distance and Calinski-Harabasz index were used to identify the gene expression patterns (Calinski & Harabasz, 1974). Promoter sequences of the core genes in starch synthesis were obtained using the Ensembl Plants database (http://plants.ensembl.org/index.html). The 2,000 bps upstream of the start codons were analysed using the online software PLANTCARE (Lescot et al., 2002). If the distance between two genes was less than 2,000 bps, the sequence between the start codon of the first gene and the termination codon of the second gene were searched for cis-regulatory elements using the PlantCARE database. Additionally, all statistical analyses were performed using the R language (R Core Team, 2019).

Real-time quantitative PCR

To verify the RNA-Seq results, 11 DEGs closely related to endosperm development and starch metabolism were selected for quantitative real-time PCR (qRT-PCR) using SYBR Green I (Bio-Rad) and a CFX96 real-time PCR detection system (Bio-Rad, Hercules, USA). The gene-specific primers used in this assay were designed using Primer Premier 5.0 software (Premier, CAN) and are listed in Table S1, and the specificity of the primers was further verified using the non-redundant (Nr) database and Primer-BLAST (GenBank, NCBI). Each sample from each time point was divided into three biological replicates. Quantification was performed using the comparative CT method with actin (gene ID: GRMZM2G082484) as the internal control. The amplification conditions were 95 °C for 15 min and 40 cycles of 95 °C for 10 s, 60 °C for 20 s and 72 °C for 30 s, and these steps were followed by a melting curve programme (95 °C for 10 s followed by an increase from 65 °C to 95 °C in 5-s increments of 0.5 °C) to confirm the specificity of the PCR amplification. The relative expression levels were calculated as previously described (Livak & Schmittgen, 2001).

Isolation of starch and determination of the amylose content

To eliminate interference from other ingredient, fat was removed using diethyl ether prior to the isolation of starch, and 200 g of corn flour was then steeped in 1,000 mL of distilled water. The mixture was evenly stirred, and the homogenate was passed three times through a high-precision stainless-steel sieve (150 mu/0.106 mm) to improve the purity of the isolated starch. The slurries were then poured into 50-mL centrifuge tubes and centrifuged at 2,000 rpm for 20 min. This procedure was repeated three times. The supernatant in the centrifuge tube was then discarded, and the sediment was suspended in 40 mL of 0.5% (w/v) sodium hydroxide solution for at least 4 h. The supernatant was discarded, and 40 mL of new 0.5% (w/v) sodium hydroxide solution was added. After the mixture was evenly stirred, the centrifuge tube was centrifuged at 4,000 rpm for 10 min. These steps were repeated three times. Subsequently, 30 mL of absolute ethanol was added to the centrifuge tube, and the tube was centrifuged at 4,000 rpm for 10 min. The final sediment was then oven-dried (40 °C) or air-dried and ground using a mortar and pestle. The amylose content of the B73 starch was determined from the size-exclusion chromatography weight distributions of the debranched starch and was calculated as the ratio of the area under the curve of the amylose branches (defined to have a degree of polymerization P 100) to the area under the curve of the entire distribution (including both the amylopectin and amylopectin branches) (Fitzgerald et al., 2009). The amylose content of the sequenced samples (SD609 and HS68) was determined using a two-wavelength iodine binding procedure, which was previously described in detail (Zhu et al., 2008). Two amylose determination methods were used to ensure high-accuracy measurements.

Overview of the RNA-Seq datasets

To obtain global transcriptome maps of the early endosperms of maize with different amylose contents, we performed an RNA-Seq analysis of the common maize SD609 (average amylose content of 27.70%) and the high-amylose maize HS68 (average amylose content of 50.76%) (Fig. 1A). Specifically, the Illumina HiSeq2000 platform was used to perform paired-end RNA-Seq analyses of the endosperm tissues from the two maize plants at 5, 10, 15 and 20 DAP. We generated approximately 0.4 billion high-quality reads, and 95.34% of the mapped reads could be uniquely mapped to the B73 reference genome (Table S2). Additionally, the visualizations of the reads mapped to the maize genome sequences varied across the four developmental stages, which indicated that the transcriptome is dynamic during the early stage of maize endosperm development and that the RNA-Seq data were reliable for further analysis (Fig. S1).

Dynamic gene expression patterns and cellular function

To gain insight into the temporal and spatial gene expression patterns during the early stage of maize endosperm development, the gene expression levels at each stage of endosperm development in SD609 and HS68 were classified into five categories based on their FPKM values (Fig. 1B, Table S3). The highest level had the lowest number of expressed genes, i.e., approximately 2,000, whereas the lowest level had more expressed genes than the other levels, i.e., approximately 8,000 (Fig. 1). In particular, the difference in the overall number of expressed genes between SD609 and HS68 was only 326, and 22,494 and 22,358 genes were expressed in the endosperm of B73 at 10 and 20 DAP, respectively. Approximately 6,000 extended (expressed) genes were detected in SD609 and HS68 at 10 (29,205 and 28,421, respectively) and 20 DAP (28,723 and 28,144, respectively) (Chen et al., 2014).

To further understand the divergence in expressed genes between SD609 and HS68 and reduce the influence of transcription noise, only genes with FPKM values ≥ 1 were used in the estimation of the number of expressed genes. In total, 21,986 expressed genes were found in the SD609 and HS68 endosperms at all four early stages of development (Fig. 1C), and a GO enrichment analysis suggested that these common genes were mainly involved in cellular carbohydrate biosynthesis, cellular protein metabolism, cellular nitrogen compound biosynthesis, and other essential processes of endosperm development (Fig. S2). A total of 390 genes were expressed only during specific stages of endosperm development in SD609 and HS68, and most of these genes were found to be related to the response to stress and the regulation of gene expression and cell development (Table S4). In particular, sugary enhancer 1 (se1; GRMZM5G842214), which is upstream of the starch metabolism pathway, was only detected at 5 DAP in HS68. However, previous genetic studies have uncovered that se1 is a recessive modifier of su1 that leads to a higher sugar content and an increased kernel moisture for longer post-harvest periods (Ferguson, Rhodes & Dickinson, 1978; Carey, Rhodes & Dickinson, 1982). In addition, increased maltose concentrations have also been detected in mature se1 kernels (Ferguson, Dickinson & Rhodes, 1979; Tadmor et al., 1995). These characteristics indicate that se1 might be a regulator that affects the ratio of amylose to amylopectin in the early stage of endosperm development in HS68.

Overview of gene expression patterns

To understand the gene expression patterns in maize with different amylose contents, pairwise comparisons of all developmental stages of SD609 and HS68 were performed. After genes with read counts lower than 10 at each stage were filtered out, a total of 15,922 DEGs were identified based on the criteria log₂fold change (log₂FC) ≥1 or ≤-1 and FDR ≤ 0.05. Of these, 1,193 DEGs were obtained from the comparison between SD609 and HS68 at the various developmental stages, 7,813 DEGs were consistently expressed in SD609 and HS68 at the different developmental stages, and 4,914 and 2,002 DEGs were DEGs that were only detected in SD609 and HS68, respectively, and are thus considered material-specific DEGs (Fig. S3, Table S5). A functional enrichment analysis indicated that 1,193 DEGs between SD609 and HS68 were mainly involved in regulation of the intracellular signalling cascade, response to hormone stimulus, system development and fatty acid biosynthetic process (Fig. S4A). In addition, the material-specific DEGs were also found to regulate the same biological processes, such as the intracellular signalling cascade, fatty acid biosynthetic process, dephosphorylation, and response to external stimuli (Fig. S4B, Table S6). Notably, only the material-specific DEGs detected in HS68 were found to be involved in the response to water deprivation, and these DEGs in HS68, which can potentially be attributed to the higher amylose content of this maize variety compared with that of common maize, likely induced the loose binding of starch grains and ultimately resulted in rapid water dispersal, as observed in these kernels at different stages of development but not in those of common maize (Kibar, Gönenç & Us, 2010; Ai & Jane, 2018) (Figs. S4C–S4D). In contrast, the material-specific DEGs of SD609 were enriched in diverse biological processes, such as cell communication, the defence response, nucleosome assembly and the response to oxidative stress (Fig. S4D). These results indicated the existence of obvious metabolic differences between common and high-amylose maize at the early stage of kernel development.

A temporal expression pattern analysis of the shared DEGs between SD609 and HS68 was conducted based on k-means clustering and the CH index, and 16 clusters were identified. The results showed that although the same genes were differentially expressed in the two maize materials, their expression patterns were not synchronized and showed a single or multiple transition points at the different development stages. For instance, clusters 6, 11 and 15 exhibited significant differences in expression transition (Fig. 2). In addition, some clusters, such as clusters 7, 8 and 12, showed synchronous fluctuations in their expression patterns but differences in their expression levels (Fig. 2). The qRT-PCR patterns for all candidate genes were almost equivalent to the RNA-Seq patterns (Fig. 3, Table S7). A stringent GO term enrichment analysis (P < 0.0001, FDR < 0.01) of the shared DEGs revealed that the gene sets of most clusters were involved in the regulation of system development, reproductive process, and multicellular organismal development (Fig. S5A, Table S6). Additionally, some clusters were highly correlated with specific biological processes, such as negative regulation of transcription (cluster 5), glucan metabolic process (cluster 6), regulation of catalytic activity (cluster 11) and nucleosome metabolic process (cluster 16) (Fig. S5A).

We further classified the material-specific DEGs and found no unique expression transition pattern for a specific material. A thorough analysis of the functional enrichment of the material-specific DEGs revealed that some unique expression patterns were closely related to key biological processes in endosperm development. For instance, cluster 9 of SD609 and cluster 3 of HS68, both of which are related to cell differentiation, showed the same expression patterns, and clusters 3 and 8 of SD609 and clusters 1 and 10 of HS68, all of which are involved in the regulation of organ development, exhibited synchronous expression transitions. In addition, cluster 4 of SD609 and cluster 2 of HS68, which are related to negative regulation of biological process, also showed the same expression pattern (Fig. S5–S6, Table S6). Interestingly, we found opposite expression patterns for clusters that are relevant to the same biological process; for example, clusters 3 and 4 of SD609, which are related to negative regulation of biological process, exhibited the opposite expression patterns (Figs. S5–S6). A similar finding was also found for some DEGs between the two different materials; for example, cluster 9 of SD609 and cluster 3 of HS68, which are involved in the regulation of growth, presented opposite expression transitions (Figs. S5–S6). These results indicated that our gene expression pattern sets provide a sufficient dynamic atlas for the identification of genes that function in key biological processes during early maize endosperm development and provide valuable insights for annotating the function of putative genes.

Identification of weighted gene coexpression network analysis (WGCNA) modules associated with starch synthesis-related enzymes

Starch metabolism is a complex physiological process due to its close connection with sucrose metabolism and its regulation by dozens of enzymes (Fig. 4). Additionally, some key enzymes play crucial roles in regulating the synthesis of amylose and amylopectin and might also be crucial factors responsible for the differences in the amylose content between SD609 and HS68. To obtain an in-depth understanding of the potential regulatory mechanisms and core gene networks responsible for the differences in the amylose content between SD609 and HS68, we performed a WGCNA using the 15,922 above-described DEGs and identified 31 WGCNA modules (Fig. 5). We searched the core genes involved in starch synthesis in all the modules based on a previous study (Qu et al., 2018). In total, 19 of the 28 core genes associated with starch synthesis were differentially expressed and were distributed among 10 modules, namely, Plum (AGPLS1), Plum_1 (SSIIa), Steel_blue (AGPLS4), Blue (SBEIIb and SSIV), Dark_magenta (su1), Dark_orange (SBEI, SBEIII and SSIIIb-b), Dark_violet (SSIIb), Grey_60 (AGPLS2, AGPSS3, PUL and SSI), Maroon (SSIIIb-a) and Orange-red_3 (AGPSS1, GBSSI, GBSSII and SSIIIa) (Fig. 6A).

A GO enrichment analysis revealed that genes in each of the 10 modules were not only significantly enriched in specific biological processes but also involved in the regulation of common pathways. For example, the Orange-red_3 module includes the most critical enzyme, GBSS, which regulates amylose synthesis, and other genes belonging to this module were specifically and significantly enriched in the lipid biosynthetic process, fatty acid biosynthetic process, and cell proliferation (Fig. 6B, Table S8). Therefore, the Orange-red_3 module might play an important role in regulating the balance between cell development and the synthesis of nutrients. The Grey_60 module also showed specific enrichment and appeared to be highly associated with the glycerolipid metabolic process, glycerophospholipid metabolic process and phospholipid metabolic process (Fig. 6B, Table S8). These biological processes and the starch biosynthesis process both utilize fructose-6-P as a substrate for further extension, which suggests that the Grey_60 module likely affects the relationships between sucrose decomposition and starch synthesis through substrate competitive inhibition. Additionally, it should be noted that the 10 modules exhibited a high correlation and formed a cascade-based community, and a functional analysis revealed that genes in the 10 modules were highly significantly enriched in the reproductive process (Fig. 6B, Fig. S7). This result indicated that the reproduction process might be the most fundamental biological process and could further affect other essential metabolism processes, including the starch synthesis process, by interacting with multiple genes belonging to different modules during the early stages of maize endosperm development.

Identification of hub genes and key networks associated with starch synthesis

Most maize agronomic traits and metabolic processes are complex due to not only the actions of single genes but also interactions between combinations of genes. A sizeable portion of the genes in each network module exhibited extremely high connectivity with other genes belonging to the 10 modules in which the core genes in starch synthesis were distributed. These highly connected genes are usually considered hub components of a network and have been shown to be biologically important in multiple previous applications. Therefore, we performed a NetworkAnalyzer analysis and found that 5% (918) of the genes in the 10 modules had a higher degree of connection and could thus be identified as hub genes for further study (Table S9). A functional enrichment analysis of these hub genes indicated that many of these genes encode proteins functioning in chromatin and protein binding, oxidoreductase activity and transferase activity (transferring hexosyl and glycosyl groups). In addition, these hub genes mainly play an important role in the response to hormone stimulus, regulation of the cellular component size, system development, heterocycle catabolic process, and response to water deprivation (Table S9).

According to the gene classification results, 77 transcription factors of 918 hub genes were identified, and these belonged to distinct families, some of which (AP2/ERF, ARF, bHLH, MYB, WRKY, NAC, bZIP, and GRAS) have been reported to play roles in cell development and the synthesis and accumulation of starch in maize, rice, hulless barley, Arabidopsis, and other plants (Iwase et al., 2011; Licausi, Ohme-Takagi & Perata, 2013; Liu et al., 2015; Huang et al., 2016; Butardo Jr et al., 2017; Cai et al., 2017; Tang et al., 2017; Xiao et al., 2018) (Table S9). In particular, only one core gene of starch synthesis (GRMZM2G158043, ZmPUL) was identified as a hub gene in the Grey_60 module, and this gene is mainly involved in the modification of excessively branched chains or the removal of improper branches of amylopectin formed by branching enzymes to maintain the cluster structure of amylopectin. More remarkably, beta-glucosidase 18 and aldose 1-epimerase (encoded by GRMZM2G100452 and GRMZM2G031628, respectively) were also found in the Grey_60 module and are mainly related to carbohydrate transport and metabolism (Table S8). These results indicated that some hub genes of the Grey_60 module might play a greater role in regulating the synthesis of amylose and amylopectin. However, this finding raises the question of whether other core genes of starch biosynthesis are key factors that determine the contents of amylose and amylopectin.

To further explore the roles of core genes of the starch synthesis pathway in regulating the contents of amylose and amylopectin, we extracted the coexpression network of 19 differentially expressed core genes of starch synthesis. The analysis revealed that starch biosynthesis is regulated by multiple cascade networks that fragment into more sub-networks, and some of the sub-networks of the core genes of starch synthesis, particularly the sub-networks centred on ZmPUL, were linked by one or more linking genes (Fig. 7). Other disjoint sub-networks might rely on more linking genes as bridges in the regulation of starch synthesis (Fig. 7). Additionally, the relationships of the genes in the community network were also further supported by those of the protein-protein interactions in the STRING database (https://string-db.org/) presented in a previous study (Qu et al., 2018). A rigorous functional enrichment analysis indicated that these linking genes play crucial roles in cell cycle regulation, replication, the MAPK cascade, energy metabolism, binding (including protein binding, nucleotide binding and ion binding), carbohydrate transport and metabolism, lipid transport and metabolic processes, amino acid transport and metabolism, and translation and posttranslational modification (Table S10). These results further indicated that the above-mentioned biological processes might be intimately linked to starch synthesis, and the specific functions of core genes of starch synthesis might be affected by regulatory genes of these metabolism processes. Nevertheless, the detailed regulatory and interactive mechanisms between core genes of starch synthesis and genes involved in other key biological processes remain to be elucidated.