New Insights of Salicylic Acid Into Stamen Abortion of Female Flowers in Tung Tree (Vernicia fordii)

Tung tree (Vernicia fordii), an economically important woody oil plant, is a monoecious and diclinous species with male and female flowers on the same inflorescence. The extremely low proportion of female flowers leads to low fruit yield in tung orchards. The female flower normally develops along with stamen abortion; otherwise sterile ovules will be produced. However, little knowledge is known about the molecular basis of the female flower development in tung tree. In this study, integrated analyses of morphological and cytological observations, endogenous phytohormone assay and RNA-seq were conducted to understand the molecular mechanism of the female flower development in tung tree. Cytological observation suggested that the abortion of stamens in female flowers (SFFs) belongs to the type of programmed cell death (PCD), which was caused by tapetum degeneration at microspore mother cell stage. A total of 1,366 differentially expressed genes (DEGs) were identified in female flowers by RNA-seq analysis, of which 279 (20.42%) DEGs were significantly enriched in phenylpropanoid biosynthesis, phenylalanine metabolism, flavonoid biosynthesis, starch and sucrose metabolism, and plant hormone signal transduction. Stage-specific transcript identification detected dynamically expressed genes of important transcription regulators in female flowers that may be involved in PCD and floral organ development. Gene expression patterns revealed that 17 anther and pollen development genes and 37 PCD-related genes might be involved in the abortion of SFF. Further analyses of phytohormone levels and co-expression networks suggested that salicylic acid (SA) accumulation could trigger PCD and inhibit the development of SFF in tung tree. This study provides new insights into the role of SA in regulating the abortion of SFF to develop normal female flowers.


INTRODUCTION
Flower development attracts great attention as a fascinating topic for studying plant development and evolution. Unisexuality is considered to be an important transition in the evolutionary history of angiosperms (Barrett, 2010). Many flowers become unisexual after floral organs are specified, but during the process of differentiation, carpel or stamen abortion or arrest occurs and organs become non-functional (Sobral et al., 2016). Plant reproduction requires the development of complex structures that interact with each other and that have an inherently limited life span. Thus, programmed cell death (PCD) is involved in shaping the sexual and non-sexual organs of the flower and in their removal once they are no longer needed (Wu and Cheun, 2000). For example, the gynoecium in male flowers is degenerated by PCD just after initiation of floral primordium in Zea mays (Cheng et al., 1983). In Actinidia deliciosa, PCD also induces male sterility in female flowers (Coimbra et al., 2004).
Researchers suggest that PCD has a close relationship with the accumulation of salicylic acid (SA), and SA controls the timing and extent of PCD in the hypersensitive response (HR) (Greenberg, 1997;Alvarez, 2000;Kovacs et al., 2016). As an endogenous signaling molecule, however, SA is suggested in connection with plant flowering. For example, SA induces the expression of FLOWERING LOCUST 2 (FT2) in the regulation of the flowering in Pharbitis nil (Yamada and Takeno, 2014). In Arabidopsis, SA controls early flower development by NON-EXPRESSOR OF PR3 (NPR3) (Shi et al., 2013). Moreover, SA plays an important role in regulating pollen viability and floret fertility in rice (Zhao et al., 2018). SA can also induce mRNA accumulation of Callose synthase 5 (Gsl5) for deposition of callose in pollens (Ostergaard et al., 2002).
Tung tree (Vernicia fordii), native to China, is widely planted in China and other countries for its ornamental purpose and tung oil production (Cui et al., 2018;Li et al., 2018). Tung oil shows strong dryness, strong adhesion, acid and alkali resistance, so it has been widely used in various industrial fields (Cahoon et al., 1999;Chen et al., 2010). In recent years, tung oil has been attracted world-wide attention due to production security, environmental concerns, and negative effect of synthetic chemical coatings on human health (Meininghaus et al., 2000;Tsakas et al., 2011;Wei et al., 2012;Yang et al., 2018). As a monoecious plant species, tung tree produces a low ratio of female to male flowers (approximately 1:27) (McCann, 1942) and functionally abnormal female flowers, which results in low fruit yield in tung orchards. Generally, female tung flowers present bisexual characteristics at early stages, and then switch to unisexual flowers at late stages along with cell death in stamens in female flowers (SFFs) (Mao et al., 2017). However, the reason causing cell death and the transition of tung flowers from bisexuality to unisexuality remains uncertain. Mao et al. (2017) found that PCD occurred during stamen abortion in female tung flowers, whereas the molecular mechanism is unclear. Toward this end, we conducted integrated analyses of morphological and cytological observations, endogenous phytohormone assay and RNA-seq of male and female flowers. Based on these analyses, we proposed a possible mechanism for female flower development in tung tree.
of Forestry and Technology (Qingping Town, Yongshun County, Hunan Province) under natural conditions. The flower was randomly collected every 5 days from March to April, 2017. Male and female flowers at stage 2 in 30 days before flowering (DBF, C1 and X1, Figures 1C,D), stage 4 in 20 DBF (C2 and X2, Figures 1G,H), stage 6 in 10 DBF (C3 and X3, Figures 1K,L) were separately collected from the same plant and used for transcriptome sequencing. Three biological replicates were analyzed for each sample with 100 flowers counted as one replicate. Controlled pollination ('Putaotong' × 'Duiniantong') was performed in April, 2017. A total of 300 female flowers and 300 bisexual flowers were used to evaluate the fruit setting with three replicates (each 100 flowers).

