Transcriptomic Analysis Reveals Key Candidate Genes Related to Seed Abortion in Chinese Jujube (Ziziphus jujuba Mill).

Background: Seed abortion is a common phenomenon in Chinese jujube that seriously hinders the process of cross-breeding. However, the molecular mechanisms of seed abortion remain unclear in jujube.
Methods: Here, we performed transcriptome sequencing using eight flower and fruit tissues at different developmental stages in Ziziphus jujuba Mill. 'Zhongqiusucui' to identify key genes related to seed abortion. Histological analysis revealed a critical developmental process of embryo abortion after fertilization.
Results: Comparisons of gene expression revealed a total of 14,012 differentially expressed genes. Functional enrichment analyses of differentially expressed genes between various sample types uncovered several important biological processes, such as embryo development, cellular metabolism, and stress response, that were potentially involved in the regulation of seed abortion. Furthermore, gene co-expression network analysis revealed a suite of potential key genes related to ovule and seed development. We focused on three types of candidate genes, agamous subfamily genes, plant ATP-binding cassette subfamily G transporters, and metacaspase enzymes, and showed that the expression profiles of some members were associated with embryo abortion.
Conclusion: This work generates a comprehensive gene expression data source for unraveling the molecular mechanisms of seed abortion and aids future cross-breeding efforts in jujube.

INTRODUCTION

Chinese jujube (Ziziphus jujuba Mill.), a member of the Rhamnaceae family, is thought to have been cultivated for more than 7,000 years and is one of the most economically important fruit trees in China [1]. It is widely distributed and rich in germplasm resources [2]; it also has high nutritional and medicinal value, as it contains several biologically active components required by the human body, including vitamin C, phenolics, flavonoids, triterpenic acids, polysaccharides, and microelements [3, 4]. Jujube fruits can be consumed fresh, dried, or processed in various forms [5]. In addition to its high economic value, the Chinese jujube also provides ecological benefits. Because of its ability to grow in diverse environments, the jujube tree has played a critical role in restructuring the agricultural industry and in the economic development of poverty-stricken regions in China. Some of the disadvantages associated with the production of Chinese jujube fruit include the severe abscission of flowers and fruits and the ease with which the fruit can crack. Although improvements in jujube varieties are essential for enhancing the quality, yield, and efficiency of jujube cultivation [6], embryo abortion is a common phenomenon in jujube, which prevents the formation of viable seeds, restricts cross-breeding efficiency, and makes obtaining hybrid progeny a major challenge. There is thus an urgent need to understand the causes of seed abortion at the molecular level, as such information is critically important for breeding new cultivars specifically and improving the jujube industry more generally.

The factors affecting plant seed abortion have been studied extensively. In model plants, genetic studies have greatly improved our understanding of embryo development, and mutations of genes related to embryo development have been shown to cause seed abortion. In Arabidopsis thaliana, a great number of genes required for seed development have been identified through large-scale mutagenesis approaches [7, 8], which have provided valuable resources for functional analyses. Biotic and abiotic stresses, the partitioning of carbohydrates, and phytohormones (e.g., auxin, gibberellin) can affect the development of seed and fruit, and even lead to seed abortion [9]. Studies on invertases, enzymes hydrolyzing sucrose into glucose and fructose, have revealed that the elevated expression of invertase and high glucose concentrations enhance embryo development and reduce seed abortion in maize and tomato [10, 11]. In fruit crops, studies on seedlessness have identified many genes that may be selected during domestication. For example, in grapevine, VvAGL11 (a homologous gene of Arabidopsis Seedstick) has been found to play a central role during the selection of seedless grapevine cultivars [12]; sequence variations in the promoter and coding regions of VvAGL11 have been shown to be correlated with seed abortion [13].

In recent years, the advent of high-throughput ‘omics’ studies has greatly facilitated our understanding of the molecular regulation of seed abortion at the transcriptional, biochemical, and metabolic levels. Based on large-scale transcriptomic and proteomic analyses, 140 differentially expressed genes (DEGs) and 41 differentially expressed proteins were identified between normal and aborted embryos in chrysanthemum [14]. Comparative proteomic analysis in longan (Dimocarpus longan Lour.) has yielded more than 1000 differentially expressed proteins [15]. In peanut, transcriptomics studies have also been used to study embryo abortion and have identified thousands of DEGs, including genes involved in cell division, the stress response, embryonic development, and signal transduction [16, 17]. These results have highlighted the complexity of the regulation of embryo development, which involves the coordination of multiple pathways; more generally, these findings have emphasized the need for a deeper understanding of the molecular functions of critical genes.

