Comparative Transcriptome Analysis Identiﬁed Di ﬀ erentially Expressed Genes between Male and Female Flowers of Zanthoxylum armatum var.

: As a traditional spicy condiment, Zanthoxylum armatum var. novemfolius is of high economical and medicinal value. Despite the long history of human cultivation, the molecular mechanisms underlying ﬂower development are still poorly understood in Z. armatum . In this study, we performed de novo transcriptome assembly and comparative analysis of female and male ﬂowers in Z. armatum . A total of 94,771 unigenes were obtained, and 50,605 unigenes were successfully annotated against the public database. Transcriptome data showed that 20,431 annotated unigenes were di ﬀ erentially expressed genes. KEGG enrichment analysis revealed that the most representative pathway was plant hormone signal transduction. Among them, 41, 16, 41, 27, 95, and 40 unigenes were involved in the biosynthesis and signaling of abscisic acid, ethylene, cytokinin, gibberellin, auxin, and jasmonic acid, respectively. Transcription factors also played crucial roles in ﬂower development, such as AGL11 , PMADS2 , and NAC . These results provided an important basis for characterizing the potential mechanism of ﬂower development and enriching the knowledge of reproduction genetics in Z. armatum .


Introduction
Zanthoxylum armatum var. novemfolius, a variant of Zanthoxylum armatum, belongs to the Rutaceae family, which has high economic value and strong resistance to barren soil and drought [1]. This species was first discovered in the Jiangjin district of Chongqing and is currently widely cultivated in Chongqing, Sichuan, and Guizhou provinces, China. Compared with other species and varieties of Zanthoxylum L., Z. armatum has large mature particles, with an obvious early flowering and fruiting period, high fruit yield, rich aroma, pure taste, and high oil content. Z. armatum is apomixis, dioecious, and the cultivated varieties are mostly pure female plants [2]. However, in recent years, we have found that the female Z. armatum plants (formerly pure female plants) exhibit male, female, and hermaphroditic flowers, leading to a sharp decline in the yield. Therefore, inhibiting the occurrence of male flowers becomes an urgent problem to be solved for the Z. armatum industry. Elucidating the molecular mechanisms underlying the occurrence of male flowers in Zanthoxylum L. is an important basis for achieving this goal.
Flower development is a complex and accurately coordinated biological and morphological process consisting of spatial regulation of a considerable number of organ-specific genes during the life cycle of higher plants [3]. The establishment of male and female flower organs is related to flower organ identity genes, which encode proteins including hormone regulators, transcription factors (TFs), microRNAs, and chromatin-modifying proteins [4,5]. The well-known ABCDE model explains the differentiation of flower organs, and most genes related to flower development are members of the MADS-box family [6,7]. AGAMOUS (AG) controls carpel, AGAMOUS-LIKE (AGL) gene AGL11 plays a role in ovule formation, and APETALA1 (AP1) participates in the transition from vegetative to reproductive growth, specifying sepals and petals. APETALA3 (AP3) and PISTILATA (PI) promote stamens [8]. These genes all play a role in flower development and may lead to the formation of male flowers. In addition to TFs, flower development is also regulated by plant hormones. According to a large number of studies, genes related to plant hormone synthesis and signal transduction pathways have effects on flower development. For example, the ACS gene can inhibit stamen development in female flowers of cucumber [9]. Loss of the tasselseed1 (ts1) gene, which encodes lipoxygenase affecting jasmonic acid synthesis during inflorescence development of maize, leads to a decrease in endogenous jasmonic acid concentration and stamen deficiency [10]. In kiwifruit, Shy Girl (SyGI), a type-C cytokinin response regulator, was identified to be a dominant suppressor of carpel development [11].
The emergence of male flowers in Z. armatum is currently considered to be due to poor pruning, and abuse of fertilizers and plant growth regulators during flower bud differentiation, but little is known about its molecular mechanisms. In this study we, therefore, analyzed male and female flower morphologies and compared the transcriptomes of male and female flowers in Z. armatum using high-throughput sequencing to explore the candidate genes regulating male flower occurrence. Our findings laid the foundation for revealing the molecular mechanism of male flower formation in Z. armatum.

