Transcriptome Analysis Reveals Important Transcription Factor Families and Reproductive Biological Processes of Flower Development in Celery ( Apium graveolens L.)

: There are few reports on the reproductive biology of celery, which produces small ﬂowers in a long ﬂowering period. Anther development was analyzed by para ﬃ n sectioning and related genes were examined by transcriptome sequencing and qPCR. The development process was divided into nine stages based on the signiﬁcant changes in the cell and tissue morphologies. These stages included: archesporial stage, sporogenous cell stage, microspore mother cell stage, dyad and tetrad stage, mononuclear microspore stage, late uninucleate microspore stage, binuclear cell stage, mature pollen stage, and dehiscence stage. A total of 1074 di ﬀ erentially expressed genes were identiﬁed by transcriptome sequencing in the early ﬂower bud, middle ﬂower bud, and early ﬂowering period. Functional annotation indicated that these genes were involved in physiological and biochemical processes such as ribosomes metabolism, sugar metabolism, and amino acid metabolism. Transcription factors such as C2H2, AP2 / ERF, bZIP, WRKY, and MYB played key regulatory roles in anther development and had di ﬀ erent regulatory capabilities at various stages. The expression patterns based on qPCR and transcriptome data of the selected transcription factor genes showed consistency, suggesting that these genes played an important role in di ﬀ erent ﬂower development stages. These results provide a theoretical basis for molecular breeding of new celery varieties with pollen abortion. Furthermore, they have enriched research on the reproductive biology of celery and the Apiaceae family.


Introduction
In flowering plants, the anther is where male gametophytes (pollen) are produced. Anthers of angiosperms are composed of four pollen sacs connected by an anther chamber. The anther wall

Cytological Study on Anthers of Celery
Fresh anthers were used for cytological study. The anthers were divided into nine stages based on their size. They were added to 8 mL of 50% FAA fixative solution (V formalin: V glacial acetic acid: V50% alcohol = 1: 1: 18) and fixed at 4 • C for 24 h, dehydrated with ethanol and then embedded in paraffin. The thickness of the slice was 5 µm; safranin O-fast green reagent was used to stain anther tissues. Morphological structure at each developmental stage was observed by optical microscope (Nikon Eclipse Ci, Tokyo, Japan) and photographed [23].

RNA Extraction and Transcriptome Sequencing
Total RNA was extracted from the samples of each stage using the total RNA kit (Tiangen, Beijing, China) according to the manufacturer's instructions. The purity and integrity were tested using the NanoPhotometer spectrophotometer and Agilent 2100 bioanalyzer, respectively. Oligo (dT) was used to enrich the mRNA with a polyA tail. The mRNA was then fragmented and the first strand of cDNA was synthesized with random oligonucleotides and reverse transcriptase. RNaseH and dNTPs were added to synthesize the second strand of cDNA in the DNA polymerase I system. The purified double-stranded cDNA was amplified by PCR and the products were further purified to obtain libraries for Illumina HiSeq 4000 sequencing (Novogene Technologies Co., Ltd., Beijing, China). The transcriptome data were submitted to NCBI under the accession number SRP250986.

De Novo Assembly
Raw data were filtered to obtain clean reads by removing the joint, N-containing, low-quality reads (the base number Q ≤ 20 accounted for more than 50% of the entire read length). Assembly of clean reads was performed by Trinity while the contig sequences of each sample were obtained by Overlap [24]. Contig sequences obtained by paired-end were further assembled to obtain the transcripts of each sample. Corset software was used to aggregate transcripts into several clusters based on shared reads among the transcripts [25]. Transcripts with differential expression among the samples were separated from the original cluster to establish a new cluster. The longest cluster sequence was selected as the unigenes for subsequent analysis.

Functional Annotation
Annotation results were obtained by blasting unigenes sequence with Nr (NCBI protein database), Nt (NCBI nucleic acid sequence database), Pfam (protein domain annotation classification system), Swiss-Prot database, KOG (lineal homology database), GO (genomic theory database), and KEGG (Kyoto gene and genome encyclopedia database). Blast2GO was used to annotate and classify GO terms.

Analysis of Transcript Abundance and Selection of Differentially Expressed Genes
Transcript abundance of all unigenes in the three stages was estimated by calculating the read density as fragments per kilobase per million (FPKM). Differentially expressed genes (DEGs) were screened by DESeq2 with the selection criteria set as |log2 (FoldChange)| > 1 and padj < 0.05 [26]. When the expression of a gene in the two groups of samples was more than double, the differential expression was determined as DEG.