Morphological and Cytological Observation
Morphology of male and female flowers at six developmental stages was observed by a 3D super depth digital microscope (ZEISS Smartzoom 5, Germany). The samples were prefixed and observed with the paraffin sectioning method described by Liao et al. (2014). The samples were sectioned into 8-10 µm sections using a Leica RM 2265 microtome (Leica Camera AG, Germany), and were stained with Safranin O and Fast Green FCF according to the Sass's method (Ruzin, 2000). Observations and photomicroscope of the sections were conducted using a microscope (Olympus BX-51, Japan).

RNA Extraction and Sequencing
Total RNAs were isolated using RNAprep pure plant kit DP441 (TIANGEN, China), according to the manufacturer's protocol. RNA sequencing was performed by Illumina Hiseq 2000 (Illumina, United States). Our reference genome sequences and annotations of tung tree were appliced for the analysis of RNA-seq data. Raw reads of the RNA-seq data are found in NCBI Sequence Read Archive database with the Accession Nos. SRX3843588; SRS3089151; SRS3089154; SRS3089148; SRS3089147; and SRS3089150.

Data Analysis
The detailed processes of de novo assembly and functional annotation were performed as described by Feng et al. (2017) The number of all mapped reads for each gene was counted and normalized into the Fragments Per Kilobase of transcript per Million fragments mapped (FPKM) (Florea et al., 2013). Next, DESeq (Anders and Huber, 2010) was used to identify differentially expressed genes (DEGs), and genes with an adjusted false discovery rate (FDR) ≤ 0.01 and log2 FoldChange absolute value of ≥ 1 were marked as significantly different between the two samples. Stage-specific genes were selected as described by Feng et al. (2017). GO and KEGG pathway analyses were performed, and for each category, a two-tailed Fisher's exact test was employed to test the enrichment of the identified protein against all database proteins. Correction for multiple-hypothesis testing was performed using standard FDR-control methods. The venn diagram, clustering, heatmap, and principal component analysis (PCA) were performed with the venn diagram function, the k-means function, pheatmap, and factoextra package in R software, respectively.

Co-expression Networks Analysis
The transcripts with an average abundance (calculated from three biological replicates) of ≥ 1 (FPKM value) in at least one of the six samples were utilized to construct a coexpression network using the weighted correlation network analysis (WGCNA) R software package (Langfelder and Horvath, 2008). Modules were constructed using the following parameters: maxBlockSize = 10000, power = 18, networkType = "unsigned, " mergeCutHeight = 0.25, minModuleSize = 30. The significant modules in female flowers were determined according to an effective threshold (≥0.85) of the Pearson's correlation coefficient (PCC) value and a p-value of ≤ 0.05. In order to identify genes which were correlated in expression with those genes involved in anther and pollen development, PCD, and SA, the top 50 genes in each significant module were used for network construction according to correlation degree. The co-expression networks were visualized using the Cytoscape v.3.5.1 program 1 (Shannon et al., 2003).