Previous studies have primarily focused on the physiology, biochemistry, and embryology aspects of the mechanisms of seed abortion in jujube [18-20]. By comparison, little is known of the molecular mechanisms. Z. jujuba Mill. ‘Zhongqiusucui’ is a newly bred cultivar with large, highly sweet fruits [21]. We have previously observed defects in the abnormal embryo sac, despite the fact that the development of microsporogenesis appeared normal, which led to incomplete embryo development [22]. Here, we used RNA-seq technology to sequence and analyze the transcriptome of jujube flowers and fruits at different developmental stages. A great number of DEGs were identified for further functional analysis, which revealed a suite of genes potentially involved in seed abortion. Generally, this work provides a foundation for understanding the regulation of seed development in jujube.

MATERIALS AND METHODS

Plant Materials

The experimental jujube cultivar was Z. jujuba Mill. ‘Zhongqiusucui’, which was selected from a bud mutation of Z. jujuba Mill. ‘Tangzao’ in Qidong County, Hunan Province, China [21]. Its flowers are hermaphroditic and are composed of the flower stalk, calyx, petal, flower disk, stamen, and pistil. It has a high seed abortion rate (approximately 95%) under natural pollination according to our previous study (unpublished data).

The experimental site was the jujube tree experimental base of Central South University of Forestry and Technology (Qidong County, Hunan Province; 112.12°E, 26.78°N). There were about 300 five-year-old jujube trees that showed good growth and developmental conditions, relatively uniform degrees of tree vigor, and no signs of diseases or pests. Flower samples were collected at five developmental stages: floral bud oblate stage (BO), floral bud-yellowing stage (BY), sepal-flattening stage (SF), stamen-wilting stage (SW), and ovary enlargement stage (OE) [23, 24]. We randomly collected whole flowers at the same developmental stage from flower-bearing branches in the middle of jujube trees and placed them into 1.5 ml centrifuge tubes; the tubes were placed into liquid nitrogen when they were full. A total of nine tubes of flowers were collected for each stage. Each stage had three biological replicates, and the mixed samples of every three tubes were considered to represent one biological replicate. Fruit samples were harvested at three developmental stages. At the FrI stage, there were fruits with a normal embryo. In the FrII stage, the jujube fruit embryo began to abort. The samples for the FrII stage were jujube fruits with light embryo abortion. The FrIII stage is the period when the embryo begins mass abortion. The jujube fruit samples for the FrIII stage exhibited a higher embryo abortion degree. We randomly collected whole fruits at the same stage from fruit-bearing branches in the middle of the jujube trees. Each fruit was cut with a razor blade to determine whether it was bearing seeds or not. The ovaries were cut and placed into a 1.5 ml centrifuge tube; tubes were placed into liquid nitrogen when they were full. A total of nine tubes of flowers were collected for each stage. Each stage had three biological replicates, and the mixed samples of every three tubes were considered to represent one replicate. All plant tissues were transferred to and stored in a -80°C ultra-low temperature freezer until analyses were performed.

Microscopic Observations

The routine paraffin section method [25] was used to generate 8-μm sections. Sections were stained using the modified Ehrlich's haematoxylin staining method [26] and mounted with neutral gum. A Leica DMi8 inverted microscope was used to observe and take photographs.

Total RNA Extraction, cDNA Library Construction, and RNA Sequencing

