Analyzing Differentially Expressed Genes and Pathways Associated with Pistil Abortion in Japanese Apricot via RNA-Seq

Reproduction is a critical stage in the flower development process, and its failure causes serious problems affecting fruit quality and yield. Pistil abortion is one of the main factors in unsuccessful reproduction and occurs in many fruit plants. In Japanese apricot, the problem of pistil abortion is very common and affects fruit quality and plant yield; however, its molecular mechanism is not clearly understood. Therefore, in the current study, we used RNA-Seq to identify the differentially expressed genes (DEGs) and pathways actively involved in pistil abortion. A total of 3882 differentially expressed genes were found after cutoff and pairwise comparison analysis. According to KEGG pathway analysis, plant hormone signaling transduction and metabolic pathways were found most significantly enriched in this study. A total of 60 transcription factor families such as MADS-box, NAC and TCP showed their role in this process. RT-qPCR assays confirmed that the expression levels were consistent with RNA-Seq results. This study provides an alternative to be considered for further studies and understanding of pistil abortion processes in Japanese apricot, and it provides a reference related to this issue for other deciduous fruit crops.


Introduction
Japanese apricot (Prunus mume Sieb. et Zucc) is an important fruit and ornamental plant. It has great economic importance, high profit and world market demand, and it is used in value-added products like jams, pickles, etc. [1,2]. Flower development is a fundamental stage in the life cycle of the plant and plays a significant role in the sexual reproduction process [3]. Pistil abortion is a widespread phenomenon occurring in fruit plants and has been discussed in different fruit crops like pomegranate, Xanthoceras sorbifolia and olive [4][5][6]. In Japanese apricot, the problem of pistil abortion is very common

RNA-Seq Data Analysis
After sequencing, the expression level of each sequence library was standardized as FPKM (Fragments Per Kilobase of Transcript per Million), and the most differentially expressed genes (DEGs) were selected for further analysis. Hierarchical clustering was performed using TBtools (v 0.667). The significance level of the transcripts was measured using the FDR (False Discovery Rate) control method [27] to rationalize p-values. Gene regulation was determined through cutoff values with an absolute log 2 fold change of ≥+1 (up-regulated) and ≤−1 (down-regulated) method with p-values less than 0.001.
To find out the putative transcription factor (TF) related to these DEGs, Getorf (EMBOSS: 6.5.7.0) was used to identify ORF (Open Reading Frame), hmmsearch (v 3.0) was used to align ORF to the TF protein domain, and then the unigenes were identified according to TF family characterized by PlantTFDB 4.0. To further characterize the function of DEGs, GO (Gene Ontology) and KEGG (Kyoto Encyclopedia of Genes and Genomes) enrichment were performed. Based on DEGs, GO is functionally categorized as follows: (1) biological process (BP), (2) cellular component (CC) and (3) molecular function (MF). On the basis of differential gene expression analysis, the biological functions of differentially expressed genes were further understood by using GO (Gene Ontology). GO functional enrichment analysis was carried out according to their functional classes using FDR ≤0.05 as the significant enrichment level. For pathway analysis, all identified DEGs were mapped to the KEGG database. Enrichment analysis of both GO and KEGG was performed using the phyper function in R software with p-value < 0.05 [28]. The data from sequencing were deposited to SRA-NCBI under accession PRJNA646597.

Endogenous Hormone Measurement
The phytohormone contents were determined using the HPLC method described earlier [29]. The samples identical to RNA-Seq were used for hormone measurements. Three biological replications were used for each sample.

Quantitative Reverse Transcription PCR for Data Validation
Twelve different genes were randomly selected to verify the RNA-Seq results. Gene sequences were extracted, and their primers were designed using Primer 3 software. RNA extraction and quality tests were performed as described above, and RT-qPCR was performed according to the method described by [30,31]. Using RPII as a reference internal gene, the expression levels of the genes were analyzed via the 2 − CT method [32].

An Overview of RNA-Seq Libraries
In the current study, six different cDNA libraries were constructed for normal and abortive pistils and sequenced using Illumina sequencing, producing an average of 65.8 M raw reads per sample. After the removal of raw reads and adopter sequences, 63.27 M clean reads per sample were obtained and mapped to the P. mume genome. After mapping to the reference genome, a total of 86.64% mapped reads were obtained from each sample, while uniquely mapped reads were 63.94% per sample (Table 1).

DEG Analysis of Normal and Abortive Pistils
Hierarchical clustering was performed for differentially expressed genes on the basis of their fold change values (log 2 fc), and the clustering is shown in Figure 1A. Pairwise comparison was used to analyze the regulation (up/down) and gene expression profiles between AP and NP. A total of 3882 unigenes (2033 up-regulated and 1849 down-regulated) were identified from normal and abortive pistil comparisons ( Figure 1B). The volcano plot shows the expression level and overall distribution of differentially expressed genes between these two groups ( Figure 1C). All identified DEGs and their expression values are listed in File S1.