Quantitative Real-Time PCR (qPCR) and Phytohormone Levels Analysis
A PrimeScript RT enzyme with a gDNA eraser (Takara, Japan) was used for cDNA synthesis. QPCR was performed on a Bio-Rad CFX96 Real Time PCR system using SYBR Premix ExTaq II (Takara, Japan). The primers in this step were listed in Supplementary Table S1. Tung Elongation Factor 1-α (EF1α) was used as the internal control (Han et al., 2012). The relative expression levels were calculated using the 2 − CT method (Livak and Schmittgen, 2001). Phytohormone levels were analyzed using the methods described by Pan et al. (2002).

Statistical Analysis
All experiments were performed with three biological repeats. The data were analyzed by one-way ANOVA procedure in SPSS 22.0 (IBM Corporation, United States). All figures showed the average value of three repeats. All data were expressed as means plus or minus standard deviations (mean ± SD).

Tapetum Degeneration at Microspore Mother Cell Stage Induces PCD and Abortion of SFF
Tung tree is a diclinous species and the functionally normal female flower develops along with stamen abortion in tung tree (Mao et al., 2017). However, the development of bisexual flowers often occurs in inflorescence with abnormal ovaries, leading to lower fruiting ratio (∼5%) than normal female flowers (∼68%) (Supplementary Figure S1). Therefore, we focused on staminode development in female flowers aiming to clarify the process of stamen abortion in tung tree. The development of female flowers could be divided into six important stages (stages 1-6) ( Figures 1A-F). Female flowers possessed obvious stamens at stages 1 and 2 which initially degenerated from stage 3 and completely disappeared until stage 6 ( Figures 1A-F), whereas no pistil was observed and only intact stamen developed in male flowers across all the development stages (Figures 1G-L).
The SFF showed delayed development in comparison with stamens in male flowers (SMFs). The SFF developed into sporogenous cell stage at 35 DBF (stage 1) (Figure 1M), while the SMF developed into early microspore mother cell (MMC) stage ( Figure 1S). When the SFF developed into early MMC stage at 30 DBF (stage 2) (Figure 1N), the SMF developed into later MMC stage ( Figure 1T). Interestingly, the SFF stopped growth at early MMC stage with tapetum cells degeneration and large vacuoles in MMC at 25 DBF (stage 3) (Figure 1O), while the SMF developed into tetrads stage ( Figure 1U). From 20 to 10 DBF (stages 4 to 6), the SFF shrank at carpel base, and the MMC was mostly occupied by vacuoles with nucleuses in MMC disappearing (Figures 1P-R), while the SMF developed into binuclear pollen stage (Figures 1V-X).
Taken together, the functionally normal female flower developed along with stamen abortion in tung tree. The SFF developed normally at stages 1 and 2, but degenerated from stages 3 to 6. Our morphological and cytological observation further confirmed that the abortion of SFF in tung tree belongs to the type of PCD. The tapetum cell degeneration should contribute to the abortion of SFF, finally leading to PCD in MMC.