Three biological replicates of each sample at each stage were used for the RNA-Seq experiments. Total RNA was extracted using the CTAB method. The quality and quantity of RNA were assessed by electrophoresis on 1% agarose gels and a Nanodrop 2500 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA) to ensure that the samples were suitable for transcriptome sequencing. RNA integrity was assessed by an Agilent2100 Bioanalyzer (Agilent, Palo Alto, CA, USA) in conjunction with an Agilent RNA 6000 Nano Kit. RNA samples were sent to Wuhan Frasergen Co., Ltd. (Donghu New and High-tech Development Zone, Wuhan, China), where the libraries were produced and sequenced. Oligo (dT) magnetic beads were used to enrich the poly(A) tails of mRNA, and the enriched mRNA was randomly interrupted by divalent cations in the NEB fragmentation buffer. First-strand cDNA was synthesized in an M- MuLV reverse transcriptase system with the fragmented mRNA as the template and random oligonucleotides as the primers. Second-strand cDNA synthesis was subsequently performed using DNA polymerase I and RNase H. The purified double-stranded cDNA was repaired, A-tailed, and connected to the sequencing connector. To select cDNA fragments ~240 bp in length, the library fragments were purified with AMPure XP beads (Beckman Coulter, Beverly, MA, USA). The cDNA fragments were amplified by PCR and purified by AMPure XP beads to obtain the library. Next, a Qubit 2.0 Fluorometer was used for preliminary quantification, and the library was diluted to 1.5 ng/ul. An Agilent 2100 Bioanalyzer was used to detect the insert size of the library. After the insert size met expectations, qRT-PCR was used to accurately quantify the effective concentration of the library to ensure its quality (the effective concentration of the library needed to be higher than 2 nm). Sequencing was performed on an Illumina Hiseq 2500 as per the manufacturer’s protocol.

RNA-Seq Analysis