Materials
Zanthoxylum armatum var. novemfolius is planted in germplasm resources of Chinese prickly ash, Rongchang District, Chongqing, China. Female flowers (FFS) and male flowers (MFS) of Z. armatum were collected at least 100 flowers (flower organs only) from different parts of the same tree during the early flowering stage of Z. armatum in March, 2019. The samples were immediately frozen in liquid nitrogen after collection, then transported back to the laboratory and stored at −80 • C until further use. The female and male branches were placed in a water bottle for preservation and used for morphological and cytological observation.

Morphological Observation
The female and male flowers of Z. armatum were photographed by a Cannon EOS 5D Mark digital camera (Tokyo, Japan) and a Leica S9i stereomicroscope (Leica Microsystems, Leica, Wetzlar, Germany) to observe the external morphological features. The samples were prepared for frozen section observation as follows [12]. The female and male flowers of Z. armatum were collected, fixed in 4% (v/v) paraformaldehyde for 48 h at 4 • C, and then transferred to 20% and 30% (w/v) sucrose solution successively for dehydration. The samples were kept in each solution for 24 h to allow them to sink to the bottom. The fixed samples were placed on a sample holder and embedded in optimal cutting temperature (OCT) compound. Tissue sections (10 µm in thickness) were cut using a cryostat, stained with 0.1% (w/v) safranin for 1 h, and dehydrated in a graded ethanol series (50%, 70%, and 80%, 1 min each step). Finally, the sections were dyed with 1% fast green FCF stain (w/v in 95% ethanol) for 1 min, sealed with the neutral gum, and observed under a standard bright field light microscope (Eclipse Ci-L, Nikon Instruments Inc., Tokyo, Japan) with a CCD camera.

RNA Extraction, Library Construction, and RNA-Seq
Total RNA of each sample was isolated using a modified hot borate method and the QIAGEN RNeasy Plant Mini kit (QIAGEN, Hilden, Germany) and used for transcriptome sequencing and qRT-PCR analysis. The quality and quantity of total RNA were assessed using 1% agarose gel electrophoreses and a Nano-Drop ND1000 spectrophotometer (NanoDrop Technologies, Wilmington, DE, USA), and RNA integrity was evaluated using an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). Poly (A) mRNA was isolated and enriched from total RNA using oligo (dT)-attached magnetic beads and then chemically broken into short fragments in fragmentation buffer. DNA fragments with suitable size were selected for PCR amplification to generate the cDNA library. Finally, the libraries were sequenced using the Illumina HiSeq™ 2500 platform at Shanghai Majorbio Bio-pharm Technology Co., Ltd (Shanghai, China).