DEGs in Male and Female Flowers
In order to elucidate the molecular mechanism responsible for the abortion of SFF in tung tree, the male and female flowers at stage 2 (X1, C1), stage 4 ( X2, C2), and stage 6 (X3, C3) were RNA-sequenced. Based on our tung tree genome sequence, a total of 22,344 genes were identified, of which 18,334 genes were expressed (FPKM value ≥ 1) across six flower samples.
Further, we identified 3,147 DEGs in male flowers (X2_vs_X1, X3_vs_X2) and 1,366 DEGs in female flowers (C2_vs_C1, C3_vs_C2) (Figure 2A). The PCA analysis of the total DEGs, revealed that X1 and C1 had high correlation ( Figure 2B), which indicating that the SFF in C1 had normal structure and function as MFF in X1. Besides, among the 1,366 DEGs in female flowers, 531 and 257 genes were significantly up-regulated in C2 and C3, respectively, and 205 and 524 genes were significantly down-regulated in C2 and C3, respectively, as compared to their preceding stage (Figure 2A). Furthermore, in C2_vs_C1 group, there were 356 specific DEGs among up-regulated genes, but 112 among down-regulated genes. In C3_vs_C2 group, 137 specific DEGs were up-regulated and 416 specific DEGs were downregulated. These results showed that there were more genes highly expressed in C2 than C1 and C3, indicating that a large number of genes might participate in the transition from bisexual to unisexual flowers formation (Figures 2C,D).
KEGG annotation showed that 496 (15.76%) DEGs in male flowers were enriched in 110 KEGG pathways ( Supplementary  Table S2), of which the starch and sucrose metabolism, phenylpropanoid biosynthesis, pentose and glucuronate interconversions, flavonoid biosynthesis, amino sugar and nucleotide sugar metabolism were significantly enriched pathways ( Figure 2E). In contrast, a total of 279 (20.42%) DEGs in female flowers were enriched in 103 KEGG pathways (Supplementary Table S2), of which the phenylpropanoid biosynthesis, phenylalanine metabolism, flavonoid biosynthesis, starch and sucrose metabolism, plant hormone signal transduction were significantly enriched pathways ( Figure 2F). It has been reported that sucrose accumulation can change hormonal levels in cut lily flowers and phenylpropanoid can provide the key precursor for SA biosynthesis (Mauchmani and Slusarenko, 1996;Arrom and Munne-Bosch, 2012). In this study, 13 DEGs were found to be significantly enriched in plant hormone signal pathway, indicating that the phytohormones may play an important role in regulating the development of female flowers in tung tree.

Identification of Candidate Genes Involved in PCD in the Abortion of SFF
Cytological observation confirmed that the abortion of SFF belongs to PCD type in tung tree. To gain insights into the role of PCD-related genes in the abortion of SFF, the putative 37 key PCD-related genes were identified based on their expression patterns which included 26 positive genes and 11 negative genes ( Figure 4A and Supplementary Table S5). Most of the positive genes exhibited high expression levels in C2 and C3, such as Leucine-rich repeat receptor-like serine/threonine/tyrosine-protein kinase (SOBIR1), AGP16, and BAG family molecular chaperone regulator 7 (BAG7) etc. (Figure 4A and Supplementary  Table S5). These genes have been reported to behavior in PCD (Gao and Showalter, 1999;Li et al., 2017;Sala et al., 2017). Likewise, negative genes showed low expression levels in C2 and C3 and high expression in C1 such as Lesion Simulating Disease 1 (LSD1), Calcium/calmodulin-regulated receptor-like kinase 2 (CRKL2), Homologies of Respiratory Burst oxidase homolog protein A (RBOHA) and so on (Torres et al., 2002;Yang et al., 2010;Guo et al., 2013) (Figure 4A and Supplementary Table S5). The qRT-PCR analysis of six genes revealed consistent expression patterns with those generated by RNA-seq data (Figures 4H-M). The expression patterns demonstrated the important role of PCD-related genes in the abortion of SFF in tung tree.

Phytohormone Levels in Male and Female Flowers
To explore the possible phytohormone regulation in the abortion of SFF, the endogenous levels of auxin, abscisic acid (ABA), gibberellin (GA), jasmonic acid (JA), cytokinin (CK), and SA in tung flowers at different stages were measured ( Table 1). The male flower exhibited the same patterns of IAA, ABA, and GA levels with the female flower. The IAA level was upregulated reaching the highest at stage 3, and it was significantly higher in X3 (8.35 ng/g) than that in C3 (1.68 ng/g) ( Table 1). The ABA levels were down-regulated across three stages, and it was significantly higher in C1 (165.68 ng/g) than that in X1 (87.64 ng/g) ( Table 1). For GAs, GA4 level was significantly higher than GA1, GA3, and GA7 in all samples. In addition, the GA4 level was the lowest in C2 (1.66 ng/g) and X2 (1.0 ng/g) ( Table 1).
The levels of JAs and CKs showed different patterns in comparison with IAA, ABA, and GA (Table 1). JA was significantly higher than MEJA across all samples. It was up-regulated in male flowers with the highest value (129.4 ng/g) in X3 (Table 1). For CKs, the TZR level was significantly higher than Zeatin, IP and IPA. The TZR level was up-regulated in both male and female flowers (Table 1). Notably, SA showed different patterns from other phytohormones. No SA was detected in C1, X1 and X2, and only extremely low levels of SA were detected in C3 (0.88 ng/g) and X3 (0.86 ng/g) ( Table 1). In contrast, SA was significantly high in C2 (16.93 ng/g) where the abortion of SFF occurred.
According to the analysis of phytohormone levels, SA is likely to play an important role in the abortion of SFF instead of other phytohormones.