DEG Functional Enrichment Analysis
The DEGs in NP vs AP comparisons and their functional classes were collected, and the most significant enriched terms are shown in Figure 2. In the BP category, cellular process, metabolic process and biological regulation were the leading terms; for the CC category, the most enriched terms were membrane, membrane part and cell; while in the MF category, binding and catalytic activity were the prominent terms.  For KEGG pathway study, all identified DEGs were mapped to the KEGG database, and further involvement of genes in different pathways was predicted. A total of 134 KEGG pathways were identified based on enrichment analysis (File S2). Plant hormone signal transduction, plant-pathogen interaction and different metabolic pathways were found most enriched related to this study. Gene regulation (up and down) and different processes involved in the pathway study are shown in

Expression Profiling of Transcription Factor Encoding Genes
Transcription factors (TFs) are the key regulatory proteins involved in different biological and physiological processes, especially in flower organogenesis and floral development [33,34]. Therefore, in the present study, we mapped all DEGs to the plant transcription factor database (PlantTFDB) to study the expression pattern of key transcription factors. A total of 1470 TF genes were identified and divided into different transcription factor families. Among all TF families, MYB, AP2-EREBP, bHLH, NAC, MADS WRKY and TCP were the most relevant to this study. Figure 4A shows the distribution and percentage of genes in each transcription factor family, and Figure 4B

DEGs in Response to Plant Hormone Signaling Transduction Pathway
Earlier studies reported that various phytohormones such as auxin, cytokinin (CK), gibberellin (GA), etc., participate in different floral development processes [35]. Therefore, the regulation and expression pattern of the DEGs related to plant hormones were analyzed. An overview of the plant hormone signaling transduction pathway is shown in Figure 5 and their values are listed in File S3. In total, 59 genes in the auxin-signaling pathway, 23 genes in the cytokinin-signaling pathway, 36 genes in the GA-signaling pathway, 21 genes in the ABA-signaling pathway, 9 genes in the Ethylene-signaling pathway, 36 genes in the BR-signaling pathway, 28 genes in the JA-signaling pathway and 12 genes in the SA-signaling pathway were identified. Gene expressions related to different plant hormones involved in this study are shown in Figure 5A The endogenous contents of these hormones were determined in both normal and abortive pistils ( Figure 5B). The contents of auxin, cytokinin and ethylene were significantly higher in normal pistils-50.03, 424.43 and 25.87 ng/g, respectively-than those of abortive pistils. The contents of gibberellin and abscisic acid were significantly higher in abortive pistils, 234.05 and 130.28 ng/g, respectively, than those of abortive pistils. Lower contents of auxin, cytokinin and ethylene-28.29, 332.78 and 19.73 ng/g-were found in abortive pistils, while lower contents of gibberellin and abscisic acid, recorded as 172.87 and 69.81 ng/g, respectively, were found in normal pistils.

DEGs in Response to Different Metabolic Pathways
In the present study, different metabolic pathways, such as phenylpropanoid biosynthesis (268 genes), diterpenoid biosynthesis (38 genes), flavonoid biosynthesis (90 genes), isoflavonoid biosynthesis (24 genes) and zeatin biosynthesis (15 genes), were identified. The expression patterns of DEGs involved in different biosynthesis-related pathways of normal and abortive Japanese apricot pistils were analyzed ( Figure S2 and File S4).

Validation of Genes Through RT-qPCR
To verify the results of RNA-Seq, twelve genes were randomly selected, and their relative expression level was verified through RT-qPCR. A detailed description of the selected genes and primers is listed in File S5. The expression levels of all genes were consistent with the results of RNA-Seq ( Figure 6).

