The profile of seed abortion in S. grosvenorii
The structural changes in the tissues at each stage were observed using TEM. There was no obvious edema in the cytoplasm and no obvious separation of the plasma wall at stage 5DAF (Figure 1a and b). The nucleus was irregular, the mitochondria were oval, the crest and rough endoplasmic reticulum were moderately dilated, the surface ribosome was locally exfoliated, and the number of intracellular starch granules was abundant. Severe edema of the cytoplasm and slight separation of the plasma wall were observed at stages 10DAF (Figure 1c and d) and 15DAF (Figure 1e and f), respectively. Severe separation of the plasma wall, local damage to the cell membrane, and thickness of the cell wall were observed at 20DAF (Figure 1g and h) and gradually recovery was observed from stages 25DAF (Figure 1i and j) to 30DAF (Figure 1k and l).
We also examined the levels of endogenous plant hormones, ABA, ZR, IAA, and GA3, during seed abortion using ELISA (Figure 2a-d). Interestingly, the levels of ABA and ZR decreased remarkably at stages 15DAF and 20DAF, respectively. The above changes in submicroscopic structure and levels of endogenous plant hormones indicated that 15–20DAF was a key transition period during seed abortion. We will mainly focus on the exploration of the transcriptomic changes at stages 15DAF and 20DAF in the subsequent analyses.
Global analysis of the RNA-Seq data
Genes with significant differential expression (|log2FC|>1 and FDR-adjusted p-value < 0.05) were identified between consecutive stages. The number of DEGs in the five comparisons is shown in Figure 3a and b. The total DEGs in 10DAF vs 5DAF were 7,692 (4,145 upregulated and 3,546 downregulated), 15DAF vs 10DAF were 3,252 (1,487 upregulated and 1,765 downregulated), 20DAF vs 15DAF were 984 (364 upregulated and 620 downregulated), 25DAF vs 20DAF were 2,124 (942 upregulated and 1,182 downregulated), and 30DAF vs 25DAF were 885 (258 upregulated and 627 downregulated).
The DEGs were further functionally classified into GO slims (Tables S1–S5; p < 0.05). DEGs in 20DAF vs. 15DAF were identified to be involved in meristem structural organization, morphogenesis of a branching structure, adaxial/abaxial axis specification, positive regulation of catalytic activity, oxidoreductase activity, and glucosidase activity (Table S3), whereas those in 25DAF vs. 20DAF were involved in trehalose metabolic process, seed germination, response to sucrose, transmembrane receptor histidine kinase activity, and so forth. (Table S4).
Seven common DEGs shared by all comparisons of adjacent stages were identified (Figure 3b), which revealed fluctuations in expression levels across all stages (Figure 3c). These DEGs were found to be associated with toxin metabolic processes, regulation of gene expression, glucan metabolic process, and external encapsulating structure organization.
Identification of temporal expression trends across S. grosvenorii seed transcriptomes
We performed STEM analysis to visualize the expression patterns of mRNAs. As a result, seven expression profiles were found to be statistically significant (Figure 4a). Profiles 0 and 19 tended to be continuously downregulated and upregulated during seed abortion, respectively, whereas profiles 1 and 18 reached an expression peak and valley at stages 15DAF and 20DAF, respectively.
KEGG pathway enrichment analysis was performed for the four expression trends (19, 0, 1, and 18). The genes in profile 0 with a continuously downregulated trend were mainly enriched in starch and sucrose metabolism, plant hormone signal transduction, linoleic acid metabolism, replication and repair, and other glycan degradation (Figure 4b), whereas those in profile 19, which were continuously upregulated, were related to amino acid metabolism; phenylpropanoid biosynthesis; ubiquinone and other terpenoid-quinone biosynthesis; flavonoid biosynthesis; stilbenoid, diarylheptanoid, and gingerol biosynthesis; diterpenoid biosynthesis; glycolysis/gluconeogenesis; and plant-pathogen interactions (Figure 4c). The genes in profile 1 with the lowest expression level at 15DAF and 20DAF were involved in RNA polymerase, butanoate metabolism, inositol phosphate metabolism, and fatty acid degradation (Figure 4d). Finally, the genes in profile 18 with the highest expression levels at stages 15DAF and 20DAF were involved in arachidonic acid metabolism, amino sugar and nucleotide sugar metabolism, one carbon pool by folate, cyanoamino acid metabolism, glutathione metabolism, isoquinoline alkaloid biosynthesis, phagosome, glycerolipid metabolism, and fructose and mannose metabolism (Figure 4e).
Co-expression network associated with seed abortion
A gene co-expression network was used to cluster 21,152 genes after filtering those with low expression levels, resulting in 16 modules with different colors (Figure 5a). Combined with the module-trait relationship and module significance, we eventually identified the brown2 module as the most phenotypically relevant module, which showed a negative correlation with the level of ABA (the correlation coefficient was -0.71 and the FDR-adjusted p-value was 9e-04), whereas the lightsteelblue module showed a positive correlation (Figure 5b). The indianred4 and skyblue3 modules were also found to be negatively and positively associated with the level of ZR, respectively, with statistical significance (Figure 5b).
We explored the expression profiles of each module across the different stages (Figure 6a). Interestingly, we found that the genes in the brown2 module were specifically, highly expressed at 15DAF and 20DAF (Figure 6a). By combining the above results, we considered the brown2 module to be negatively associated with the changes in ABA levels at 15DAF. According to the results of GO enrichment analysis, the genes in the brown2 module were enriched in monoterpene metabolic processes and terpene metabolic processes (Figure 6b). Additionally, KEGG pathway analysis revealed that the genes were involved in sesquiterpenoid and triterpenoid biosynthesis and amino sugar and nucleotide sugar metabolism (Figure 6c). These findings indicate that these genes, which are related to terpene metabolism, are highly expressed at stage 15DAF and may play an important role in the decrease of ABA levels, thereby leading to seed abortion.
Finally, the co-expression networks of the brown2, indianred4, lightsteelblue, and skyblue3 modules were filtered by a weight value greater than 0.15 and visualized (Figure 7a-d). The genes with the highest betweenness centrality score in the brown2 module were evm.TU.tig00004479.11 and evm.TU.tig00153447.44, which are non-specific lipid-transfer protein 2-like (LTP) and protein root initiation defective 3-like, respectively (Figure 7a). The genes with the highest betweenness centrality score in the indianred4 module included CCR2, PCBER, and CYP75B2, which were cinnamoyl-CoA reductase 2-like, isoflavone reductase-like protein IRL, and cytochrome P450 71A1-like, respectively (Figure7b). The genes with the highest betweenness centrality scores in the lightsteelblue module were beta-galactosidase 10 (Os01g0875500-2, Os01g0875500-3, and Os01g0875500-1) (Figure7c). Finally, the genes with the highest betweenness centrality score in the skyblue3 module were calcium-transporting ATPase 12 plasma membrane-type-like, ACA12-1 and ACA12-2, and the ethylene-responsive transcription factor, ERF105 (ERF5) (Figure 7d).