Analysis of SA Synthesis and Signaling Pathway
To uncover the regulatory mechanism of SA involved in the abortion of SFF, we analyzed in detail the pathway of SA synthesis and signaling in tung flowers. Two pathways have been reported to be involved in plant SA synthesis (Chen et al., 2009). One is regulated by phenylalanine ammonia lyase (PAL) and chorismate mutase (CM), and the other by isochorismate synthase (ICS). In Arabidopsis, there are four homologs encoding PAL (PAL1-4), three homologs encoding CM (CM1-3), and two homologs encoding ICS (ICS1/2) . In tung tree genome, we found three CM members (CM1-3), one PAL and one ICS (ICS2). The three CMs displayed different expression patterns. Generally, FPKM of CM1 was higher than CM2 and CM3 across all samples. CM1 displayed different expression profiles in male and female flowers, and was up-regulated with the development of female flowers ( Figure 5A). CM2 was highly expressed in C2 of female flowers, and down-regulated in male flowers ( Figure 5A and Supplementary Table S6). CM3 was down-regulated in both male and female flowers (Figure 5A and Supplementary Table S6). PAL was up-regulated in both male and female flowers. ICS2 was highly expressed in C2 of female flowers, and down-regulated in male flowers (Figure 5A and Supplementary Table S6).
In Arabidopsis SA signaling pathway, three receptors are reported, including NON-EXPRESSOR OF PR1 (NPR), transcription factor TGA (TGA), PATHOGENESIS-RELATED (PR) (Kaltdorf and Naseem, 2013). Three homologs of NPR (NPR1/3/5), four TGAs (TGA1-3, TGA7), and two PRs (PR-1, PR-4B) were found in tung tree genome (Figure 5A and Supplementary Table S6). The expression level of NPR3 was higher than NPR1 and NPR5, and down-regulated in male and female flowers of tung tree. NPR1 was down-regulated in female flowers, but up-regulated in male flowers. NPR5 was highly expressed in C2 of female flowers, and up-regulated in male flowers ( Figure 5A and Supplementary Table S6). Among the TGA family, TGA3 showed higher expression level than the other three members. It was highly expressed in C2 of female flowers, and down-regulated in male flowers (Figure 5A and Supplementary Table S6). Interestingly, two members of the PR family (PR-1, PR-4B) were highly expressed in C2 of female flowers (Figure 5A and Supplementary Table S6).
Correlation analysis revealed that expression levels of six genes among the selected 14 genes were highly correlated with SA, including ICS2, CM2, NPR5, TGA3, PR-1, and PR-4B, while only two genes (CM3 and NPR3) were significantly correlated with MESA ( Figure 5B). Besides, the qRT-PCR analysis of six genes revealed consistent expression patterns with those generated by RNA-seq data (Figures 5C-H). Based on the above analyses, ICS2, CM2, NPR5, TGA3, PR-1, and PR-4B should be the important regulators in the pathways of SA synthesis and signaling and play a role in the abortion of SFF in tung tree.