Discussion
Sequential regulation of genes plays a significant role in plant growth and development, and it is important to have detailed information on a gene to understand the molecular mechanism underlying any developmental process [36]. In Japanese apricot, the problem of pistil abortion also occurs. The molecular mechanism of pistil abortion in Japanese apricot is not clear. Therefore, this study aimed to analyze the regulation of genes involved in this process. Numerous studies have already reported on genome-wide transcriptional programs during reproductive development [37,38]. In the present study, RNA-Seq was performed using Illumina sequencing, and a total of 3882 genes were generated after cutoff. Cluster analysis of the identified genes was performed to uncover the expression pattern and changes during pistil abortion. Furthermore, a large number of TF families were also recognized, and some of them were strongly related to the stated problem. In general, our study provides an alternative way to further study and understand the expression profiles of molecular mechanisms and their gene regulation during pistil abortion in Japanese apricot.
TFs play an important role in flower development and regulation processes [39]. In this study, we found a large number of TF families, and most of them were closely related to flower development. From these TF families, MYB, AP2-EREBP, bHLH, NAC, MADS, WRKY and bZIP were highly expressed and have been previously reported in various plant species during the reproductive process [37]. In Arabidopsis thaliana, the bHLH TF family is involved in various flower development processes [40], while bZIP has also been identified in many plant species as the key regulator in flower formation and development [41]. Some other TFs such as YABBY and GRF have also been reported in many studies to have an important role in the reproductive process, and studies have shown that decreasing the expression level of GRF TFs has an effect on normal pistil development. In our study, we also noted a decreased level of GRF TFs, confirming a previous study on Arabidopsis [42] that this TF might have a specific role in pistil abortion in Japanese apricot.
Plant hormones are signaling molecules that play an important role in plant physiological processes such as plant metabolism, morphogenesis and growth [43], and numerous studies have shown that these hormones regulate the pistil development process [44]. In the current study, a total of 131 DEGs related to plant hormones were identified, and their expression patterns were analyzed. Auxin is an important plant hormone regulating plant growth and development [45], and the DEGs involved in this pathway were associated with AUX1, TIR1, ARF, GH3 and SAUR sub pathways. In Arabidopsis, ARF6 and ARF8 are highly expressed in the sepals, petals and stamen and promote floral transition [46]. In this study, most of the DEGs related to AUX/IAA, GH3, SAUR, TIR1 and ARF were up-regulated and considered to promote cell expansion, plant growth and development, which is consistent with a previous study [47]. At the same time, the content of auxin was significantly higher in normal pistils as compared to abortive pistils. Thus, auxin and its transport may be required for proper pistil development. In pistils, cytokinin is an important hormone, and an increase in CK content can increase pistil morphology and ovule numbers [48]. In this study, the content of cytokinin was significantly higher in normal pistils than abortive pistils. In the CK metabolic pathway, the genes related to A-ARR were down-regulated in abortive pistils, which might be a positive regulatory response to decrease the CK content in pistils, and are considered one of the most important causes for pistil abortion. In addition, GA 3 plays an important role in flower development [49], regulates sex differentiation in plants and inhibits pistil development at an appropriate level [49]. Our results showed that the content of gibberellin was significantly higher in abortive pistils than that of normal pistils. Earlier reports in Arabidopsis showed that gal-3, a GA-deficient mutant, restricted floral organ growth, and its phenotype could be rescued through exogenous application of GA [50]. In the current study, GID1 (GA metabolic pathway) related genes were highly expressed, and it might be said that it reduces the biological activity of GA.
ABA is also an important hormone due to its significant role in plant growth and development, including flower development, seed development and environmental stress responses [51]. In Arabidopsis, ABA promotes flower bud formation through regulating flowering and photoperiod responsive genes [52], while in apples, ABA inhibits flower formation [53]. In our study, genes related to the ABA metabolic pathway showed significant expression changes, and genes related to sub pathways (ABF, PP2C and SnRK2) reported an active role during abiotic stress conditions and could adopt to any adverse environment [54,55]. In this study, the content of ABA was significantly higher in abortive pistils than in normal ones, and the genes related to ABA also showed higher expression in abortive pistils as compared to normal pistils. Ethylene is also an important hormone for floral organ development, especially pistil development, and in tobacco ethylene is an upstream regulator for ovule development [56]. In this study, the content of ethylene was significantly higher in normal pistils as compared to abortive pistils. The expression level of ethylene-related genes was also higher in normal pistils, indicating that ethylene might be a key regulator involved in pistil abortion in Japanese apricot.
Moreover, genes regulating different metabolic pathways, including phenylpropanoid biosynthesis, flavonoid and isoflavonoid biosynthesis, zeatin biosynthesis and their intermediates, had lower expression levels in abortive pistils, indicating that they were involved in the developmental process. In rice, most metabolic pathways are found enriched during floral organ development [37], as with many other plant species [57,58]. Flavonoid synthesis is required for pollen function and is very common in several floral organs [59].

Conclusions
In the current study, RNA-Seq was used to determine the regulation and expression pattern of genes in normal and abortive pistils of Japanese apricot. A total of 3882 DEGs were identified, and a comparison analysis was performed between normal and abortive pistils. We mainly focused on the expression pattern/trends of genes related to plant hormones and metabolic processes. Changes in the activity of hormones such as auxin, CK, GA and ethylene play an important role in pistil abortion. Related genes are involved in hormone synthesis and expression to regulate hormone content and promote pistil abortion under adverse conditions. Our study provides a basis for understanding the molecular mechanism of pistil abortion, and it provides a reference for future studies related to this issue in other fruit crops.
Supplementary Materials: The following are available online at http://www.mdpi.com/2073-4425/11/9/1079/s1, Figure S1: An overview of the present study (A) The overall experimental design process (B) The route map and analysis for RNA-Seq data, Figure S2: Expression analysis of the DEGs related to different metabolic pathways involved in this study, File S1: Total number of differentially expressed genes in this study, File S2: Significantly enriched pathways and their DEGs identified in this study, File S3: Number of DEGs and their expression involved in the plant hormone signaling transduction pathway, File S4: Number of DEGs and their expression involved in different metabolic pathways, File S5: Primers used for RT-qPCR.