Dramatic morphology and reserve changes during seed development
Whole developing soybean seeds at four sequential time points (20, 25, 30, and 40 DAF) covering the major stages of seed development were sampled in triplicate (Fig. 1A). This period of seed development is coincident with seed growth stages from R5 (beginning seed) to R6 (full-size seed) of soybean development [23], which are morphologically comparable to early and mid-maturation stages as previously described [7]. This period also co-occurs with the rapid accumulation of lipids and storage proteins [20, 24]; thus, the seed size and seed weight increased significantly as it grew toward a fully developed seed (Fig. 1A). Two major storage reserves in soybean seeds, namely, protein and oil, were quantified during this period. As shown in Fig. 1B, the accumulation of stored oil was approximately 14.11% at D20 and continuously increased to 17.85% at D40 with the greatest increase of 1.85% occurring in the period from D20 to D25 relative to the later periods of D25-D30 by 1.11% and D30-D40 by 0.78%. In contrast, as the seeds matured, the stored protein content was 50.87% at D20, which decreased to 45.83% (5.04% decrease) at D25 and continuously declined to 42.26% at D30 and to 41.77% at D40. Thus, the oil content in the developing seeds increased more dramatically in the early period and less dramatically in the late period of seed filling, and vice versa for the protein decrease. The reduction in the relative percentage of protein content over the period is due to its negative correlation with the oil content in soybean seeds [25]. In contrast, the actual protein levels increased rapidly between 10 and 40 DAF [26]. Hence, this period of seed development presents rapid morphological changes (such as size and weight) and major metabolic changes (such as oil and protein), which should be governed by a variety of cellular processes/gene activities in seeds [17, 20].
Transcriptomic abundance shifts during seed development
An in-depth understanding of biological processes associated with seed development and reserve accumulation in soybean seeds is important for improving seed quality and yield potential. Given the multiplicity of genes controlling seed development coupled with seed filling, time-course transcriptomic profiling of four-stage seeds that covered the representative period of seed development was performed with the aim of obtaining a global view of the regulatory mechanism underlying this complex development process. In total, ~ 313.4 million sequencing reads from all sampled seeds were generated, with an average of 26.1 million reads per sample (Table S1). On average, over 97.43% (25.4 million) of reads were kept after trimming the low-quality reads and used for read alignment. Uniquely mapped reads were used to quantify gene expression abundance as normalized by counts per million (cpm). After removing reads/transcripts < 5, we determined that, on average, 35,000 (62.5%) out of 56,044 gene models were expressed in developing seeds. We observed a gradient reduction (approximately 700 genes, 2% of the total) in the number of expressed genes over the course of seed development. The reduction could be, in part, related to the shutdown of many transcription processes as the seeds grow toward dormancy [27]. A previous study has reported a programmed increase in methylation levels during soybean seed development [7]. Thus, it would be interesting to determine whether the increased methylation level is related to the reduction of the seed transcriptome.
Both clustering dendrogram analysis and multidimensional scale analysis indicated that gene expression levels among the replicated samples were generally highly related to those between time-point samples (Fig. 2A, B), indicating gene expression specificity to the time points. For the transcript abundance, approximately 61.07% of expressed genes were in the range of 0–10 cpm, 36.41% were in the range of 10–100 cpm, and only 2.53% of genes were highly expressed, at cpm ≥ 100 (Fig. 2C). As the seeds enlarged, the number of genes with an abundance ≥10 cpm generally decreased, and the number of genes with an abundance ≤10 cpm gradually increased. Among the expressed genes, 32,326 were expressed across all four stages of developing seeds. There were 942, 329, 200, and 546 genes specifically expressed at D20, D25, D30, and D40, respectively (Fig. 2D). These results suggest a temporal shift in transcript abundance for many highly expressed genes as seeds developed, which might reflect the foregoing morphological and/or metabolic changes during seed development.
A cluster dendrogram analysis showed that replicates from each condition were closely clustered as supported by a similar separation pattern in a multidimensional scaling (MDS) plot of log-FC values. We further performed qRT-PCR on 16 randomly selected DEGs (Table S4) for validation, and a high correlation (r = 0.87) between the RNA-seq and qRT-PCR results indicated the robustness of our RNA-seq results.
Genes with high transcript abundance during seed development and filling
We next characterized the genes that were constitutively highly expressed across all four time points. In total, 443 genes with a mRNA abundance cpm ≥ 100 were identified [28] (Table S3). GO term enrichment analysis indicated that enriched biosynthetic processes (BPs) essential for nutrient accumulation were amino acid metabolism, carbohydrate metabolism, glycolysis, and hexose metabolism (Table S3). This gene set was also enriched for the molecular function (MF) term nutrient reservoir activity and cellular component (CC) terms related to lipid storage bodies, endoplasmic reticulum (ER), and ribosomes, sites where fatty acids are primarily synthesized and stored [29, 30]. When we relaxed the threshold to cpm ≥ 50, the number of newly defined highly expressed genes increased by threefold to 1386 genes (Table S3). Accordingly, the number of expressed genes in the majority of enriched terms was approximately doubled or tripled, such as hexose metabolic process (8 to 22), glycolysis (7 to 18), cellular nitrogen compound metabolic process (18 to 36), cellular amino acid biosynthetic process (15 to 32), and ribosomal proteins (47 to 132). It should be noted that no change in the number of genes with terms related to monolayer-surrounded lipid storage bodies (7), nutrient reservoir activity (11), and the sulfur amino acid metabolic process (6 to 7) was observed, suggesting the importance of maintaining high transcript abundances of these genes for seed filling. Additional enriched GO terms for the 1386 genes identified were as follows: protein localization, glutamine family amino acid biosynthetic process, intracellular transport, gene silencing by RNA, and vesicle-related components (COP-II-, COP-I, coated vesicles, transport vesicles, ER to Golgi transport vesicles, and membrane-bound vesicles).
A detailed investigation of the 443 genes (cpm ≥ 100) revealed several with demonstrated roles in the accumulation of storage reserves, such as oil, protein, and carbohydrates. Not surprisingly, some genes participating in the metabolism of protein and oil were identified. However, several of the highly abundant genes detected during seed filling have rarely been reported (Table S3), such as the oil accumulation-related genes, FAD2 (467.83–2000.69) catalyzing oleic acid accumulation [31], OLEOSIN proteins (373.03–6190.57 cpm) encoding the oil body surrounding lipid storage bodies [30, 32], and metabolic genes for the biosynthesis of leucine, thiamin (vitamin B1), and ribosomal proteins. Genes encoding storage proteins, such as cupin family proteins, showed the highest transcript abundance, ranging from 34,000–100,000 cpm, in developing seeds followed by xylem bark cysteine peptidase 3 (Glyma.08G116300 with 28,399.81 cpm), trypsin inhibitor 1 (222.79 ~ 18,559.51 cpm), and oleosin family proteins, which were highly abundant (373.03–6190.57 cpm). These proteins in soybeans appear to be redundant and serve as a surfactant for the oil body [33, 34]. A recent study on oleosin on chr20 [32] has indicated that it is likely regulated by LEC2 [15], by which oil production is enhanced through increasing oil body turnover. Investigation of the cis-regulatory sequences in the promoter region could be helpful in using this set of genes to achieve oil content improvement. In contrast, a set of highly abundant trypsin inhibitor genes (222.79–18,559.51 cpm) could be silenced by gene editing to inhibit their expression or mutate them to achieve better seed quality because trypsin inhibitors have an antinutritional effect [35].
More attention should be given to nutrient transporters that mediate nutrient flux into the developing embryo, such as the SWEET sugar transporters (274.98–943.19 cpm), proton pump interactor 1 and SUCROSE SYNTHASEs (121.5–1,186.21 cpm) controlling starch synthesis. Recent studies have highlighted the essential role of SWEET proteins in transporting sugar from maternal to embryonic tissues to support seed development and oil accumulation [12, 36, 37], indicating the agriculturally critical role of causing high oil accumulation in cultivated soybean seeds relative to the low oil content in wild soybeans [11]. High abundances of these transporters, together with related studies, suggest that the rate of sugar transport and sugar accumulation is important for nutrient reserve accumulation in seeds. Further investigation of downstream affected steps/components may help in identifying the major route of carbon flow for oil accumulation or carbon partitioning between oil and proteins in soybeans and, perhaps, other plant species, such as maize and rice [38, 39]. Several LOX1 genes (541.27–16,226.05 cpm) had high transcript abundances in the present study, but homologous proteins have not been detected during seed filling in wheat and maize endosperms [40, 41], suggesting a lineage-specific role of LOX proteins in soybean seeds, possibly related to seed storage proteins and lipids, such as unsaturated fatty acids [42, 43], or special tastes and odors [44].
Transcription factors, such as ABA INSENSITIVE 3 (ABI3; Glyma.08G357600; Glyma.18G176100) and its interacting regulator, bZIP67, were identified in high abundance. The high abundance of both regulators may be due to the increasing need for their involvement in cell growth and the accumulation of storage reserves as seeds enlarge. Indeed, both regulators have common functions and are essential for initiating gene expression changes for biological processes critical for seed development, such as gibberellic acid signaling, photosynthesis, and maturation-related processes [16, 17, 45, 46]. LEC2s (cpm = 0) were undetectable in this study, possibly because they are highly specific to certain seed compartments and barely detected in whole seeds [15]. In addition, well-documented genes homologous to LEC1, which is a central regulator of seed development [17], were not in high abundance but showed decreased expression over time. Other high-abundance regulators with demonstrated roles in Arabidopsis that are less explored in soybean seeds were also identified, such as PLASTID TRANSCRIPTIONALLY ACTIVE 9 (PTAC9), TANDEM CCCH ZINC FINGER PROTEIN 4 (TCZF4) and circadian rhythm regulators, including TIME FOR COFFEE (TIC), LATE ELONGATED HYPOCOTYL 1 (LHY1), and CONSTANS-LIKE 4 (COL4). Determination of the functions of the regulators and the cognate cis-regulatory elements may facilitate the identification of seed regulatory networks that have contributing roles in the regulation of genes temporally and spatially in developing seeds [17, 18, 47]. All highly expressed genes and the functional annotations are listed in Table S2.
Genes are temporally regulated during seed development
To comprehensively identify genes involved in developmental processes during seed filling, we conducted a pairwise comparison of transcriptomic profiles among the different developmental stages of seeds. We observed a gradient increase in the number of DEGs between the two neighboring time points as the seeds developed, with 523 in D25 vs. D20, 1486 in D30 vs. D25, and 2485 in D40 vs. D30 (Fig. 3). When comparisons were made across nonadjacent time points (over 10-day intervals), including 4079 in D30 vs. D20, 5617 in D40 vs. D25, and 7670 in D40 vs. D20, the numbers of DEGs increased. In total, the pairwise comparison identified 8925 nonredundant DEGs that are required to program soybean seed development from early to full seed development (Table S2).
Clustering analysis was performed for all 8925 DEGs based on expression and is shown in a heat map (Fig. 4). The analysis classified the DEGs into seven groups (A to G) exhibiting different gene expression patterns. The difference in the temporal expression patterns suggested stage-specific roles of these grouped DEGs throughout seed development. Groups A and G represented the two largest clusters (2182 and 4945 DEGs, respectively) while exhibiting contrasting expression patterns. The DEGs in cluster A displayed the predominant and highest expression levels at D40, and those in group G showed a continuous decrease in gene expression as seeds matured. Enrichment analysis for cluster A genes showed that D40-specific genes were enriched for cell communication (GO:0007154) and transcription (GO:0006351), and the enriched biological processes occurred in the cell periphery and plasma membrane (Fig. 4B). This stage represented the later stage of reserve production and storage as well as the active period of the shift to maturation in which multiple programming might be initiated, such as cessation of organ differentiation, cessation of cell division, and preparation for seed desiccation [48]. A representative gene with specifically high expression at D40 is ABI3 (Fig. 4C) [16], which is an essential regulator in abscisic acid (ABA) signaling primarily participating in seed maturation and dormancy in Arabidopsis [49]. The high abundance of AB13 during D20-D30 and dramatic increase toward D40 suggested its critical role during seed development, especially in the later stage of maturation-related processes. ABI3 has recently been demonstrated to be a key transcription factor modulating the transition from morphogenesis to maturation by combinational interaction with other seed transcription factors [17]. Further study of the coexpressed genes or binding sites [17, 21] would provide an opportunity to understand regulation in soybean seeds. Identification of the key seed development regulators (such as ABI3) indicated the robustness of our analysis and suggested that other seed development-associated genes identified in our study deserve further attention. For example, a similar expression pattern was also observed for two sucrose synthases (SUSs) involved in cellulose and starch biosynthesis in seeds [50].
Genes classified in cluster B (810) were highly expressed at D40 and had a similar expression pattern to that in cluster A, except for those genes that were upregulated beginning at D30. Enriched GO terms for cluster B were transmembrane transport (GO:0055085) and nutrient reservoir activity (GO:0045735), suggesting that biological processes associated with nutrient transport, such as sugar transporters, amino acid transporters, and nutrient accumulation, were active beginning at D30 and maintained at D40. The developmental period of D30-D40 is approximately coincident with the middle-maturation stage, which represents the major period of accumulation for storage, such as that of protein and lipids [7]. Expression pattern analysis showed that despite the varying expression patterns observed, genes involved in the aforementioned biological processes showed obvious upregulation during seed development and filling (Fig. 4C), such as cupin proteins, oleosins, amino acid permeases (APPs), FAD2, nitrate transporters, and seed maturation proteins (SMPs). Such stark upregulation in expression undoubtedly reflects the critical roles of the related biological processes in storage reserve accumulation. Transcription factors involved in seed development were also upregulated in cluster B, such as the AUS/IAA transcription factor, YABBY, and uncharacterized PLATZ. Interestingly, the genes related to flowering (TFLs, Tof12, and PHYA) showed an upregulation/downregulation shift during the course, suggesting expanded roles during seed development, but this idea has yet to be explored. This finding is partially supported by a recent study demonstrating that Arabidopsis TFL1 not only functions as a shoot identity gene but also serves as a signaling molecule essential in determining seed size, and the tfl1–20 loss-of-function mutant has been shown to produce larger developing seeds than wild type through regulation of endosperm cellularization [51]. The expression patterns for these DEGs during seed development are illustrated in Fig. 4C.
In contrast, the top GO terms enriched for the genes in cluster G are mainly related to photosynthesis (GO:0015979) and microtubule-based processes (GO:0007017). Downregulation of these gene sets mainly reflect the reduction of the photosynthetic process as seeds developed, which as evidenced by downregulation of the bulk of genes encoding varying processes of photosynthesis, such as chlorophyll A-B binding family proteins, high-chlorophyll fluorescent proteins, and light-harvesting complex photosystem II (Table S1). Other family genes exhibiting downregulated expression patterns were also identified, such as GA2OXs, motor proteins, LACS9s, stay green-like genes, and LEC1-like genes (Table S1). The decrease in expression of these genes with diverse roles associated with phytohormones (such as GA) [52], chlorophyll breakdown during senescence [53], and seed developmental processes [47] might correspond to multiple processes of gradual cessation of cell division or commencement of maturation as seeds develop and mature. Clusters C-F comprised relatively fewer DEGs that displayed relatively high expression levels at two of the four selected time points, and no GO term was significantly enriched.
DEGs exhibit spatial expression patterns in seed compartments and nonseed tissues
We next examined the spatial expression profiles of all DEGs in different developmental tissues using the Harada-Goldberg soybean RNA-seq dataset (http://seedgenenetwork.net). This dataset contains 10 tissues of Williams 82, which is the same accession used in our study, including 6 reproductive tissues [floral buds (FB) and five stages of developing seeds, including whole-mount globular (GLOB), heart (HRT), cotyledon (COT), early-maturation (EM), and dry seed (DS)] and four vegetative tissues (leaf (LF), root (RT), stem (STM), and seedling (SDL)). Although this study was performed independently, it was complementary to ours, and integration of our results with spatial analysis using multiple tissues may provide additional insight into the mechanism.
The DEGs exhibited clear tissue-specific (TS) patterns, as indicated by the higher expression in certain tissues compared to others. This analysis grouped the DEGs into seven large (TS1-TS7) groups with each containing over 400 genes, and the remaining three groups (not shown) comprised less than 30 genes (Table S2). The seven large groups, which were further investigated, are shown in Fig. 5A. Enrichment analysis for the seven TS groups revealed that no common top five enriched GO terms, indicating possible certain roles in the highly expressed tissues. Interestingly, genes that were highly expressed in developing seeds were divided into three groups: TS1 grouped genes highly expressed in three early-stage developing seeds, such as GLOB, HRT, and COT; TS2 grouped genes highly expressed in the EM-stage seeds; and TS3 grouped genes highly expressed in the later-stage seeds and DS-stage seeds. Approximately 2647 (29.7%) of the DEGs were clustered into TS1, consisting of genes primarily expressed in GLOB-, HRT-, and COT-stage seeds, and these genes showed continuous expression with either upregulation or downregulation. The majority (79.1%, 2094 of 2647) of TS1 was classified into the G group (downregulation group) (Fig. 5A, C), exhibiting downregulation as seeds grew. Genes (585, 6.6% of total) in TS2 were predominantly expressed in the EM-stage seeds. In total, 970 (10.9%) DEGs were exclusively expressed in DS-state seeds in TS7. This spatial difference in expression among the different stages of seeds suggested the possible different roles of the two gene sets in seeds. Indeed, TS1 genes are mainly involved in processes related to the cell cycle and DNA replication, which might account for early seed compartment morphogenesis and growth (Fig. 5B), and TS2 (EM) was more specialized in metabolic processes, such as fatty acid biosynthesis and metabolism of carboxylic acids and organic acids. In contrast, genes highly expressed in DS-stage seeds participate in the response to environmental stimuli, such as temperature, heat, alcohol, and plant hormones, which regulate plant growth, such as abscisic acid, and over half (67.7%) of this gene set was upregulated (Fig. 5C). For example, LEA14 (LATE EMBRYOGENESIS ABUNDANT 14) from TS7 is involved in wounding and light stress as well as possible protection against desiccation [54, 55]. Thus, the DEGs highly expressed in DS-stage seeds might be partially related to preparing transcripts for dormancy and, possibly, seed germination under favorable conditions.
Recent studies have indicated that the expression of many of these seed-expressed genes is highly spatially specific in seeds, and this specificity represents functional relevance. For example, seed coat-specific SWEET39 is involved in sugar delivery from maternal plants to developing embryos [56] and endosperm-specific TFL1 in endosperm cellularization, which affects seed size [51]. In addition, endosperm-exclusive SWEET15 participates in seed development by regulating sugar delivery into developing embryos [12]. These studies suggest that many of the genes highly expressed in seed compartments, as identified here (Table S2), can be prioritized for detailed investigations to gain a more comprehensive understanding, including YABBA, GA3OX1, NF-YB6, NF-YA9, FAD8, KAS I/III, SWEETs, OLEOSINs, AATs (amino acid transporters), and uncharacterized LEC1-LIKE (Glyma.07G268100, Glyma.17G005600), which serves as a key regulator of fatty acid biosynthesis in Arabidopsis [57].
In total, 1357 (15.2%) DEGs clustered into TS3 and were mainly expressed in leaves; they are involved in many aspects of photosynthesis with over 80% of the DEGs downregulated as seeds matured (Fig. 5C). These results indicated that this gene set maintain photosynthetic roles and is also associated with seed developmental processes. It is likely that the expression of many of these photosynthesis-related genes in seeds might be influenced by mosaic interactions of different transcription factors, such as LEC1, ABI3, bZIP67, and AREB3 [17, 47]. Genes that were highly expressed in the RT and STM shared many DEGs in common and clustered in TS5, which agreed with the tissue functionality of both being involved in nutrient translocation [58]. These genes are mainly involved in carbohydrate and sugar metabolism as well as nutrient and ion transport to meet the large demand for nutrient uptake and transport (Fig. 5B). TS6 was the smallest group (400), with genes highly expressed in SDLs and participating in processes such as oxidation-reduction. These results indicated that genes that have relatively high expression in nonseed tissues, such as LF, FB, and RT, might also maintain tissue-represented functions in developing seeds (Fig. 5B); however, their expression is relatively low in developing seeds, suggesting their indispensable roles in seed development. For example, AAP2 from TS6 (stem and leaf) affects the levels of nitrogen and carbon assimilation as well as seed yield by altering xylem-phloem transfer in amino acid transfer [59], and FRD3 (MATE efflux family protein) from TS6 functions in the root xylem for efficient iron uptake from the xylem into leaf cells [60]. Therefore, the identification of these genes with nonseed tissue specificities in developing seeds was reflective of the necessity of their basic functionalities in maintaining the needed processes to support seed development. Further investigation of these genes may provide deep insight into their supportive roles during seed development.
As indicated in this and other studies [17, 61], multiple layers of gene activities are involved in seed development, and the expression of the genes is also subjected to environmental effects. Despite promising genes identified from integrative analysis, it is ideal to compare the transcriptomes in different seed compartments and nonseed tissues at each of the time points in the same environment, which would offer a more comprehensive understanding of gene activities in developing seeds over the period.
Coexpression network-inferred seed regulatory network
We next extracted the coexpression connections of all DEGs from a global SoyNet database (http://www.inetbio.org/soynet/) constructed with transcriptomes from different seed compartments at different developmental stages, along with other tissues [47, 62]. Only a small fraction of the DEGs could be reconstructed into modules. In total, 134 modules were identified with the numbers of nodes (genes) ranging from 2 to 151. After color assignment based on the aforementioned TS groups, 11 modules consisting of relatively tissue-specific patterns were identified (Fig. 5D). It is interesting that many genes from TS1 (GLOB, HRT, and COT) and TS3 (LF) were clustered together as individual TS-specific modules compared to the other TS groups. Two modules specific for TS1 were identified. TCX2 (TESMIN-LIKE CXC2) controls stem cell division by regulating stem cell type-specific networks [63]. TCX2 slightly deviated from the main module, suggesting its role in orchestrating the connected genes of this module for coordinated division of different cell types for seed development, which is in agreement with the function of GLOB/HRT/COT, where ontogenesis of different seed compartments occurs during this stage (http://seedgenenetwork.net/). Similarly, TDM1 plays an important role in meiosis termination [64], and it was strongly associated with genes in the TS1-dominated modules. We determined that uncharacterized TKL2 (TRANSKETOLASE 2) might be a key connector bridging two modules comprising LF and GLOB/HRT/COT genes. It is likely that the link between RBCS3B (RUBISCO SMALL SUBUNIT 3B) and BLH1 (EMBRYO SAC DEVELOPMENT ARREST 29) connected the two individual modules. These results suggest a functional relevance of these modules in seed development, and these genes that are underexplored in soybeans deserve more attention.