Co-expression Networks Analysis of Expression Genes Involved in the Abortion of SFF
To gain more insights into the regulatory relationships of anther and pollen development genes, SA-and PCD-related genes in female flowers of tung tree, we performed a WGCNA of transcript expression in male and female flowers (FPKM values ≥ 1). Consequently, 21 co-expression modules were identified for each sample (Figure 6A). An effective PCC threshold of ≥ 0.85 and a p-value of ≤ 0.05 were trained to generate the significant modules. Finally, we found that each sample of female flowers had one significant module, namely C1 with MEblue module (PCC = 0.91, p-value = 0.01), C2 with MEpink module (PCC = 0.88, p-value = 0.02), and C3 with MEbrown module (PCC = 0.85, p-value = 0.03) (Figure 6A). The MEblue module of the co-expression network, representing the development of SFF, contained 3469 genes. Of the 3469 genes, a number of genes were identified as anther and pollen development genes (TAZ1, MS1, MS2), PCD-related genes (RING1, CBSX5, and POB1), and SA-related genes (CM3, NPR3, RGA1, RGA4) (Supplementary Table S7 (Figure 6B and Supplementary Table S7). Furthermore, WRKYs were identified as having high connectivity between SA-and PCD-related genes, or SA-related genes and anther and pollen development genes, i.e., WRKY41 (3132 edges), WRKY32 (3178), WRKY39 (2992), suggesting their important biological functions in regulating anther and pollen development of female flowers in tung tree ( Figure 6B and Supplementary Table S7).
The MEbrown module, contained 2453 genes and represented the abortion of SFF in tung tree (Supplementary Table S7). Interestingly, only PCD-and SA-related genes were detected in the network, such as AGP16, RAN1, BAG7, FLAVIN-DEPENDENT MONOXYGENASE1 (FMO1), ASGP1, HVA22C, CM1, TGA7, and PAL, etc. (Supplementary Table S7). Furthermore, there were high connectivity among (1) the PCD-related genes, including AGP16 (1872 edges), RAN1 (1509), BAG7 (1925), and FMO1 (1010), (2) the SA-related genes, including CM1 (2114), and PAL (1660) (Figure 6D and Supplementary Table S7). This result suggests that these genes play a key roles in regulating the abortion of SFF in tung tree. Taken together, the analyses of gene co-expression networks suggested that SA might regulate the anther and pollen development and the PCD, resulting in the abortion of SFF in tung tree.

DISCUSSION
The production of unisexual flowers in flowering plants has evolved more than 1000 times (Renner and Feil, 1993). Unisexual flowering systems promote outbreeding and are considered as driving forces in plant evolution. Although unisexual flowering systems have been studied in A. deliciosa, Quercus sube, and Diospyros kaki (Rocheta et al., 2014;Akagi et al., 2016Akagi et al., , 2018, the regulatory mechanisms are unclear.

Anther and Pollen Development Genes Regulate the Abortion of SFF
In tung tree, we found that the SFF was arrested in C2 before the formation of tetrads in microspore meiosis (Figure 1), which was similar with Asparagus officinalis (Caporali et al., 1994). More importantly, the abortion of SFF began with tapetum degeneration. Throughout the abortion process of SFF, the degeneration of tapetum may be a key factor causing the MMC degeneration in SFF. The tapetum plays an important role in pollen grain development which serves as a nutritive tissue, providing metabolites, nutrients, and cell wall precursors for the development of pollen grains, therefore any obstruction of the tapetum development will lead to male sterility in plants (Pacini, 1997;Chen et al., 2018). In tung tree, the tapetum degeneration may result in the abortion of SFF.
The use of classical genetic screens has uncovered a large number of genes involved in anthers development and in the production of viable pollen. In our study, 17 anther and pollen development genes were highly expressed in C1, but lowly expressed or silenced in C2 and C3 (Figure 4). It is worth mentioning that SPL/NZZ, and EXS1/EMS1 have important roles in early anther development and differentiation, while AMS, MS1, MYB33, MYB65, and MYB103 are involved in tapetum and pollen wall development (Pearce et al., 2015). Thus, low expression or silence of anther and pollen development genes, especially the genes of tapetum development, may result in the abortion of SFF in C2 and C3 of female flowers.

PCD Triggers the Abortion of SFF
Programmed cell death is a controlled cellular suicide which plays an important role in plant development. The nuclear DNA degradation, vacuolization, or the activation of specific proteases belongs to specific morphological and biochemical features of PCD (Hautegem et al., 2015). In tung tree, SFF were degenerated when stamens developed into early MMC, and MMC were vacuolize in C2 of female flowers (Figure 1). In Thymelaea and Meliaceae, the gynoecium in male flowers is aborted by PCD at meiosis stage (Caporali et al., 2006;Gouvea et al., 2008). Hence, PCD is essential to many aspects of plant morphogenesis and is a normal component of flower development (Rogers, 2006).
Programmed cell death is controlled by both positive regulators and negative regulators. In this study, positive regulators were highly expressed in C2 or C3 of female flowers, like SOBIR1, ASPG16, and BAG7, which are well-known to positively regulate PCD (Gao and Showalter, 1999;Li et al., 2017;Lu et al., 2017;Sala et al., 2017). Negative regulators showed low expression in C2 and C3, such as RBOHs, which are the source of reactive oxygen species in the oxidative burst and negatively regulate cell death (Torres et al., 2002). These results suggested that PCD-related genes rapidly caused the abortion of SFF in tung tree.  (Boualem et al., 2008;Tao et al., 2018). CK promotes female flowers development in Populus tomentosa, Jatropha curcas, and Castanea henryi Chen et al., 2014;Fan et al., 2017). In tung tree, contents of GA4, IAA, TZR and JA were significantly higher in X3 of male flowers than C3 of female flowers ( Table 1). Only SA was specifically and highly expressed in C2, the initial stage of the abortion of SFF (Table 1). Therefore, the SA accumulation may suppress development of SFF in tung tree.
NPR5, TAG3, PR-1, and PR-4B in SA pathway showed high expression in C2 of female flowers. Interestingly, ICS2 also exhibited the highest expression level in C2 among all samples. In Arabidopsis, isochorismate pathway is the main source of SA accumulation (Wildermuth et al., 2001). Further, the isochorismate pathway also has been shown to be active in tomato and tobacco (Uppalapati et al., 2007). Together, isochorismate pathway may be also the main source of large SA accumulation in C2 of female flowers in tung tree.

SA Mediates the Abortion of SFF by Anther and Pollen Development Genes and PCD-Related Genes
Salicylic acid mediates PCD in plants under biotic stress (Alvarez, 2000). In our study, co-expression networks analysis indicated that anther and pollen development genes, SA-and PCD-related genes showed high connectivity in the network, and that WRKYs were the important connectors between SA-related genes and anther and pollen development genes, and SA-and PCD-related genes (Figure 6). In Arabidopsis, members of WRKYs have interaction with SA-related genes. For example, WRKY58 interacts with NPR1 to regulate SA accumulation in plant defense (Wang et al., 2006). The lower level of SA limits NPR1-NPR3 interaction, enabling NPR1 to restrict the spread of PCD (Fu et al., 2012). SA-dependent signal potentiation loop has been proposed to be negatively regulated by protein LSD1 (Aviv et al., 2002), which also explains the runaway cell death phenotype of the lsd1 mutant. In addition, FMO1 is not only necessary for systemic accumulation of SA and downstream signaling after pathogen infection but also regulates PCD in Arabidopsis (Xu and Brosche, 2014). In rice, SA is suggested to regulate pollen viability and floret fertility (Zhao et al., 2018). Herein, we proposed a possible signaling network where anther and pollen development genes and PCD-related genes work synergistically with SA-related genes to induce the abortion of SFF in tung tree (Figure 6).
Additionally, we constructed a hypothetical model for the development of normal female flowers in tung tree. At the early stage, SFF develops normally in female flowers with low level of SA (Figure 7). With SA accumulation, PCD-related genes are activated and highly expressed, and subsequently anther and pollen development genes, especially tapetum genes are suppressed (Figure 7). Consequently, the abortion of SFF occurs and the normal female flower develops.

CONCLUSION
We confirmed that the abortion of SFF in tung tree belongs to the PCD type and demonstrated that tapetum degeneration at MMC stage is the major reason causing the abortion of SFF. Furthermore, we constructed a model for the abortion of SFF in tung tree based on integrated analyses of morphological and cytological observations, endogenous phytohormone assay and RNA-seq. During the middle stage in female flowers, SA accumulation triggers PCD activation and arrests anther and pollen development, which ultimately results in the abortion of SFF in tung tree. This study provides valuable information for better understanding the development of female flowers in tung tree.

AUTHOR CONTRIBUTIONS
LZ, XT, and ML conceived and designed the experiments. ML, WL, XF, and GZ performed the experiments. ML, XF, HL, MS, GZ, YF, and LZ analyzed the data. ML, WL, GZ, XT, HL, YF, and LZ contributed reagents, materials, and analysis tools. ML and LZ wrote the manuscript. All authors discussed results and commented on the manuscript.