The software Trimmomatic (v0.4.5) [27] was used to filter the raw reads in the FASTQ format by removing low-quality reads, which consisted of reads in which more than 20% of the bases had a quality lower than 10, reads containing adaptors, or reads with more than 5% of unknown bases. Next, the quality of clean reads was verified using FastQC, including the Q20, Q30, and GC content of the clean data. All downstream analyses were based on clean data of high quality. Clean data were aligned to the Z. jujuba genome (https://www.ncbi.nlm.nih.gov/genome/?term=Zizyphus+jujuba) using Tophat2 [28] and bowtie2 [29], with the mapping parameters “-library-type fr-unstranded -p 6 -G”. The gene expression levels of all samples were estimated by RSEM [30] for each sample, and fragments per kilobase of transcript per million mapped reads (FPKM) were used to estimate gene expression levels.

Differential Gene Expression Analysis

Comparisons were made among seven groups: BO-vs-BY, BY-vs-SF, SF-vs-SW, SW-vs-OE, OE-vs-FrI, FrI-vs-FrII, and FrII-vs-FrIII. The genes of three biological replicates were pooled in each group. DEGs between two groups were identified using the R package DEGseq2 [31] with the following threshold values: |log2 (fold change) | > 1 and FDR < 0.05. Gene Ontology (GO) annotation and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses were performed for the DEGs. Hypergeometric tests of GO analysis and pathway functional enrichments of the DEGs were performed with the Stats R package (version 3.7.0). GO functional enrichment analysis and KEGG pathway functional enrichment analysis were performed via phyper, a function in R. The false discovery rate (FDR) for each p-value was calculated; in general, terms with an FDR < 0.05 were defined as significantly enriched.

WGCNA

WGCNA (weighted gene co-expression network analysis) is a systems biology method for describing the correlations among genes across multiple samples. This method finds clusters (modules) of highly correlated genes and relates modules to external sample traits. Co-expression networks were constructed using the WGCNA (v1.47) package in R [32]. After filtering genes not expressed in more than half of the samples, gene expression values were imported into WGCNA to construct co-expression modules using the automatic network construction function ‘blockwiseModules’ with default settings, except that the power was 15; TOM Type was unsigned, ‘mergeCutHeight’ was 0.25, and ‘minModuleSize’ was 50. Genes were clustered into 12 correlated modules. To identify biologically significant modules, module eigengenes were used to calculate the correlation coefficients with samples. Intramodular connectivity (K.in) and the module correlation degree (MM) of each gene were calculated by the WGCNA R package, and genes with high connectivity tended to be hub genes that might have important functions. For genes in each module, GO and KEGG pathway enrichment analyses were conducted to analyze the biological functions of modules.

Quantitative Real-time PCR (qRT-PCR) Analysis

To validate the accuracy of the RNA-seq results, 15 DEGs were selected for quantitative real-time PCR (qRT-PCR) analysis, and ZjEF1α (GeneID: 107422504) was used as the housekeeping gene to normalize the expression levels. Gene-specific primers were designed by Primer Premier 6.0 software (http://www.PremierBiosoft.com) (Supplementary file ). The cDNA was synthesized from 500 ng of the total RNA in a 10-μL reaction system by PrimeScript™ RT Master Mix Kit (Takara, Beijing, China). The qRT-PCR reactions were carried out using a TB Green™ Premix Ex Taq ™ II (Tli RNaseH Plus) Kit (TaKaRa, Beijing, China). The 20-μL reaction volume contained 5 μL of TB Green Premix Ex Taq II (Tli RNaseH Plus) (2×), 0.4 μL of ROX reference Dye II (50 ×), 0.8 μL of each primer (10 μM), 2 μL of cDNA, and 6 μL of sterile distilled water. The PCR procedure was 95°C for 30 s, 40 cycles of 95°C for 5 s, and 60°C for 30 s, followed by a dissociation stage. The qRT-PCR reaction was performed on an ABI 7500 Real-time PCR System. The relative expression level of each gene was calculated using the 2-ΔΔCt method. Each reaction was performed with three biological replicates.

RESULTS

Morphological Characteristics of Jujube Seeds at Different Developmental Stages

Studies of early embryology have identified ovule abortion to be one of the main causes of seed abortion in Z. jujuba Mill. ‘Zhongqiusucui’ [22]. The characteristic signs of embryo abortion in jujube were the cessation of development and the gradual shrinkage of kernels or seeds, leaving only seed traces (Fig. ). Normally well-developed jujube seeds grew as the fruit developed and filled each ventricle of the ovary (Fig. ). In contrast, abortive seeds ceased growing in the middle, despite the continued development of the fruit. Abortive seeds then gradually withered and shriveled as the jujube fruit matured. We hypothesized that seed abortion in Z. jujuba Mill. ‘Zhongqiusucui’ was caused by ovule abortion; this leads to the prediction that crucial genes mediating the process of abortion should be misexpressed. Therefore, we sequenced the transcriptomes of jujube flowers (BO, BY, SF, SW, and OE; Fig. ) and fruits (FrI, FrII, and FrIII; Fig. ) at different developmental stages to explore the relationship between the expression of genes during ovule development and embryo or seed abortion.

RNA-seq, Reads Mapping, and DEG Analyses

To reveal the molecular mechanisms underlying seed abortion in jujube, a total of 24 samples (including three biological replicates) were sequenced using the Illumina Hiseq 2500 platform. After quality control, 81.12 Gb of high-quality clean reads were obtained. The number of clean reads was between 21,441,794 and 30,551,956. Q20 ranged from 97.4% to 97.6%, Q30 from 93.5% to 93.9%, and GC content from 44.1% to 44.9%. Statistical results are shown in Supplementary file . Clean data were aligned to the Z. jujuba genome (https://www.ncbi.nlm.nih.gov/genome/?term=Zizyphus). The mapping ratio of all clean reads in the sample to the reference genome was greater than 75%, indicating that the 24 samples had a high degree of matching with the Z. jujuba genome (https://www.ncbi.nlm.nih.gov/genome/?term=Zizyphus) (Supplementary file ). The high mapping ratio also demonstrated the reliability of the sequencing results. Thus, transcriptome sequencing data in this study were accurate, and the sequencing quality was high.

We calculated the FPKM of each sample to explore differences in gene expression among different jujube tissues. To obtain the DEGs, seven comparison groups were analyzed, including BO-vs-BY, BY-vs-SF, SF-vs-SW, SW-vs-OE, OE-vs-FrI, FrI-vs-FrII, and FrII-vs-FrIII. We used ‘FDR < 0.05 & |log2FC| > 1’ as the criterion for assessing significant differences in gene expression. When log2FC > 1, the DEG was considered up-regulated. In contrast, when log2FC < -1, the DEG was considered down-regulated. We conducted up- and down-regulated analyses of DEGs for the seven comparison groups, and the results are shown in Fig. (-). Overall, 5,696, 4,722, 7,651, 1,972, 1,254, 5,422, and 762 DEGs were detected for BO-vs-BY, BY-vs-SF, SF-vs-SW, SW-vs-OE, OE-vs-FrI, FrI-vs-FrII, and FrII-vs-FrIII, respectively. We found that SF-vs-SW had the largest number of DEGs, indicating that this is a critical developmental event during floral development. There were overlapped DEGs among the seven comparisons. After removing duplicate overlapping genes, a total of 14,012 DEGs were obtained.

Functional Characterization of DEGs by GO and KEGG Enrichment Analyses

To identify the molecular functions of DEGs at different developmental stages in jujube flowers and fruits, the DEGs in seven comparison groups were analyzed by GO enrichment. The top 50 GO terms enriched in each comparison are listed in Fig. (), Supplementary file and . We analyzed the GO terms enriched in BO-vs-BY and FrI-vs-FRII in detail. According to the classification of GO terms in biological process, most DEGs were enriched in ‘metabolic process (GO:0008152)’, ‘cellular process (GO:0009987)’, and ‘single-organism process (GO:0044699)’ (Figs. and ). These biological processes were likely related to flower and seed development. According to the classification of GO terms in molecular functions, most DEGs were enriched in ‘catalytic activity (GO:0003824)’, ‘binding (GO:0005488)’, ‘transporter activity (GO:0005215)’, and ‘nucleic acid binding transcription factor activity (GO:0001071)’. Among cellular component GO terms, DEGs were primarily enriched in ‘membrane (GO:0016020)’, ‘plasma membrane part (GO:0044459)’, and ‘membrane part (GO:0044425)’ (Figs. and ). All of the above GO terms were closely related to the regulation of biological processes and the synthesis and metabolism of functional substances. The numbers of DEGs enriched in these GO terms in each comparison were large, and both up-regulated and down-regulated DEGs were contained in these GO terms (Supplementary file ).

KEGG pathway enrichment analysis was used to identify major biochemical pathways in which DEGs participated. KEGG enrichment analysis was carried out on DEGs of the seven comparison groups. The top 20 significantly enriched pathways in each comparison group are shown in Figs. ( and ), (Supplementary file and ). The primary significantly enriched pathways (p < 0.05) in the seven control groups were ‘starch and sucrose metabolism (ko00500)’, ‘phenylpropanoid biosynthesis (ko00940)’, ‘flavonoid biosynthesis (ko00941)’, ‘carotenoid biosynthesis (ko00906)’, ‘plant hormone signal transduction (ko04075)’, ‘MAPK signaling pathway plant (ko04016)’, ‘nitrogen metabolism (ko00910)’, ‘protein processing in the endoplasmic reticulum (ko04141)’, and ‘ABC transporters (ko02010)’ (Fig. ). These pathways were significantly enriched in at least three comparisons. Among them, ‘phenylpropanoid biosynthesis (ko00940)’ and ‘flavonoid biosynthesis (ko00941)’ were significantly enriched (p < 0.05) in all seven groups. Both up-regulated and down-regulated DEGs were detected in the pathway ‘phenylpropanoid biosynthesis (ko00940)’ in BO-vs-BY, BY-vs-SF, SF-vs-SW, SW-vs-OE, and FrI-vs-FrII (Supplementary file ). The pathway ‘phenylpropanoid biosynthesis (ko00940)’ was down-regulated in OE-vs-FrI and up-regulated in FrII-vs-FrIII (Supplementary file ). The pathway ‘flavonoid biosynthesis (ko00941)’ was up-regulated in BY-vs-SF and OE-vs-FrI and down-regulated in BO-vs-BY, SW-vs-OE, and FrI-vs-FrII (Supplementary file ). There were both up-regulated and down-regulated DEGs in the pathway ‘flavonoid biosynthesis (ko00941)’ in SF-vs-SW and FrII-vs-FrIII (Supplementary file ). These enriched pathways were associated with a variety of metabolic biological processes, including the synthesis and degradation of hormones, amino acids, sugars, lipids, and secondary metabolites, and were closely related to ovule development.

Gene Co-expression Network Analysis

Gene co-expression networks were constructed based on pairwise correlations between genes using WGCNA [33]. Before WGCNA, we filtered the selected genes to remove low-quality genes or samples that had an unstable impact on the results to improve the accuracy of network construction. A total of 13,155 genes were filtered, and there were 20,037 effective genes after filtering. Modules were defined as clusters of highly interconnected genes; genes within the same cluster had high correlation coefficients. A total of 12 distinct modules (labeled by different colors) were identified and illustrated by a dendrogram (Fig. ). The size of these modules depended on the number of genes that they contained (Fig. ). Among them, the number of genes clustered in the dark green module was the highest (9,394); the coral 2 module contained the lowest number of genes (130). The grey module contained a set of genes that could not be classified into any other modules (Fig. ).

The module eigengene is the first principal component of a given module and can be considered to represent the module’s gene expression profile. A heat map of the expression patterns of samples was performed based on the 12 module eigengenes for the 12 distinct modules (Fig. ). We also performed a gene expression trend analysis for each module separately (Supplementary file ). GO and KEGG pathway enrichment analyses on the genes of these modules were performed (Supplementary files and ). These 12 modules were screened based on the expression patterns and the results of GO and KEGG pathway enrichment analyses of each module.

Dark green and yellow-green modules were composed of genes that were highly expressed in the FrII and FrIII seed abortion stages. The genes clustered in dark orange 2, medium orchid, and saddle brown modules were highly expressed during the key period of seed abortion (FrI). Indiana red 4, plum, and light pink 4 modules contained genes that were highly expressed during the critical stages of the embryo sac (BY and SF) but were expressed at low levels during seed abortion stages (FrII and FrIII). The expression patterns of these seven modules had periodic correlations, with high or low expression levels during the critical period of ovule or seed development (Supplementary file ). The dark green module contained the largest number of DEGs and was significantly enriched for a larger number of biological process GO terms (Supplementary file ). Many biological processes related to development (14 terms), flower development (7 terms), embryo development (7 terms), and reproductive process (5 terms) were significantly enriched, such as ‘developmental cell growth (GO:0048588)’, ‘flower development (GO:0009908)’, ‘embryo development (GO:0009790)’, and ‘embryo sac development (GO:0009553)’ (Supplementary file ). GO terms enriched in the yellow-green module were primarily related to the synthesis, metabolism, and transportation process of cell biosynthesis (10 terms), proteins (28 terms), saccharides (10 terms), and vitamins (6 terms), as well as transport (10 terms) (Supplementary file ). In addition, GO terms related to cell wall formation were significantly enriched, such as ‘cell wall polysaccharide metabolic process (GO:0010383)’ and ‘cell wall organization or biogenesis (GO:0071554)’ (Supplementary file ). The above analysis indicated that key modules were involved in multiple biological processes, including flower development, embryo and seed development, and metabolism, that might be essential for the regulation of seed abortion.

Expression Profiling of Key Genes Related to Ovule and Seed Development

The seed abortion of jujube occurs during the process of ovule and seed development. The regulatory pathway of ovule and seed development has been clarified in the model plant A. thaliana. A series of important genes affecting either ovule or seed development have been identified in A. thaliana [34], providing a valuable reference for similar studies on other plant species. Based on previously collected data [34-37], our study used the ovule and seed development in A. thaliana as a reference, including related genes as precursors, and identified homologous genes in jujube based on sequence similarity searches (BioEdit version 7.0.9.0) (Supplementary file ). A total of 38 homology genes in jujube were selected (Supplementary file ). There were 7 modules (dark green, plum, indian red 4, medium orchid, light green, dark orange 2, and saddle brown modules) for all of the homologous genes; most of them were in the module, and the pink module had the second highest number of homologous genes Supplementary file ). We then used transcriptome data to analyze the expression patterns of these genes during ovule and seed development in jujube. The expression profiles of genes related to jujube ovule and seed development in different functional categories were evaluated, including ovule development in the early stage, seed coat development, and endosperm development (Fig. ).

Many ovule identity genes have been identified in A. thaliana, such as SEEDSTICK(STK)/AG11, SHATTERPROOF1 (SHP1), SHP2, AGAMOUS (AG), BEL1, and SEPALLATA(SEP) [38-40], which regulate ovule development primarily through formatting complexes [34, 41]. The expression levels of ovule identity genes (ZjAG, ZjSEP1/2, ZjSEP3, ZjBEL1, and ZjSHP1/2) were higher during the jujube flower stages (BO, BY, SF, SW, and OE) than during the seed development stages (FrI, FrII, and FrIII) (Fig. ). However, the lowest transcript levels of the STK/AGL11 homolog were detected in the BO stage, and the highest levels were observed in the FrI stage (the early developmental stage of ovules); the expression of the STK/AGL11 homolog was low in the FrII and FrIII stages (Fig. ). These results suggest that STK/AGL11 plays an important role in the regulation of ovule development in jujube.

Potential Genes Related to Ovule and Seed Abortion in Jujube

The seed abortion process has been extensively studied in various plant species, and the number of potential related genes is vast. The MADS-box gene is one of the most important transcription regulators in plants and plays an important role in regulating growth, development, and signal transduction [42, 43]. As previously discussed, some ovule identity genes, such as STK/AGL11, SHP1/AGL1, and SHP2/AGL5 identified in A. thaliana, belong to the MADS-box family of genes [39]. The AG family consists of class-D MADS-box genes that play important roles in ovule development [44]. The ABCG subfamily represents the largest group of ABC transporters, which are widely distributed in various plants and play an important role in growth and development [45]. Studies of A. thaliana, Oryza sativa L., V. vinifera, and other plants have shown that ABCG subfamily members are primarily involved in the formation of the vascular system, lipid transport, pollen wall development, and epidermal cuticle formation [46-49], which may be related to ovule or seed development. Plant metacaspases (MCs) are a family of cysteine proteases that are structurally related to caspase enzymes but differ in the specific sites of their substrates [50]. MC family genes have been shown to play a critical role in the process of programmed cell death (PCD) related to embryo development [51-53]. The abortion of jujube ovules or seeds may also be caused by a similar PCD process in normal embryo development. Thus, we investigated these three types of genes in jujube, including the AG subfamily of MADS-box genes, ABCG transporters, and MC enzymes. Eight AG subfamily genes that may affect ovule or seed development were identified. We found that two AGL61-like, one AGL29-like, one AGL18-like, and one AGL6-like transcripts displayed similar expression patterns with peak expression in the BO sample (Fig. ); the expression of the AGL29-like transcript peaked in the OE sample (Fig. ). The expression of the AGL23-like gene was down-regulated during the developmental stage of the embryo sac (from SF to SW) and then was significantly up-regulated during the seed abortion stages (FrII and FrIII) (Fig. ). Based on sequence similarity analysis, we uncovered 22 ABCG-like transcripts, and five groups with distinctive expression profiles were identified (Fig. ). Four genes, including GeneID:107412473, GeneID:107411471, GeneID:107414661, and GeneID:107412484, were highly expressed during the embryo development stages (FrII and FrIII), suggesting that they play critical regulatory roles (Fig. ). We identified seven MC genes in jujube (Fig. ), and expression analysis revealed that five genes were up-regulated during the late stage of embryo development (Fig. ).

To confirm the accuracy of our RNA-seq data, 15 genes were randomly selected, and their expression patterns were validated by qRT-PCR. The results showed that the expression patterns of these genes were essentially consistent with the expression patterns observed in the RNA-seq data, indicating that our RNA-seq data were highly reliable and reproducible (Fig. , Supplementary file ). Further mining and functional characterization of potential causal genes are needed to elucidate the mechanism underlying seed abortion in jujube.

Fig. (9)

qRT-PCR verification results of 15 genes. The left axis indicates the relative expression values from qRT-PCR experiments; the right axis indicates the RNA-seq expression values in FPKM; r represents the correlation coefficient between qRT-PCR and RNA-seq expression data. (A higher resolution / colour version of this figure is available in the electronic copy of the article).

DISCUSSION

Genes Related to Ovule and Seed Development are Essential Regulators of Seed Abortion

To explore the potential mechanisms underlying ovule and seed abortion in jujube, we exploited the potential regulatory genes based on three aspects: ovule development during the early stages, seed coat development, and endosperm development (Fig. ). The ovule identity genes (STK/AGL11, SHP1, SHP2, and SEP) primarily regulate ovule development through the formation of a complex in A. thaliana [34, 41, 54]. BEL1 interacts with the AG-SEP complex to form the BEL1-AG-SEP complex, which plays an important regulatory role [38]. Similarly, the ovule has been shown to develop abnormally and change into carpellate structures in stk, shp1, shp2, and bel1 mutants [38, 39]. In this study, we found that during the abortion of ovules (FrI-vs-FrII), homologs of STK/AGL11, INO, and BEL1 showed significantly down-regulated expression patterns in jujube (Fig. ). These results suggested that the functions of the ovule identity genes were conserved, and suppression of their expression levels might be a molecular signature of embryo abortion during the early developmental stages.

The down-regulation of VvAGL11 expression in Vitis vinifera L. fruits can promote the formation of seedless fruit [55, 56]. Moreover, in Punica granatum L., the expression level of PgAGL11 initially increased and then decreased during the development of fertile flowers but showed the opposite pattern in abortive flowers [57]. Transcript levels of ZjSTK/AGL11 were low during the seed abortion stages (FrII and FrIII) and fluctuated during the flower developmental stages (Fig. ). Therefore, ZjAGL11 likely played a role in the regulation of ovule and seed development and might be involved in jujube seed abortion. The up-regulated or down-regulated expressions of these ovule identity genes during the critical periods of ovule and seed abortion suggested that these genes might be related to ovule and seed abortion in jujube.

ABCG Transporters are involved in Seed Abortion

Plant ATP-binding cassette subfamily G (ABCG) members are membrane-localized transporters that play an important role in transporting various molecules [58]. Members of ABCG are primarily involved in the formation of the vascular system, lipid transport, pollen wall development, and epidermal cuticle formation [46-49]. In A. thaliana, AtABCG11 and AtABCG12 have been shown to be involved in the formation of the epidermal cuticle [49, 59-62], specifically the transportation of epidermal lipid precursors from epidermal cells to the surface of plant organs [49, 59-62]. Although mutations of ABCGs often lead to male fertility and pollen defects [63], ABCG genes have been found to be involved in seed abortion. VvABCG20 has been shown to be specifically expressed in ovules of V. vinifera L. and participate in the ovule abortion process in seedless grapes; its expression level was significantly higher in seeded grapes than in the seedless variety [64]. When the LeABCG20 gene was silenced in Lycopersicum esculentum, the number of seeds decreased significantly [64]. We previously have shown that microsporogenesis can proceed normally and that embryo abortion occurred during the ovule development stages [22]. In this study, the putative ZjABCG20 (GeneID:107435601) showed up-regulated expression during ‘FrI-vs-FrII’ and had relatively higher transcript levels during the seed abortion stages (FrII and FrIII) (Fig. ). A putative ZjABCG26 (GeneID:107416341) had an expression pattern similar to the putative ZjABCG20 (GeneID:107435601) (Fig. ). These results indicated that these members could be important regulators for seed development. However, additional studies are needed to establish the biochemical functions of ZjABCG26 and ZjABCG20, including their roles in ovule and seed abortion in jujube.

Expression Profiles of MCs are a Molecular Signature of PCD and Seed Abortion

PCD is a major mechanism of seed abortion in various plants, which is regulated by various internal and external factors [65]. For example, in A. thaliana, salt stress can induce PCD at multiple stages of gametophyte development and cause seed abortion [66]. The MC family genes play an important role in the process of PCD, and MC-dependent PCD is one of the most important processes in plant embryo development [51-53]. The abundance and enzymatic activity of MCs are one of the most characteristic molecular signs of PCD [51]. In a previous study in grape, Zhang et al. [67] found that the MC family genes VvMC1, VvMC3, and VvMC4 were differentially expressed during ovule development, and their expression patterns coincided with the endosperm abortion period of seedless grape. In this study, we identified seven MC-like transcripts that were highly expressed through transcriptomic analysis (Fig. ). We found that five MC-like genes displayed similar expression patterns: expression was low during early stages of ovule development (SW, OE, and FrI) and significantly up-regulated during ‘FrI-vs-FrII’; expression levels of these genes were highest during the seed abortion stages (FrII and FrIII) (Fig. ). There were multiple annotated copies of MC genes corresponding to the Arabidopsis reference (Fig. ), which indicated that the MC gene family might be expanded in jujube. The differential expression patterns of copies suggested that the functions of MCs were diversified during flower and seed development. Overall, we speculate that some MC genes regulate seed abortion by mediating the PCD pathway.

CONCLUSION

Here, we conducted an integrative transcriptomic analysis of ovule and seed development in jujube to explore the mechanisms underlying ovule and seed abortion. RNA-seq data were used to identify DEGs as well as GO terms and pathways during different stages of ovule and seed development. Our study revealed a suite of potential key genes related to ovule and seed development, as well as their regulatory patterns, and mapped the networks of potential genes related to jujube ovule and seed development in three aspects: ovule development during the early stages, seed coat development, and endosperm development. Our data will aid future efforts to explore the mechanisms underlying ovule and seed abortion. We also screened related genes from the AG subfamily as well as ABCG and MC family genes. Seed development is a complex process with many layers of regulation. This RNA-Seq dataset promises to provide insight into the molecular mechanisms affecting ovule and seed abortion, contribute to future genetic studies, and enhance the quality of breeding in jujube.