Analysis of Transcription Factors for Differentially Expressed Genes
ITAK software and PlantTFDB (http://planttfdb.cbi.pku.edu.cn) were used to identify celery anther transcription factors and screen different gene transcription factor families for classification [27]. Furthermore, a cluster heat map was drawn to represent the expression of these genes.

Verifying Gene Expression Levels by qPCR
Transcription factor families related to flower development were selected for further analysis of their expression levels at different developmental stages based on their results in DEGs. Double-stranded cDNA samples of celery flower at S1, S2, and S3 stage produced for transcriptome sequencing were used as templates for qPCR. Primer Premier 6.0 was used for designing specific primers. The experiments were performed on a BIO-RAD CFX96 quantitative PCR instrument (Bio-rad, Hercules, CA, USA) according to the method of the SYBR Premix Ex Taq kit (TransGene, Beijing, China). TUB was selected as the reference gene; the relative expression level of genes was calculated using the method of 2 −∆∆CT [28]. SPSS23.0 was used for significance analysis and GraphPad Prism 8.0 was used for plotting.

Characteristics of Anther Development in Celery
Anther development in celery was divided into nine stages. The names of each stage and the corresponding cell morphological changes are shown in Table 1 and Figure 1. Between stages 2 and 3, the epidermis differentiates to form four anther chambers. Archesporial cells differentiate into primary parietal cells and primary sporophytic cells. The former differentiates into the endothecium (fibrous layer), the middle layer, and the tapetum, whereas the latter divides to form secondary sporoblasts, which can develop into microspore mother cells. At the period of the microspore mother cell, the four-cell layers of the anther which are the tapetum, middle layer, fibrous layer, and epidermis can be seen. Moreover, callose deposition gradually began between the plasma membrane and the primary wall, and further extended from the corner until it surrounded the entire mother cell. With the development of anthers, microspore mother cells underwent meiosis at the fourth and fifth stages, forming dyad and tetrad. At this point, the tapetum had reached its maximum and most active stage while the callose gradually degraded and microspores dissociated. If the callose does not degrade during this period, it will lead to abnormal meiosis, which will affect the release of microspores, resulting in abnormal development of the anther and subsequent formation of male-sterile plants. Following the late uninucleate microspores period, the cells undergo mitosis twice to form a large vegetative cell and two small sperm cells. The pollen sac contained three cells with mature pollen. The buds are tender yellow and uncracked. The length of the petals is 0.05-0.1 mm 1 Archesporial stage The shape of the anther was oval, and the epidermis has formed.
The buds are tender yellow and uncracked. The length of the petals is 0.07-0.1 mm 2 Sporogenous cell stage The anther epidermis gradually differentiated into 4 anther chambers, the archesporial cells differentiated into primary parietal cells and primary sporophytic cells.
The buds are yellow-green and uncracked. The length of the petals is 0.13-0.16 mm 3 Microspore mother cell stage From the surface layer of the anther, the tapetum, middle layer, fiber layer, and epidermis can be clearly seen from the inside to the outside. As the microspore mother cell develops, thick callose wall on the surface of the cell continued to thicken and the original cellulose wall gradually degraded.
The buds are turquoise and uncracked. The length of the petals is 0.17-0.20 mm 4 Dyad and tetrad stage Microspore mother cells underwent meiosis to form dyads and tetrads. The middle cells were squeezed by the tapetum and outer layers and gradually disappeared. The tapetum cells were the largest and most active.
The buds are pale green and uncracked. The length of the petals is 0.18-0.20 mm 5 Mononuclear microspore stage The callose wall around the tetrad degraded, and the microspores are freed out and their nucleus is centered with their nuclei centered. The middle layer cells were almost absent.
The buds are pale green and uncracked. The length of the petals is 0.20-0.22 mm 6 Late uninucleate microspores stage The vacuoles in the microspore cells gradually increased to form larger vacuoles which pushed the nucleus to the edge of the cell. The tapetum gradually disintegrated.
The buds are half-cracked, the petals are pale white and the length is 0.22-0.25 mm 7 Binuclear cell stage The nucleus on the side divided into two nuclei by mitosis and formed two daughter cells of different sizes. The larger cells were vegetative cells and the smaller germ cells. The tapetum continued to decompose and this was the peak of tapetum decomposition.
The bud is fully open, the petals are white and the length is 0.27-0.29 mm 8 Mature pollen stage The germ cells of the microspore formed a trinuclear pollen containing one vegetative cell and two sperm cells through mitosis.
The space between the upper and lower chambers deteriorates, forming two chambers. The tapetum was very small or disappears completely.
The bud is fully open, the petals are white and the length is 0.27-0.29 mm 9 Dehiscence stage The epidermal cells in the stomata area degenerated and the cells opened causing mature pollen grains to be released.

Transcriptome Sequencing and Data Assembly
To understand the molecular mechanisms behind the progression of flower development in celery, gene expression was analyzed in three flower stages: early flower bud stage (S1), middle flower bud stage (S2), and early flowering stage (S3). A total of 173,094 transcripts were obtained using Trinity and 82,225 unigenes were selected by clustering the transcripts via Corset ( Table 2). The length of unigenes ranged from 301 bp to 11,982 bp, the average length was 1117 bp, and the length of N50 was 1679 bp. Among all unigenes, the range of 300-500 bp accounted for the largest proportion of length, at 32.30%. Unigenes longer than 2,000 bp accounted for the smallest proportion of length, at 15.44%.

Transcriptome Sequencing and Data Assembly
To understand the molecular mechanisms behind the progression of flower development in celery, gene expression was analyzed in three flower stages: early flower bud stage (S1), middle flower bud stage (S2), and early flowering stage (S3). A total of 173,094 transcripts were obtained using Trinity and 82,225 unigenes were selected by clustering the transcripts via Corset ( Table 2). The length of unigenes ranged from 301 bp to 11,982 bp, the average length was 1117 bp, and the length of N50 was 1679 bp. Among all unigenes, the range of 300-500 bp accounted for the largest proportion of length, at 32.30%. Unigenes longer than 2000 bp accounted for the smallest proportion of length, at 15.44%.

Functional Annotation for Unigenes
All unigenes were successfully annotated against seven databases: Nr, Nt, Pfam, KOG, Swiss-Prot, KEGG, and GO. The ratio of successful annotation of each database was obtained which provided a reference for other studies (Table S1) and the functional annotation information of all databases is summarized in Table S2. A total of 39,699 unigenes (48.28% of total genes) could be characterized by GO annotation ( Figure S2): "Biological process" ontology could be assigned to the largest number of genes, followed by "Cellular component" and "Molecular function". In the biological process, the number of genes in cellular processes accounted for the largest proportion, followed by metabolic processes and the single-organism process, and the lowest was cell aggregation. A total of 19,068 (23.19% of total genes) unigenes had annotation information in KOG database I with 26 categories. These included general function prediction, signal transduction mechanisms, RNA processing and modification, defense mechanisms transcription, and other physiological processes ( Figure S3). A total of 17,276 (21.01% of total genes) unigenes were successfully annotated in the KEGG database with the most frequently represented pathways including Cellular Processes, Environmental Information Processing, and Genetic Information Processing ( Figure S4).

Identification and Selection of the Differentially Expressed Genes (DEG)
Using DESeq2 software, 39,218 DEGs were identified in at least one comparison (S1 vs. S2, S1 vs. S3, S2 vs. S3) ( Figure 2A and Table S3). The number of DEGs in S1 vs. S3 was the largest, with a total of 19,849, followed by S2 vs. S3, then S1 vs. S2, with the number of genes being 15,631 and 3738, respectively. There were also 1074 genes that were differentially expressed in the three combinations. All these results suggested that the transcriptome in celery flowers was increasingly different as development progressed. Comparing each combination of the DEGs (including upregulation and downregulation) showed ( Figure 2B) that the number of upregulated DEGs was 1,637 and downregulated DEGs were 2,101 at S1 vs. S2. However, the upregulated numbers of S1 vs. S3 and S2 vs. S3 genes were 17,965 and 15,177, respectively, and downregulated numbers were 1,884 and 454 respectively. This further suggested that differential genes expression was apparent during flower development and, in particular, that an initial reorganization of gene expression was followed by the upregulation of a great number of transcripts in later stages.

Functional Annotation for Unigenes
All unigenes were successfully annotated against seven databases: Nr, Nt, Pfam, KOG, Swiss-Prot, KEGG, and GO. The ratio of successful annotation of each database was obtained which provided a reference for other studies (Table S1) and the functional annotation information of all databases is summarized in Table S2. A total of 39,699 unigenes (48.28% of total genes) could be characterized by GO annotation ( Figure S2): "Biological process" ontology could be assigned to the largest number of genes, followed by "Cellular component" and "Molecular function". In the biological process, the number of genes in cellular processes accounted for the largest proportion, followed by metabolic processes and the single-organism process, and the lowest was cell aggregation. A total of 19,068 (23.19% of total genes) unigenes had annotation information in KOG database I with 26 categories. These included general function prediction, signal transduction mechanisms, RNA processing and modification, defense mechanisms transcription, and other physiological processes ( Figure S3). A total of 17,276 (21.01% of total genes) unigenes were successfully annotated in the KEGG database with the most frequently represented pathways including Cellular Processes, Environmental Information Processing, and Genetic Information Processing ( Figure S4).

Identification and Selection of the Differentially Expressed Genes (DEG)
Using DESeq2 software, 39,218 DEGs were identified in at least one comparison (S1 vs. S2, S1 vs. S3, S2 vs. S3) ( Figure 2A and Table S3). The number of DEGs in S1 vs. S3 was the largest, with a total of 19,849, followed by S2 vs. S3, then S1 vs. S2, with the number of genes being 15,631 and 3738, respectively. There were also 1074 genes that were differentially expressed in the three combinations. All these results suggested that the transcriptome in celery flowers was increasingly different as development progressed. Comparing each combination of the DEGs (including upregulation and downregulation) showed ( Figure 2B) that the number of upregulated DEGs was 1,637 and downregulated DEGs were 2,101 at S1 vs. S2. However, the upregulated numbers of S1 vs. S3 and S2 vs. S3 genes were 17,965 and 15,177, respectively, and downregulated numbers were 1,884 and 454 respectively. This further suggested that differential genes expression was apparent during flower development and, in particular, that an initial reorganization of gene expression was followed by the upregulation of a great number of transcripts in later stages.  The trend seen in Figure 2B was confirmed when the expression levels of the DEGs were clustered ( Figure 3A). Specifically, the expression level of most DEGs was higher in S3 stage than in the S1 and S2 stages. As shown in Figure 3B, the expression level of 1,034 genes steadily increased between the three stages (subcluster 1). These genes were upregulated from S1 to S3 and the increase from S1 to S2 was larger than that from S2 to S3. A total of 4,361 genes were upregulated between S2 and S3 stages (subcluster 2), while the expression levels of 11,339 genes had a slight upward trend (subcluster 3). Moreover, 420 genes showed a slight downward trend from S1 to S2, but were significantly downregulated at S3 (subcluster 4). A total of 229 genes were significantly downregulated from S1 to S2, but not significantly at S3 (subcluster 5). In addition, 5,023 genes were significantly downregulated during S1 to S2 while upregulated during S3 (subcluster 6). These results further indicate the complexity of gene regulation during celery flower development. The trend seen in Figure 2B was confirmed when the expression levels of the DEGs were clustered ( Figure 3A). Specifically, the expression level of most DEGs was higher in S3 stage than in the S1 and S2 stages. As shown in Figure 3B, the expression level of 1,034 genes steadily increased between the three stages (subcluster 1). These genes were upregulated from S1 to S3 and the increase from S1 to S2 was larger than that from S2 to S3. A total of 4,361 genes were upregulated between S2 and S3 stages (subcluster 2), while the expression levels of 11,339 genes had a slight upward trend (subcluster 3). Moreover, 420 genes showed a slight downward trend from S1 to S2, but were significantly downregulated at S3 (subcluster 4). A total of 229 genes were significantly downregulated from S1 to S2, but not significantly at S3 (subcluster 5). In addition, 5,023 genes were significantly downregulated during S1 to S2 while upregulated during S3 (subcluster 6). These results further indicate the complexity of gene regulation during celery flower development.

Functional Annotation of Differentially Expressed Genes
To further study the biological significance of gene variation between flower development stages, GO and KEGG database analyses were performed. Figure 4 shows that in the biological process, the DEGs were enriched in chitin and glucosamine-containing compound metabolic processes. DEGs were significantly enriched in S1 vs. S3 and S2 vs. S3 were associated with the amino sugar metabolic process. DEGs that were significantly enriched in S1 vs. S2 and S1 vs. S3 were linked to the ribonucleoprotein complex and ribosome biosynthesis, while DEGs significantly enriched in S1 vs. S3 participated in the lipid transport process. DEGs that were significantly enriched in S1 vs. S3 and S2 vs. S3 were associated with the cytoplasm and cytoplasmic part. DEGs found to be significantly enriched in S1 vs. S2 and S1 vs. S3 were involved in extracellular regions, ribonucleoprotein complexes, ribosomes, and extracellular regions. DEGs found to be significantly enriched in S2 vs. S3 were involved in threonine-type endopeptidase activity and threonine-type peptidase activity, while those significantly enriched in S1 vs. S2, S1 vs. S3, and S2 vs. S3 were associated with the structural constituent of cuticle, structural molecular activity, and chitin binding. To further study the biological significance of gene variation between flower development stages, GO and KEGG database analyses were performed. Figure 4 shows that in the biological process, the DEGs were enriched in chitin and glucosamine-containing compound metabolic processes. DEGs were significantly enriched in S1 vs. S3 and S2 vs. S3 were associated with the amino sugar metabolic process. DEGs that were significantly enriched in S1 vs. S2 and S1 vs. S3 were linked to the ribonucleoprotein complex and ribosome biosynthesis, while DEGs significantly enriched in S1 vs. S3 participated in the lipid transport process. DEGs that were significantly enriched in S1 vs. S3 and S2 vs. S3 were associated with the cytoplasm and cytoplasmic part. DEGs found to be significantly enriched in S1 vs. S2 and S1 vs. S3 were involved in extracellular regions, ribonucleoprotein complexes, ribosomes, and extracellular regions. DEGs found to be significantly enriched in S2 vs. S3 were involved in threonine-type endopeptidase activity and threonine-type peptidase activity, while those significantly enriched in S1 vs. S2, S1 vs. S3, and S2 vs. S3 were associated with the structural constituent of cuticle, structural molecular activity, and chitin binding. Corrected P-Value < 0.05 represents enrichment, the smaller the P-Value, the closer the color is to red and the more significant the enrichment. Figure 5 reveals that KEGG annotations for the three pairs (S1 vs. S2, S1 vs. S3, and S2 vs. S3) were significantly enriched in proteasome and ribosome pathways. DEGs that were significantly enriched in S1 vs. S2 were involved in oxidative phosphorylation, carbon fixation in photosynthetic organisms, pentose and glucuronate interconversions, plant-pathogen interaction, pentose phosphate pathway, glycolysis/gluconeogenesis, and other metabolic processes. We speculate that these processes regulated nutrients levels during anther development. There was no obvious enrichment in these processes in S1 vs. S3 and S2 vs. S3.In contrast, DEGs that were not enriched in S1 vs. S3 and S2 vs. S3 were associated with processes that regulate nutrients levels and were enriched in the phagosome and spliceosome processes. DEGs found to be only enriched in S2 vs. S3  5 reveals that KEGG annotations for the three pairs (S1 vs. S2, S1 vs. S3, and S2 vs. S3) were significantly enriched in proteasome and ribosome pathways. DEGs that were significantly enriched in S1 vs. S2 were involved in oxidative phosphorylation, carbon fixation in photosynthetic organisms, pentose and glucuronate interconversions, plant-pathogen interaction, pentose phosphate pathway, glycolysis/gluconeogenesis, and other metabolic processes. We speculate that these processes regulated nutrients levels during anther development. There was no obvious enrichment in these processes in S1 vs. S3 and S2 vs. S3. In contrast, DEGs that were not enriched in S1 vs. S3 and S2 vs. S3 were associated with processes that regulate nutrients levels and were enriched in the phagosome and spliceosome processes. DEGs found to be only enriched in S2 vs. S3 were linked to the N-glycan biosynthesis, mRNA monitoring pathway, protein processing in endoplasmic reticulum, purine metabolism, pyrimidine metabolism, DNA replication, and other processes.
Agronomy 2020, 10, x FOR PEER REVIEW 10 of 19 were linked to the N-glycan biosynthesis, mRNA monitoring pathway, protein processing in endoplasmic reticulum, purine metabolism, pyrimidine metabolism, DNA replication, and other processes.

Analysis of Transcription Factors Regulated by the DEGs
In this study, we analyzed the transcription factors for the DEGs and named them using the transcription factor family information found in the plant transcription factor database. A total of 558 differential TFs were classified into 34 transcription factor families: C2H2, AP2/ERF, bHLH, bZIP, WRKY, MYB, MADS-M-type, etc. ( Figure 6A). A heat map drawn from the gene expression levels of each transcription factor family member ( Figure 6B) shows that AP2/ERF and GRAS were highly expressed in the three stages, while C2H2 and TCP were generally lowly expressed in the three stages. The expression levels of WRKY, MYB, AP2/ERF, and some other transcription families at S1 stage were relatively high, but their expression level patterns showed a downward trend at S2 and S3 stages. At S1 stage, the expression levels of bZIP, MADS-MIKC, and MADS-M-type were relatively low but the level of these transcription factors generally increased in S2 and S3 stages. Notably, PHD and NF-YC family genes were not significantly expressed at stages S1 and S2, but they were upregulated at stage S3. These results indicated that each stage of flower development was associated with different TFs. In addition, it was also found that the number of members in PHD, C2H2, HB-other, AP2/ERF, and C3H in celery accounted for a higher proportion of the corresponding TF family members. The smaller the p-Value, the closer the color is to red and the more significant the enrichment.

Analysis of Transcription Factors Regulated by the DEGs
In this study, we analyzed the transcription factors for the DEGs and named them using the transcription factor family information found in the plant transcription factor database. A total of 558 differential TFs were classified into 34 transcription factor families: C2H2, AP2/ERF, bHLH, bZIP, WRKY, MYB, MADS-M-type, etc. ( Figure 6A). A heat map drawn from the gene expression levels of each transcription factor family member ( Figure 6B) shows that AP2/ERF and GRAS were highly expressed in the three stages, while C2H2 and TCP were generally lowly expressed in the three stages. The expression levels of WRKY, MYB, AP2/ERF, and some other transcription families at S1 stage were relatively high, but their expression level patterns showed a downward trend at S2 and S3 stages. At S1 stage, the expression levels of bZIP, MADS-MIKC, and MADS-M-type were relatively low but the level of these transcription factors generally increased in S2 and S3 stages. Notably, PHD and NF-YC family genes were not significantly expressed at stages S1 and S2, but they were upregulated at stage S3. These results indicated that each stage of flower development was associated with different TFs. In addition, it was also found that the number of members in PHD, C2H2, HB-other, AP2/ERF, and C3H in celery accounted for a higher proportion of the corresponding TF family members.
Notably, the transcriptional abundance and expression of Agr14035 were higher in S1, lower in S2 and S3, while those of Agr25217 were downregulated across the three stages. The expression level of Agr13446 was downregulated by at least four times from S1 to S2, and was upregulated at least two times from S2 to S3.
The expression levels of WRKY family genes (Agr14189, Agr36822, and Agr19214) were significantly different at three stages of flower development. Of note, they were significantly altered in the early flower bud stage, but not significantly altered in between middle flower bud and early flowering stage. The expression levels of Agr14189 decreased by seven-fold from early flower bud to middle flower bud stage, while that of Agr36822 and Agr19214 decreased at least 20 times and 4 times, respectively.
The expression of Agr05322 was similar across the three stages but Agr01062 tended to decline across the three stages. The expression of Agr39039 decreased from S1 to S2 and then increased from S2 to S3, which matched with changes in the corresponding transcriptional level. Agronomy 2020, 10, x FOR PEER REVIEW 12 of 19 Notably, the transcriptional abundance and expression of Agr14035 were higher in S1, lower in S2 and S3, while those of Agr25217 were downregulated across the three stages. The expression level of Agr13446 was downregulated by at least four times from S1 to S2, and was upregulated at least two times from S2 to S3.
The expression levels of WRKY family genes (Agr14189, Agr36822, and Agr19214) were significantly different at three stages of flower development. Of note, they were significantly altered in the early flower bud stage, but not significantly altered in between middle flower bud and early flowering stage. The expression levels of Agr14189 decreased by seven-fold from early flower bud to middle flower bud stage, while that of Agr36822 and Agr19214 decreased at least 20 times and 4 times, respectively.
The expression of Agr05322 was similar across the three stages but Agr01062 tended to decline across the three stages. The expression of Agr39039 decreased from S1 to S2 and then increased from S2 to S3, which matched with changes in the corresponding transcriptional level. The expression level of Agr08149 in the MADS-M-type family was upregulated by seven times in S1 compared to S2, reaching the highest expression level at S2, but it was downregulated by at least five times at S3 compared to S2. The expression level of Agr37458 was upregulated by three times in S1 compared to S2, being highest in S2, and was downregulated by five times in S3 compared to S2. The expression levels of these two genes changed significantly from early flower bud stage to early flowering stage implying that they have important roles in celery flower development.
The expression level of the three genes in the C2H2 transcription family were significantly upregulated from S1 to S2, of which Agr04312 was upregulated by approximately 85-fold, and Agr28404 and Agr30688 were upregulated by at least 50 and 60 times, respectively. At S2 and S3, all of their expression levels were very low and did not alter significantly.

Characteristics of Changes during Anther Development
In angiosperms, flower development is the most important reproductive growth process of plants. The anther is the major pollen-producing part of the plant whose development greatly relies on its fertility. Development of plant anthers is a continuous process related to corresponding cellular processes, such as meiosis, mitosis, the formation of anther walls, development of tapetum cells, and dynamic changes of callose [29,30]. When tapetum cells of the anther develop to a certain stage, they need to be degraded by programmed cell death to provide nutrition and materials needed for the development of microspores. If the degradation of tapetum cells is advanced or delayed, it would lead to male sterility [31]. Previous studies have shown that impaired callose synthesis, accumulation, and degradation would lead to the failure of meiosis and ultimately affect plant fertility [32]. It is of great importance to distinguish the fine structure and features of the anther development period, which can help accurately pinpoint the time and space in which fertility-related features appear. Wheat anthers were divided into 15 stages from stamen formation to anther senescence. The stages included precallose, central callose, meiotic, tetrad, young microspore, vacuolar microspore, among others [33]. The color of flower buds and state of flower dehiscence were observed in celery in comparison with wheat. This can be combined with anther length for cytological analysis. Anther development of Phyllostachys edulis can be classified into six stages: primary sporogenous cell, secondary sporogenous cell, microsporocyte, spore pull over period, mononuclear microspore, and mature pollen [34]. These stages of development are similar to those of celery, but the microspore cells of Phyllostachys edulis only undergo one mitosis, and mature pollen grains are mostly two-celled; while the microspore cells of celery undergo two mitoses, and mature pollen grains are three-celled.

Genes that Regulate Physiological and Biochemical Processes during Flower Development
Anther development is a complex process regulated by multiple genes. Several genes related to anther development were isolated from Arabidopsis, rice, tobacco, and other plant species. Their regulatory functions were identified as well. Studies in Arabidopsis found that TDF1 (DEFECTIVE IN TAPETAL DEVELOPMENT AND FUNCTION 1) and AMS (ABORTED MICROSPORE), which encodes a R2R3 MYB factor and a basic helix-loop-helix (bHLH) factor, respectively, can form a feed-forward loop (FFL) to continuously activate the expression of their regulatory targets [35]. Furthermore, GRF (GROWTH-REGULATING FACTOR) family members such as GRF1, GRF2, GRF3, and GRF5 play a key role in the regulation of carpel and anther development [35,36]. BnaC.CP20.1 plays an important role in tapetum degeneration and pollen wall formation in Arabidopsis [32]. HTH1(HOTHEAD-Like) is involved in keratin biosynthesis, anther development, and pollen fertility of rice [37]. EAT1 (ETERNAL TAPETUM1), belonging to bHLH, is a key regulatory factor that can trigger mitotic phase phasiRNA biogenesis in rice anthers, while NsylCBL10 promoter in tobacco is active in mature anthers [38,39]. In this study, numerous DEGs related to anthers development were screened during flower development by transcriptome sequencing technology, and a regulatory model for the anther development was constructed through functional annotation and expression analysis ( Figure 8). In the GO annotation, early stages of flower development involved ribosomal-related biological processes. Following the completion of anther meiosis, the tapetum would begin to form precursors of sporopollenin, leading to ribosome synthesis of sporopollenin [40]. On the other hand, later stages of flower development involved biological processes of amino sugar metabolism, cytoplasm, threonine-type endopeptidase activity, and threonine-type peptidase activity, which is speculated to be related to the maturation of pollen grains [41]. KEGG annotation also showed differences at different developmental stages. Early stages of flower development largely included oxidative phosphorylation, TCA cycle, carbon fixation in photosynthetic organisms, pentose and glucuronate interconversions, glycolysis/gluconeogenesis, and other metabolic processes. These processes were mainly linked to plant photosynthesis, starch and sucrose metabolism. This is consistent with the evidence that anther development is associated with a decrease in photosynthesis rate and reduced chlorophyll content in anthers [41]. Accordingly, the rate of sucrose and starch metabolism also decreased with flower development [41]. In later stages, DEGs related to anther development were markedly enriched in some processes such as spliceosome, mRNA monitoring, protein processing, purine metabolism, and pyrimidine metabolism. Purines and pyrimidines are raw materials for nucleotide synthesis [42]. Studies have shown that microspores are closely associated with nucleotide metabolism and spliceosome is involved in the formation of pollen walls [43,44]. These results indicate that the maturation of pollen grains occurred at the later stage of flower development. Combined with the ribosomal-related biological processes and cytoplasm in GO annotation, carbon fixation in photosynthetic organisms, and purine and pyrimidine metabolism of KEGG, and the results of anther cytology research, it can be speculated that the early flower bud stage to the middle early flower bud stage may be related to the 3-5 stages of the anther stage, and the middle early flower bud stage to early flowering stage may be associated with the process of 5-9 stages. In summary, anthers of celery have a diverse appearance and diverse internal physiological and biochemical processes at various stages of development.
endopeptidase activity, and threonine-type peptidase activity, which is speculated to be related to the maturation of pollen grains [41]. KEGG annotation also showed differences at different developmental stages. Early stages of flower development largely included oxidative phosphorylation, TCA cycle, carbon fixation in photosynthetic organisms, pentose and glucuronate interconversions, glycolysis/gluconeogenesis, and other metabolic processes. These processes were mainly linked to plant photosynthesis, starch and sucrose metabolism. This is consistent with the evidence that anther development is associated with a decrease in photosynthesis rate and reduced chlorophyll content in anthers [41]. Accordingly, the rate of sucrose and starch metabolism also decreased with flower development [41]. In later stages, DEGs related to anther development were markedly enriched in some processes such as spliceosome, mRNA monitoring, protein processing, purine metabolism, and pyrimidine metabolism. Purines and pyrimidines are raw materials for nucleotide synthesis [42]. Studies have shown that microspores are closely associated with nucleotide metabolism and spliceosome is involved in the formation of pollen walls [43,44]. These results indicate that the maturation of pollen grains occurred at the later stage of flower development. Combined with the ribosomal-related biological processes and cytoplasm in GO annotation, carbon fixation in photosynthetic organisms, and purine and pyrimidine metabolism of KEGG, and the results of anther cytology research, it can be speculated that the early flower bud stage to the middle early flower bud stage may be related to the 3-5 stages of the anther stage, and the middle early flower bud stage to early flowering stage may be associated with the process of 5-9 stages. In summary, anthers of celery have a diverse appearance and diverse internal physiological and biochemical processes at various stages of development.

Transcription Factor Families Show Different Regulation of Flower Development
The regulatory effects of transcription factors at various stages were dissimilar (Figure 8). Of note, bHLH and bZIP were mainly downregulated in the early flower bud stage but upregulated in the early flowering stage. The expression levels of MADS-M-type genes were high in the early stage but low in the late stage. To further analyze the regulatory results of genes, AP2/ERF, WRKY, bZIP, and MADS-M-type were selected for verification based on the abundance of transcription factors and functional annotation results of differential genes. The majority of results obtained by qPCR were consistent with the transcriptome analysis, while discrepancies observed in some results may have been due to differences in sensitivity between the two techniques.
It was reported that transcription factors like AP2/ERF, WRKY, bZIP, and MADS-M-type are involved in the regulation of flower development. AP2 plays an extensive role by regulating the expression of numerous genes in flower organs [45]. For instance, the AP2 gene regulates the formation of perianth in Arabidopsis during the early stages of flower development [46]. In cotton, many GhWRKY genes show high expression levels during anther development, especially in the tetrad pollen (TTP) and uninucleate pollen (UNP) stages [47,48]. In this study, we found that Agr13446, Agr00157, and Agr14035 in the AP2/ERF family and Agr36822, Agr19214, and Agr14189 in the WRKY family had higher expression levels at the S1 stage, suggesting that these genes could play an important role in the early stages of anther development. TGA9 and TGA10, members of bZIP, played an important role in regulating the formation of tapetum and anther dehiscence in Arabidopsis [49]. The expression of Agr39039 of the bZIP family increased then later decreased in the three stages of flower development in celery. Ultimately, there was little change in expression levels in the three stages. It is speculated that Agr39039 plays a regulatory role in the process of anther development. MADS-M-type was a type I MADS-box gene. Previous studies showed that these genes are expressed in male and female gametophytes as well as during endosperm and embryo early development stages. However, they mainly regulate female gametophyte, endosperm, and embryo development [50][51][52]. AGL80, which belongs to MADS-M-type, plays an important regulatory role in the development of central cells and early endosperm, but it has no clear regulatory effect on male gametophytes [53,54]. AGL23 can regulate the formation of female gametophytes and organelles during embryonic development in Arabidopsis [55]. Studies have shown that transcription factors have important regulatory roles in the development of anthers, but the specific regulatory functions of some transcription factors were unclear. For example, MADS-M-type mainly regulates the development process of female gametes, and it had no obvious regulatory effect on male gametes. In this study, some genes in the MADS-M-type family (Agr08149, Agr37458) were highly expressed in the middle flower bud period. They may play an important role in the anther development regulatory network. SUP (SUPERMAN) in Arabidopsis encodes a C2H2 protein, participates in the process of controlling stamen cell proliferation during flower development, and has discrete upstream promoter elements required for early expression in stamen primordium. Studies on the anthers of petunia found that silencing TAZ1 (TAPETUM DEVELOPMENT ZINC FINGER PROTEIN1) leads to abnormal development and premature degeneration of the tapetum, resulting in microspore infertility [56]. In this study, the three genes (Agr04312, Agr28404, Agr30688) from C2H2 family all expressed at the highest level in S1 compared with S2 and S3, indicating that these genes may play a key role in the early stage of anther development in celery [57]. It was reported that the cause of celery male sterility is abnormal morphology of the tapetum before meiosis, and it is mainly due to the immediate degradation of the abnormal tapetum after tetrad stage [21]. The highest expression of C2H2 in celery at S1 may be related to the abnormality of the tapetum and the male sterility of celery, while other transcription factor families like MYB and bHLH are also related to the formation and degradation of the tapetum. AtMYB103 is necessary for the development of tapetum in Arabidopsis, and bHLH transcription factor ETERNAL TAPETUM1 (EAT1) can regulate programmed cell death in the tapetum [38,58]. This study analyzed the development of celery anthers at histological and cytological levels, and linked it to changes in gene expression. The results provide a theoretical basis and excellent gene resources for germplasm improvement of celery and selection of new varieties. As for the genetic regulation mechanism of celery infertility, we could conduct in-depth research combined with the genome sequence of celery.

Conclusions
The purpose of this study was to carry out a cytological study and transcriptome analysis of celery anthers and identify the genes related to anther development. A total of 1074 differentially expressed genes were identified for flower development and could provide a reference for further research on specific genes related to anther development. In particular, transcription factors such as C2H2, AP2/ERF, bZIP, WRKY, and MYB could play a key regulatory role in anther development. In addition, anther development involved physiological and biochemical processes such as ribosome metabolism, glucose metabolism and amino acid metabolism. In summary, these results provide a theoretical basis for molecular breeding of new celery varieties with pollen abortion and theoretical references for reproductive biology research of celery and Apiaceae family.
Supplementary Materials: The following are available online at http://www.mdpi.com/2073-4395/10/5/653/s1, Figure S1: Morphological characteristics of flowers at different developmental stages. Figure S2: GO function annotation of unigenes. Figure S3: KOG function annotation of unigenes. Figure S4: KEGG function annotation of unigenes. Table S1: Statistical results of success rate of gene annotation were obtained from seven databases. Table S2: The comprehensive list of genes with all the annotations from the relevant databases. Table S3: The complist of detailed information for screening of differentially expressed genes of in the three samples.