De Novo Assembly and Functional Annotation
Filtering was performed through the standard bioinformatics procedure established by sequencing company to obtain high-quality clean readings for de novo assembly. Trinity software (https: //github.com/trinityrnaseq/trinityrnaseq/wiki) was employed to pool and assemble all the clean reads [13]. To obtain a functional annotation, the assembled contigs and unigenes of Z. armatum were aligned against six public databases (Nr, KEGG, GO, COG, Swiss-Prot, and Pfam) using BLASTx (e-value < 0.00001).

Analysis of Differentially Expressed Genes (DEGs)
Transcripts per million (TPM) was used to quantify the expression levels of unigenes [14]. The DEGseq R package was employed to identify DEGs from transcriptome data [15]. FDR (false discovery rate) was set as < 0.01, and absolute log2 fold change (FC) ≥1 was set as the threshold for screening the DEGs from the two samples [16].

Identification of Transcription Factors (TFs)
The Transcription Factor Prediction module (http://planttfdb.cbi.pku.edu.cn/) was used to search for the predicted peptide sequences of all unigenes against the TFs in PlantTFDB and to identify the TFs expressed during male flower formation of Z. armatum. This search includes specific information about transcription factor and statistics of transcription factor families.

qRT-PCR Verification
cDNAs were generated using the RevertAid™ First Strand cDNA Synthesis Kit (Thermo Scientific, Waltham, MA, USA). qRT-PCR was performed using a Fast SYBR Mixture (CWBIO, Beijing, China) on a real-time thermal cycler (qTower2.2, Analytic Jena, Leipzig, Germany). ZaUBC and ZaGAPDH were selected as the reference gene [17]. Gene-specific primers were designed using the online primer-design software (http://www.genscript.com.cn/technology-support/on-line-tools), and the primers used in qRT-PCR are shown in Table S1. The 2-∆∆Ct method was used to calculate the relative expression levels of the candidate genes [18]. Three technical replicates and three biological replicates were examined for each sample.

Morphological Analysis of FFS and MFS in Z. armatum
The morphologies of FFS and MFS in Z. armatum were compared by naked eye observation ( Figure 1A-D) and light microscopy observation ( Figure 1E-J). FFS of Z. armatum are small with approximately 3-4 mm in diameter. There were only sepals and pistils in FFS. The sepal is light green in color. The pistil is pale yellowish green in color, which has a superior ovary and two or three carpels ( Figure 1A,C,E-G). MFS is similar in size or slightly larger than FFS. In MFS, there were only sepals and stamens, and the pistils are degenerate. The sepal is green. The stamens are on the inner side of the sepal, ranging in number from 4 to 10 pieces (mainly six pieces), and the anthers are oval and appear bright yellow when the pollen is mature. The anther is butterfly-shaped and divided into two parts by the drug compartment, each of them has two pollen sacs ( Figure 1B,D,H-J). Each ovary of the FFS has a ventricle with two ovules. This structure is loss in MFS (degenerate pistil). The cross-cut structure of receptacle in MFS is similar to that in FFS ( Figure 1E,H). The above results indicated that the pistil breaks down higher up on the receptacle and the stamen occurs under complex control in Z. armatum.
Agronomy 2020, 10, x FOR PEER REVIEW 4 of 21 sepals and stamens, and the pistils are degenerate. The sepal is green. The stamens are on the inner side of the sepal, ranging in number from 4 to 10 pieces (mainly six pieces), and the anthers are oval and appear bright yellow when the pollen is mature. The anther is butterfly-shaped and divided into two parts by the drug compartment, each of them has two pollen sacs ( Figure 1B,D,H-J). Each ovary of the FFS has a ventricle with two ovules. This structure is loss in MFS (degenerate pistil). The crosscut structure of receptacle in MFS is similar to that in FFS ( Figure 1E,H). The above results indicated that the pistil breaks down

Sequencing Analysis and De Novo Assembly
There is unavailable whole genomic information for Z. armatum currently. For exploring the molecular mechanism of flower development, two cDNA libraries for transcriptome analysis were constructed from FFS and MFS of Z. armatum. The two libraries produced 65,533,834 (FFS) and 77,155,592 (MFS) clean reads, respectively. The percentage of Q30 reached more than 93% (Table 1).

Sequencing Analysis and De Novo Assembly
There is unavailable whole genomic information for Z. armatum currently. For exploring the molecular mechanism of flower development, two cDNA libraries for transcriptome analysis were constructed from FFS and MFS of Z. armatum. The two libraries produced 65,533,834 (FFS) and 77,155,592 (MFS) clean reads, respectively. The percentage of Q30 reached more than 93% (Table 1). 94,771 unigenes were found with an average length of 826.88 nt, of which 21,898 (23.11%) were longer than 1000 bp (Table 1, Figure S1). The transcriptome data were stored in the NCBI Sequence Read Archive (Accession Number SRR10376326, SRR10376325). These results indicated that the transcriptome data were reliable for further analysis.
The unigenes were subjected to the COG database to further homologize the gene products. All unigenes were functionally classified, and 8,348 of them were annotated to 24 COG functional categories. As shown in Figure 2D, "translation, ribosomal structure, and biogenesis" had the most gene distribution (508, 6.09%), followed by "posttranslational modification, protein turnover, chaperones" (439, 5.26%) and "carbohydrate transport and biosynthesis" (332, 3.98%). Only a small number of unigenes were divided into "cell motility" and "nuclear structure." Surprisingly, a large part of the unigenes were categorized without specific allocation, such as "general function prediction only" (462, 5.53%) and "function unknown" (179, 2.14%).

DEG Analysis
As shown in Figure 3A, 41,047 unigenes were found commonly expressed in the FFS and MFS of Z. armatum. According to the relative expression levels between FFS and MFS, the DEGs were decomposed into up-regulated and down-regulated genes. A total of 20,431 DEGs were identified between FFS and MFS, of which 11,457 and 8974 were up-regulated and down-regulated, respectively, in MFS compared with those in FFS ( Figure 3B, Figure S2). A total of 3435 DEGs were successfully assigned to 124 KEGG pathways ( Figure 3C). In particular, "plant hormone signal transduction" (map04075), "phenylpropanoid biosynthesis" (map00940), "ribosome" (map03010), "starch and sucrose biosynthesis" (map00500), "plant-pathogen interaction" (map04626), "Oxidative phosphorylation" (map00190), and "MAPK signaling pathway-plant" (map04016) were significantly enriched between FFS and MFS. From these analyses, we obtained the beneficial resources to characterize male flower formation in Z. armatum.

DEG Analysis
As shown in Figure 3A, 41,047 unigenes were found commonly expressed in the FFS and MFS of Z. armatum. According to the relative expression levels between FFS and MFS, the DEGs were decomposed into up-regulated and down-regulated genes. A total of 20,431 DEGs were identified between FFS and MFS, of which 11,457 and 8974 were up-regulated and down-regulated, respectively, in MFS compared with those in FFS ( Figure 3B, Figure S2). A total of 3435 DEGs were successfully assigned to 124 KEGG pathways ( Figure 3C). In particular, "plant hormone signal transduction" (map04075), "phenylpropanoid biosynthesis" (map00940), "ribosome" (map03010), "starch and sucrose biosynthesis" (map00500), "plant-pathogen interaction" (map04626), "Oxidative phosphorylation" (map00190), and "MAPK signaling pathway-plant" (map04016) were significantly enriched between FFS and MFS. From these analys

Simple Sequence Repeats (SSRs) and Single Nucleotide Polymorphisms (SNPs) Detection
High-throughput transcriptome sequencing is helpful for the discovery of economical and effective polymorphic genetic markers (SSRs and SNPs) [19]. SSRs are a series of simple, repetitive short sequences that are important for genome organization and phenotypic diversity and are one of the most popular molecular markers [20,21]. Among the unigenes, 23,579 SSR markers were identified by using MISA. According to the motif repetition, these microsatellites were further divided into mononucleotide repeats (14,450, 61.28%), followed by di-(4346, 18.43%), tri- (4,194,17.79%), tetra-(341, 1.45%), hexanucleotides motifs (157, 0.67%), and penta-(91, 0.39%) ( Figure 4A, Table S2). SNPs refer to the variation of a single nucleotide in a base sequence, including substitution and transversion [22]. 205,139 SNPs were identified by using SAM tools and BCF tools within all unigenes, of which 57,422 were Homo SNPs, 147,717 were Hete SNPs, and most were below 30 bp. A/G occurred most frequently than other variations ( Figure 4B,C). This experiment greatly enriched the development of SSRs and SNPs in the transcriptome of Z. armatum.
Agronomy 2020, 10, x FOR PEER REVIEW 10 of 21 Rich factor, the greater the degree of enrichment. The size of the point indicated how many unigenes were in the path, and the color of the point corresponded to a different p-value range.

Simple Sequence Repeats (SSRs) and Single Nucleotide Polymorphisms (SNPs) Detection
High-throughput transcriptome sequencing is helpful for the discovery of economical and effective polymorphic genetic markers (SSRs and SNPs) [19]. SSRs are a series of simple, repetitive short sequences that are important for genome organization and phenotypic diversity and are one of the most popular molecular markers [20,21]. Among the unigenes, 23,579 SSR markers were identified by using MISA. According to the motif repetition, these microsatellites were further divided into mononucleotide repeats (14,450, 61.28%), followed by di-(4346, 18.43%), tri- (4,194,17.79%), tetra-(341, 1.45%), hexanucleotides motifs (157, 0.67%), and penta-(91, 0.39%) ( Figure  4A, Table S2). SNPs refer to the variation of a single nucleotide in a base sequence, including substitution and transversion [22]. 205,139 SNPs were identified by using SAM tools and BCF tools within all unigenes, of which 57,422 were Homo SNPs, 147,717 were Hete SNPs, and most were below 30 bp. A/G occurred most frequently than other variations ( Figure 4B,C). This exper

Phytohormone-Related DEGs
Phytohormone is an essential factor for flower sexual determination in various plant species [23]. Some hormone-related genes were detected to investigate the roles of hormones during flower development ( Figure 5, Table S3). 9-cis-epoxycarotenoid dioxygenases (NCED), xanthoxin dehydrogenase (ABA2), and abscisic-aldehyde oxidase (AAO3), key genes involved in abscisic acid biosynthesis pathway, were down-regulated in FFS ( Figure 5A,B). Based on the transcriptome data, 16 key genes were identified in the ethylene pathway. Among them, S-adenosylmethionine synthetase (metK), ACC synthases (ACS), and ethylene-responsive transcription factor (ERF) were up-regulated in MFS compared with those in FFS, while ACC oxidases (ACO) was down-regulated. In particular, two transcripts encoding ethylene receptor ETR homolog (ZaETR1 and ZaETR2) were up-regulated in MFS compared with those in FFS, but one transcript (ETR3) showed remarkable down-regulation ( Figure 5C,D). Seventeen DEGs were related to the gibberellin biosynthesis pathway, and ten DEGs were associated with the gibberellin signal transduction pathway. Among them, gibberellin 20 oxidase (GA20ox) had a higher expression in MFS than in FFS. However, ent-kaurene synthase (KS) and gibberellin receptor (GID1) had higher expressions in FFS than in MFS ( Figure 5E,F). Indole-3-acetaldehyde oxidase (AAO1), a key enzyme involved in auxin biosynthesis pathway, exhibited remarkable up-regulation, but TIR1, transport inhibitor response of auxin signaling, underwent down-regulation ( Figure 5G,H). From the Z. armatum transcriptome data, cytokinin hydroxylase (CYP735A) and histidine-containing phosphotransfer peotein (AHP) in the cytokinin pathway were significantly highly expressed in MFS. However, Arabidopsis histidine kinase (AHK) was highly expressed in FFS ( Figure 5I,J). 12-oxophytodienoic acid reductase (OPR), enoyl-CoA hydratase (MFP2), and acetyl-CoA acyltransferase 1 (ACAA1), involved in the jasmonic acid biosynthesis pathway, were up-regulated; Jasmonate-amino acid synthetase (JAR1) and transcription factor MYC2, key genes involved in jasmonic acid signal transduction pathway, were down-regulated in MFS ( Figure 5K,L). These results showed that genes related to plant hormones may be key genes in flower development of Z. armatum.

Expression Profiling of Candidate Transcription Factors
TFs are involved in flower development and sex differentiation [24]. In the Z. armatum  Table S4). These results were similar to the Ribes diacanthum pall transcriptome data, where bZIP, MYB, WRKY, bHLH, zinc finger, and PHD finger proteins were the largest TF groups [25], and to another transcriptome data where the largest TF groups are MYB, bHLH, NAC, AP2/ERF, and MADS-box [26]. These results implied that MYB, AP2/ERF, bHLH, and WRKY were TF superfamilies in plants. The AP2/ERF and MADS-box families substantially involved in regulation of plant flower development.
The MADS-box gene family is a key transcription factor with highly conserved MADS-box DNA binding domain, which is of great importance for the regulation of plant growth and development [27,28]. Here, six putative MADS-box genes related to flower development were identified by qRT-PCR, including AGAMOUS-like, which regulates stamen and carpel development; putative gene encoding floral-binding protein 24, which defines flower development; putative gene encodes floral homeotic protein PMADS2, which relates to flower organ formation; putative gene encoding floral homeotic protein APETALA1, which is mainly responsible for the development of flower buds and petals. Some of these candidate genes such as ZaPMADS2, ZaAP1, and ZaAGL20 were lowly expressed in FFS. However, others such as ZaAGL11 were highly expressed in FFS (Figure 7, Table S4).  In addition to the MADS-box gene, other transcription factor genes such as MYB, WRKY, and NAC were also differentially expressed between FFS and MFS. DN39886_c2_g2 encoding TRB4, DN32791_c2_g2 encoding transcription factor MYB44, and DN30645_c0_g1 encoding NAC56 were up-regulated in MFS compared with those in FFS. However, MYB12 and WRKY70 had higher transcript levels in FFS than in MFS (Figure 7, Table S4).

Validation of Candidate Genes by qRT-PCR
To confirm the reliability of transcriptome data, we randomly selected 20 genes for qRT-PCR validation ( Figure 7A, Table S1). The expression patterns of the tested genes were mostly consistent with the transcriptome data. Meanwhile, correlation analysis between qRT-PCR and RNA-Seq results showed that the two were highly correlated (R 2 = 0.7724, Figure 7B). These results indicated that the transcriptome data were reliable and accurate.

Discussion
Z. armatum is a traditional condiment mainly consumed in China. However, genomic information on this species is scarce. Comparative transcriptome analysis of FFS and MFS found that 94,771 unigenes with a mean size of 826.88 nt were generated and assembled, and 50,605 unigenes were successfully annotated against the public database. The species distribution of the annotated unigenes showed that Z. armatum and Citrus sinensis had the highest similarity (43.68%) and belonged to the same family (Rutaceae). These results suggested that transcriptome assembly and annotation were reliable and provided a basis for screening genes associated with flower development of Z. armatum.
According to the annotations of DEGs between FFS and MFS in Z. armatum, the candidate regulators participating in male flower formation were screened. For example, DN35970_c5_g1, has similarity to putative Male Sterility 1 (MS1) in Arabidopsis, which controls anther and pollen development and plays a key role in regulating male gametogenesis [29]. DN31825_c0_g1 and DN31825_c0_g2 are similar to a putative pollen-specific leucine-rich repeat extensin-like protein 4 (PEX4), which is specifically expressed in flowers, stamen, pollen, and pollinated carpels of Arabidopsis [30]. DN31265_c2_g3 and DN31265_c2_g1 have similarity to FIL1, which are observed to be stamen-specific protein [31]. DN33279_c2_g2, similar to Aborted Microspores (AMS) which is essential for tapetum and pollen development in Arabidopsis [32]. These unigenes are specifically detected in male flowers of Z. armatum and might also be involved in regulating flower development Z. armatum.
Duan et al. [33] found that the genes involved in Pinus bungeana male and female conelet development were related to tryptophan biosynthesis, zeatin biosynthesis, and cysteine and methionine biosynthesis, the phenylpropanoid biosynthesis pathways and some important plant hormone pathways. Carotenoids are the precursor of abscisic acid and coordinate with gibberellin to regulate plant growth and development. The present study showed that some genes participating in the ABA biosynthetic and signaling pathway were significantly differentially expressed. Among them, ZaNCEDs, ZaABA2, and ZaAAO3 were remarkably up-regulated in MFS, suggesting that ABA might play a potential role in male flower development. As an inducing signal, gibberellin regulates the growth and development of higher plants [34]. Exogenous GA is often applied to regulate plant sexual differentiation. Researchers have different opinions on the effect of GA on plant gender expression. Chen et al. [35] believed that GA seriously affects the development of male flowers in cucumber. However, Zhao et al. [36] found that GA could promote the development of female flowers in Momordica charantia. Our data showed that the expression levels of genes associated with GA biosynthetic and signaling were notably altered in MFS, compared with those in FFS, such as GA20oxs, GA2oxs, GID1s and DELLAs, suggesting the involvements of GA in modulating flower development. Two ZaGA20ox genes (DN36887_c1_g5 and DN37407_c0_g1), key enzyme related to GA biosynthesis, were obviously up-regulated in MFS, indicating GA might promote male flower occurrence in Z. armatum. Ethylene plays a positive feminizing role in cucumber [37], melon [38], and bitter gourd [39]. We found that ZametK, ZaACS, and ZaERF were significantly down-regulated in FFS compared with those in MFS, but ZaACO was up-regulated. This result showed that the effect of ethylene on flower development in Z. armatum might be different from that in cucumber, melon, or cactus, and detailed mechanism needs further investigation. Auxin plays a vital role in organ construction and cell differentiation. When its synthesis, transport, and signal transduction are disturbed, defects in the pistil are generated. Auxin has certain female effects on Trichosanthes kirilowii [40] and melon [41]. GH3 genes, encoding IAA-amido synthetases, is important to maintain auxin homeostasis [42]. In the present study, several genes, such as ZaGH3-6 and ZaGH3-7, were up-regulated in FFS. The auxin responsive AUX/IAA protein family genes are repressors of early auxin response genes at low auxin concentrations [43]. In tomato, Indole Acetic Acid 9 (IAA9) is a negative auxin response regulator belonging to the AUX/IAA transcription factor gene family whose downregulation triggers fruit set before pollination, thus giving rise to parthenocarpy, and a tissue-specific gradient of IAA9 expression is established during flower development, the release of which upon pollination triggers the initiation of fruit development [44]. In mature embryo sacs, overexpression of the auxin-producing enzyme YUC1 and down-regulation of auxin-dependent transcription factors (ARFs) alter embryo sac cell identity [45]. Similarly, DN30775_c0_g5 and DN29744_c0_g4 annotated as IAA9 and DN38175_c0_g5 annotated as ARF were up-regulated in FFS in Z. armatum. Therefore, a large number of genes involved in auxin pathway may potentially lead to female flower development in Z. armatum. Cytokinin is a positive regulator exerting a feminizing effect in cucumber [35] and Vitis amurensis [46]. In this study, we observed that the large proportion of differentially expressed ZaCKXs, genes involved in oxidative cytokinin degradation, exhibited outstanding downregulation in FFS compared with those in MFS, implying that the suppression of cytokinin degradation might contribute to promote female development in Z. armatum. The tasselseed1 (ts1) encodes lipoxygenase that synthesizes jasmonic acid and its absence leads to a decrease in endogenous jasmonic acid concentration and stamen loss in maize [10]. opr7, opr8 display strong developmental defects, including formation of a feminized tassel, initiation of female reproductive buds at each node, and extreme elongation of ear shanks; these defects are rescued by exogenous jasmonic acid [47]. In Z. armatum, DN37710_c0_g1 (LOX2S) and DN24178_c0_g1 (OPR) were up-regulated in MFS, indicating that jasmonic acid may promote male flower formation.
TFs are important regulators in plant development. One class of TFs is MADS-box, now known as the homeotic genes of ABCDE model for flower organ identification, which explains how the combination of different genes determines the identification of four flower organs. MADS-box genes regulate the early stage of flower development and help determine the characteristics of floral organs [48,49]. Some MADS-box genes were identified in the transcriptome data. To explore the evolutionary relationships of the Z. armatum MADS-box genes, a phylogenetic tree was constructed using conserved cDNAs and published genes from other species ( Figure 6C, Table S4). Our results showed that the putative MADS-box gene MADS3 (DN40845_c0_g2) especially expressed in male flower development in Z. armatum, similar to PhPMADS3, which is expressed exclusively in stamens and carpels of wild-type petunia plants [50]. AGL11, FBP24 and PMADS2 are highly and specifically expressed during ovule formation in petunia and Arabidopsis [51]. Similarly, ZaAGL11 (DN29192_c1_g1 and DN37221_c0_g7) and ZaFBP24 (DN27652_c1_g2) showed up-regulated expression in FFS, but ZaPMADS2 (DN39331_c0_g2) showed down-regulated expression in FFS. The reason for this difference may be the variation of species. Further research is needed to solve this problem. Other TFs, such as MYB, WRKY, and NAC, were also involved in the regulation of flower development [52][53][54]. NAC-domain containing TFs are essential for normal plant morphogenesis which representatives are the Cup-Shaped Cotyledon proteins encoded by CUC2 (DN31092_c0_g9) genes more expressed in male flowers. In Silene latifolia, CUC2 is specifically expressed in male organs before the pistil stop development [55]. The role of MYB transcription factor superfamily in plant reproductive development is unequivocal. MYB proteins include a conserved domain, the MYB DNA-binding domain regulating a variety of plant-specific processes and they may play roles in different plant species. For example, R2R3-MYB and MYB88 regulate female flower development in Arabidopsis thaliana [56]. Importance in seed coat formation, trichome morphogenesis and mucilage secretion of MYB5 is also described by several studies [57]. In our investigations, R2R3-MYB (DN38469_c0_g4 and DN25343_c0_g2) and MYB5 (DN34758_c2_g1 and DN39253_c0_g1) showed unique expression in female flowers, thereby indicating that TFs play regulatory roles in flower development. Although many genes related to flower development are screened out, the molecular mechanism and expression characteristics of these candidate genes are still uncertain and need further study. Moreover, miRNAs and DNA methylation also play an essential role in the flower development and require in-depth investigation [58].

Conclusions
In this study, transcriptome data were analyzed, and several candidate genes related to male flower formation were screened in Z. armatum. Our results showed that the flowers of Z. armatum might begin to develop in the form of unisexual flower (female flower) until the ovary disappeared, and the male flower gradually appeared, thereby transforming from a female flower to a male flower. Differential expression analysis also provided evidence that the expression of genes such as AGL11, MYB12, and WRKY70 suddenly decreased, whereas that of PMADS2 sharply increased in male flowers, illustrating that these candidate genes might be involved in the male flower formation. The functional annotation of high-throughput transcriptome data obtained in this study may be used as genetic resources for formulating related plant improvement strategies such as for Z. armatum. Further study of selected candidate genes may enable researchers and breeders to select or control specific sexual expression patterns in Z. armatum.