Transcriptome Analysis of the Development of Pedicel Abscission Zone in Tomato

: Plant organ abscission is a common phenomenon that occurs at a speciﬁc position called the abscission zone (AZ). The differentiation and development of the pedicel AZ play important roles in ﬂower and fruit abscission, which are of great signiﬁcance for abscission in tomatoes before harvest. Previous studies have reported some genes involved in AZ differentiation; however, the genes regulating pedicel AZ cell development in tomatoes after AZ differentiation remain poorly understood. In this study, transcriptome analyses of tomato pedicel AZ samples were performed at 0, 5, 15, and 30 days post-anthesis (DPA). Pedicel AZ growth was mainly observed before 15 DPA. A principal component analysis and a correlation analysis were carried out in order to compare the repeatability and reliability for different samples. We observed 38 up-regulated and 31 down-regulated genes that were signiﬁcantly altered during 0 to 5 DPA, 5 to 15 DPA, and 0 to 15 DPA, which may play key roles in AZ cell enlargement. GO and KEGG enrichment analyses of the selected DEGs under all three different comparisons were conducted. Auxin-signaling-related genes were analyzed, as well as AUX/IAA, GH3, and small auxin up-regulated RNA (SAUR) gene expression patterns. The presented results provide information on pedicel AZ development, which might help in regulating ﬂower or fruit pedicel abscission in tomato production facilities.


Introduction
The tomato (Solanum lycopersicum L.) is a common and popular vegetable due to its rich nutrients, the various ways to eat it, and its pleasing flavor. Flower and fruit abscission before harvest are two of the main factors affecting the tomato yield in production facilities. The location where a flower or fruit separates from the main plant is called the abscission zone (AZ). The AZ comprises layers of functionally specialized cells with morphologically distinct features, such as smaller, square-shaped cells that are interconnected by branched plasmodesmata and dense cytoplasm [1][2][3][4]. The development of a functional AZ is a prerequisite for plant organ abscission. The development of a pedicel AZ is of great significance for flower and fruit pedicel abscission in tomatoes.
The pedicel cells are initially found in the epidermal and cortical region. With the progression of a flower bud's development, the cell division gradually spreads to the vascular bundle region and, finally, to the central parenchyma region [5]. A number of genes were reported to be involved in the differentiation and development of AZs [6][7][8][9][10][11][12][13]. Among these genes, the transcription factor MADS-box plays an important role in the differentiation and development of plant AZs. One of the members of the MADS-box family, jointless, was first isolated as a functional gene involved in the development of plant AZs.
Since tomato plants with the jointless mutation fail to develop AZs on their pedicels, the abscission of flowers or fruit does not occur normally [6]. Another MADS-box, MBP21, has shown a similar function in AZ formation; the loss of function of MBP21 led to a jointless-2 phenotype in tomatoes [8,13]. Macrocalyx (MC) is also a member of the MADS-box family that has been confirmed to play a key role in pedicel AZ development in tomatoes [10]. In Arabidopsis, seedstick (STK) is required for the normal development of the funiculus-an umbilical-cord-like structure that connects the developing seed to the fruit-and for the dispersal of seeds when the fruit matures [7]. Furthermore, blade-on-petiole genes have been verified to play an essential role in abscission formation in Arabidopsis; the loss of function of BOP1 and BOP2 resulted in no differentiation of floral and leaf AZs, suggesting that BOP proteins are essential for the establishment of AZ cells [9,11,12].
Pedicel AZ development includes cell division and cell expansion. After the cells have differentiated into AZ cells, one of the main reasons for pedicel growth is cell expansion. The plant hormone auxin presents various functions in plant growth and development, regulating the fundamental cellular processes of expansion, division, and differentiation [14][15][16]. One of the most striking effects of auxin is rapidly mediating changes in cell expansion [17]. The acid growth theory supposes that auxin activates plasma membrane (PM) H + ATPases, resulting in proton efflux. The decreased pH of the cell wall matrix solution alters the activity of proteins that modify the cell wall, leading to changes in wall extensibility. Furthermore, elevated PM H + ATPase activity hyperpolarizes the PM and increases the energy required for solute uptake (which is necessary for the maintenance of water uptake and, therefore, the turgor pressure), forcing the expansion of the wall [18][19][20]. Auxin-regulated cell expansion has also been found to play a crucial role in leaf epinasty and storage organ expansion [21,22]. SAUR, which is a member of the auxin signaling pathway, has been reported to be involved in cell expansion through the regulation of PM H + ATPase activity [23,24]. Transmembrane kinase (TMK) auxin-signaling proteins interact with plasma membrane H + ATPases, inducing their phosphorylation, which is required for auxin-induced H + -ATPase activation, apoplastic acidification, and cell expansion [25]. All of the above results indicate that auxin plays a key role in plant cell wall expansion. In addition to auxin, other hormones are also involved in cell expansion, cell elongation, and cell division. Brassinosteroids (BR) were reported to stimulate the pollen tube growth and hypocotyl elongation of pak choi (Brassica rapa chinensis) through cell and cell wall expansion [26][27][28]. Cytokinin regulated cell proliferation by influencing cell division and/or differentiation [29]. There was shown to be crosstalk between gibberellin (GA) and auxin in cell elongation; auxin promoted the degradation of Della in root cells in response to GA, which is a prerequisite for GA-induced root elongation [30].
For this study, we performed transcriptome analyses of pedicel AZs in tomatoes and aimed to determine the genes related to the development of the pedicel AZ. The results provide further information regarding the regulation of flower and fruit abscissions in tomato production.

Plant Materials and Sample Preparation
Tomato plants (Solanum lycopersicum L. cv. "micro-Tom") were grown in a greenhouse at 25/18 • C with a 16/8 h light cycle (day/night). To each flower, a small tag-which was labeled with the flower open date (i.e., the beginning of the anthesis stage)-was attached at 10:00 a.m. every day. Samples were collected at 0, 5, 15, and 30 days post-anthesis (DPA) from different plants that grew uniformly. The plants were grown at the same time and under the same conditions. At each stage, the pedicel abscission zones were cut into small segments (about 3 mm) using a sharp blade, and they were then immediately frozen in liquid nitrogen. The segments of pedicel AZs were weighed with a balance. Each sample weighted at least 1 g. Three replicates of each stage were considered.

AZ Diameter Measurement
The diameter of the AZ was measured using a vernier caliper at 10:30 a.m. each day.
At least 30 AZs were tested at each time point. Each measurement was recorded in an Excel spreadsheet. The statistical significance was determined using one-way ANOVA and is indicated by asterisks (*) in the diagram.

mRNA Sequencing by Illumina HiSeq
The total RNA of each sample was extracted using TRIzol Reagent (Invitrogen), an RNeasy Mini Kit (Qiagen). The total RNA of each sample was quantified and checked for quality using an Agilent 2100 Bioanalyzer (Agilent Technologies, Palo Alto, CA, USA), NanoDrop (Thermo Fisher Scientific Inc, Waltham, MA, USA) and 1% agarose gel. For the libraries' construction, please refer to the Illumina sequencing method [31]. The sequences were processed on an Illumina HiSeq™ 2500 platform and analyzed by GENEWIZ (GE-NEWIZ, Inc. 115 Corporate Boulevard, South Plainfield, NJ, USA). The quality control was processed using Trimmomatic (v0.30). The sequencing error rate distribution, GC content distribution, and sequencing data filtering were done successively. In order to obtain the clean data, the adapter reads, bases containing N, and low-quality reads were sequentially removed. For reads mapping, the reference genome sequences and gene model annotation files of relative species were first downloaded from the genome website (https://solgenomics.net, accessed on 15 April 2019) using Current Tomato Genome version SL4.0 and Annotation ITAG4.0. Next, Hisat2 (v2.0.1) was used to index the reference genome sequence. Finally, the clean data were aligned to reference the genome via the software Hisat2 (v2.0.1) [32].

Differential Expression Analysis
For the differential expression analysis, we used the DESeq Bioconductor package-a model based on the negative binomial distribution. After an adjustment using the approach of Benjamini and Hochberg for controlling the false discovery rate, the P-values for the genes were determined [33,34].

Principal Component Analysis
A principal component analysis provides a way to visualize sample-to-sample distances. In this ordination method, the data points (here, the samples) are projected onto a 2D plane and spread out in the two directions that explain most of the variation in the data; the x-axis is the direction that separates the data points the greatest amount (the values of the samples in this direction are in terms of PC1), while the y-axis denotes the direction (which must be orthogonal to the first direction) that separates the data the second-greatest amount (the values of the samples in this direction are in terms of PC2). The percentage of the total variance that is contained in each direction is printed on the axis label. Note that these percentages generally do not add to 100%, as there are more dimensions that explain the remaining variance (although each of these remaining dimensions will explain fewer than the two that are presented). This analysis was performed using the OmicShare tool, a free online platform for data analysis (https://www.omicshare.com/tools, accessed on 8 June 2022).

Correlation Analysis
The correlation of each of the two samples was calculated using the Pearson method. In this method, the correlation values range between −1 and 1. If the correlation value is between −1 and 0, the two samples have a negative correlation. On the contrary, if the correlation value is between 0 and 1, the correlation between the two samples is positive. If the correlation value is 0, there is no correlation between the two samples. The closer that the absolute value of the correlation coefficient is to 1, the stronger that the correlation is between the two variables; meanwhile, the closer that the correlation coefficient is to zero, the weaker the correlation is between the two variates. This analysis was also performed using the OmicShare tool (https://www.omicshare.com/tools, accessed on 10 June 2022).

GO and KEGG Enrichment Analysis
GO Term Finder was used to identify gene ontology (GO) terms in order to annotate a list of enriched genes with significant P-values (i.e., less than 0.05). KEGG (Kyoto Encyclopedia of Genes and Genomes) is a collection of databases dealing with genomes, biological pathways, diseases, drugs, and chemical substances (http://en.wikipedia.org/ wiki/KEGG, accessed on 13 June 2022). We used in-house scripts to enrich significantly differentially expressed genes in KEGG pathways [31].

Quantitative RT-PCR
Tomato pedicel AZs at 0 DPA, 5 DPA, 15 DPA and 30 DPA were collected for total RNA isolation. The total RNAs were extracted using the TRIzol reagent, following the protocol provided by the manufacturer (Invitrogen), and treated with DNase I (Thermo Scientific). The total RNAs from each sample were used for first-strand cDNA synthesis in a final volume of 20 µL. For real-time quantitative PCR (qRT-PCR), a method with two steps of Takara real-time PCR was used to amplify the representative auxin-related genes. Each experiment was performed independently three times (technical replicates) with three biological samples. Each biological sample's Ct value was the average of the three technical replicates. Each sample's Ct value was the average of the three biological Ct values. The data analysis for qRT-PCR was performed following the ∆∆Ct model [35]. The relative expression of qRT-PCR was calculated according to the formula 2 [Ct(actin) − Ct(gene)] . The error bar was the standard deviation of the three biological replicates. The primers used in this study are listed in Supplemental File Table S7. The raw Ct values are provided in Supplemental File Table S8.

Development of Pedicel AZ in Tomato
The different development stages of tomato pedicels are shown in Figure 1a. The pedicel AZ size was measured at 0 (the day of flower opening), 5, 15, and 30 days postanthesis (DPA). About 30 pedicels of tomato plants were measured using Vernier calipers. The diagram shows that the pedicel AZs significantly changed between the day of flower opening and 15 DPA (Figure 1b). The diameter of the AZ was about 1.2, 2.1, 3.2, and 3.3 mm at 0, 5, 15, and 30 DPA, respectively. There were significant differences between 0 and 5 DPA and 5 and 15 DPA, but there was no significant difference between 15 and 30 DPA. The results indicate that the main period of pedicel AZ growth or size expansion is from 0 to 15 DPA, and the AZ diameter almost reached its peak at 15 DPA.

Transcriptome Analysis of Pedicel AZ at Different Development Stages in Tomatoes
The tomato pedicel AZs of 0, 5, 15, and 30 DPA were collected for transcriptome sequencing. The principal component analysis showed that the distance between the three replicates of each group was relatively small, which meant that the interclass replicates were reliable ( Figure 1c). The samples at 0, 5, and 15 DPA showed large distances between the groups; however, the samples at 15 and 30 DPA were very close to each other ( Figure 1c). This result indicates that the pedicel AZs at 15 DPA were similar to those at 30 DPA. A correlation analysis of the 12 samples also showed a similar result (Figure 1d). The correlation coefficients of the interclass replicates were all very close to one, which meant that strong correlations existed in the two samples. The correlation coefficients between the samples at 15 and 30 DPA were also very close to one (Figure 1d

Differentially Expressed Genes That Were Selected under Three Comparisons
With the main growth stage separated into three parts, as above, we wished to determine the DEGs that had all changed under the three comparisons.  (Figure 3b). Genes were discarded if their absolute values of fold change were less than five. For the up-regulated genes, if a gene's fpkm value was less than 1 at 15 DPA, it was abandoned; however, for the down-regulated genes, if a gene's fpkm value was less than 1 at 0 DPA, it was also abandoned, as it may barely play a role during the development of the pedicel AZ in tomatoes. Based on these criteria, we selected 38 up-regulated (Table 1) and 31 down-regulated (Table 2) genes. The ID, fpkm values, and descriptions of these genes are provided in Tables 1 and 2. genes, 125 DEGs were changed under all three comparisons, while 1504, 1339, and 389 DEGs were only changed for 0-5 DPA, 5-15 DPA, and 0-15 DPA, respectively (Figure 3a). For the down-regulated genes, 167 DEGs were changed under all three comparisons, while 1103, 1710, and 422 DEGs were only changed for 0-5 DPA, 5-15 DPA, and 0-15 DPA, respectively (Figure 3b). Genes were discarded if their absolute values of fold change were less than five. For the up-regulated genes, if a gene's fpkm value was less than 1 at 15 DPA, it was abandoned; however, for the down-regulated genes, if a gene's fpkm value was less than 1 at 0 DPA, it was also abandoned, as it may barely play a role during the development of the pedicel AZ in tomatoes. Based on these criteria, we selected 38 up-regulated (Table 1) and 31 down-regulated (Table 2) genes. The ID, fpkm values, and descriptions of these genes are provided in Tables 1 and 2.

GO and KEGG Enrichment Analysis of the Key DEGs under Three Comparisons
The biological functions of the selected DEGs were annotated with GO terms under all three different developmental stages. For the up-regulated genes, 14 and 7 GO terms were in the categories associated with molecular functions and biological processes, respectively. The top three enriched GO terms in the molecular function category were "oxidoreductase activity", "DNA binding", and "protein binding", while those in the biological process category were "oxidation-reduction process", "lipid metabolic process", and "metabolic process" (Figure 4a). For the down-regulated genes, 3, 16, and 8 GO terms were part of the cellular component, molecular function, and biological process categories, respectively. The top three enriched GO terms in the cellular component category were "membrane", "extracellular region", and "integral component of membrane"; those in the molecular function category were "protein dimerization activity", "oxidoreductase activity", and "strictosidine synthase activity"; and those in the biological process category were "oxidation-reduction process", "transmembrane transport", and "biosynthetic process" (Figure 4b). The information of the DEGs that were annotated in each of the GO terms was showed in Table S9. The KEGG enrichment analysis showed that the DEGs under all three comparisons were enriched in 14 KEGG pathways (Figure 4c). The up-regulated genes were involved in nine KEGG pathways, where the top three pathways were "phenylpropanoid biosynthe-sis", "metabolic pathways", and "biosynthesis of secondary metabolites". Meanwhile, the down-regulated genes were involved in five KEGG pathways, including "amino sugar and nucleotide sugar metabolism", "diterpenoid biosynthesis", "fatty acid elongation", "metabolic pathways", and "biosynthesis of secondary metabolites" (Figure 4c).

Analysis of DEGs Involved in the Auxin Signaling Pathway
Plant hormones are well-known to be involved in plant growth development. An enrichment analysis for all hormone-related genes involved in pedicel AZ development is shown in Figure S1. Among the hormone pathways, the DEGs were significantly enriched in the auxin signaling pathway. Auxin is one of the main factors affecting cell expansion. The auxin signaling pathway is depicted in Figure 5a. A number of auxin-related genes were identified in the three different comparisons. The enrichment analysis of the DEGs belonging to different classes of auxin response genes was showed in figure S2. The numbers of DEGs belonging to AUX/IAA (K14484) and SAUR (K14488) were higher than those related to the other three families (Figures 5b and S1). In the AUX1 (13946) family, there were no up-regulated genes in any of the three stages. In the ARF (14486) family, only two genes were up-regulated during 0-5 DPA, while one gene was down-regulated during 5-15 DPA and 0-15 DPA. In the GH3 (14487) family, three and four genes were up-regulated during 0-5 DPA and 0-15 DPA, respectively, while one gene was downregulated during 5-15 DPA (Figure 5b).

Analysis of DEGs Involved in the Auxin Signaling Pathway
Plant hormones are well-known to be involved in plant growth development. An enrichment analysis for all hormone-related genes involved in pedicel AZ development is shown in Figure S1. Among the hormone pathways, the DEGs were significantly enriched in the auxin signaling pathway. Auxin is one of the main factors affecting cell expansion. The auxin signaling pathway is depicted in Figure 5a. A number of auxin-related genes were identified in the three different comparisons. The enrichment analysis of the DEGs belonging to different classes of auxin response genes was showed in Figure S2. The numbers of DEGs belonging to AUX/IAA (K14484) and SAUR (K14488) were higher than those related to the other three families (Figures 5b and S1). In the AUX1 (13946) family, there were no up-regulated genes in any of the three stages. In the ARF (14486) family, only two genes were up-regulated during 0-5 DPA, while one gene was down-regulated during

Expression Pattern of Auxin-Related Genes during All Three Development Stages
The number of auxin-related genes significantly differed in the three different stages, and the expression patterns of these genes are shown in Figure 6. In the four AUX/IAA family members, the expression of three genes-Solyc12g096980 (IAA11), Solyc06g066020 (IAA36), and Solyc03g120380 (IAA19)-showed an increasing and then decreasing trend, with a peak at 5 DPA. The other gene, Solyc08g021820 (IAA29), also showed an increasing and then decreasing trend, but its highest expression was at 15 DPA (Figure 6a). The only GH3 family member, Solyc02g092820, also showed an increasing and then decreasing trend, with an expression pattern very similar to that of IAA11 during the pedicle AZ development in tomatoes (Figure 6b). In the SAUR family, the expression of two genes-Solyc08g079150 and Solyc01g110580-showed a similar trend, with both showing their highest expression at 0 DPA. The other gene, Solyc06g053260, showed a trend opposite to these two genes, with a peak at 30 DPA (Figure 6c). The expression pattern of these genes

Expression Pattern of Auxin-Related Genes during All Three Development Stages
The number of auxin-related genes significantly differed in the three different stages, and the expression patterns of these genes are shown in Figure 6. In the four AUX/IAA family members, the expression of three genes-Solyc12g096980 (IAA11), Solyc06g066020 (IAA36), and Solyc03g120380 (IAA19)-showed an increasing and then decreasing trend, with a peak at 5 DPA. The other gene, Solyc08g021820 (IAA29), also showed an increasing and then decreasing trend, but its highest expression was at 15 DPA (Figure 6a). The only GH3 family member, Solyc02g092820, also showed an increasing and then decreasing trend, with an expression pattern very similar to that of IAA11 during the pedicle AZ development in tomatoes (Figure 6b). In the SAUR family, the expression of two genes-Solyc08g079150 and Solyc01g110580-showed a similar trend, with both showing their highest expression at 0 DPA. The other gene, Solyc06g053260, showed a trend opposite to these two genes, with a peak at 30 DPA (Figure 6c). The expression pattern of these genes were verified by qRT-PCR (Figure 7). All of these genes showed a similar expression pattern with the RNA-Seq data ( Figure 6). Two AUX/IAA genes-Solyc12g096980 (IAA11) and Solyc06g066020 (IAA36) (Figure 7a)-and the GH3 member-Solyc02g092820 (Figure 7b)-showed an increasing and then decreasing trend, but the peak of these genes' expression were at 15DPA instead of the 5DPA shown by the RNA-Seq data. were verified by qRT-PCR ( Figure 7). All of these genes showed a similar expression pattern with the RNA-Seq data ( Figure 6). Two AUX/IAA genes-Solyc12g096980 (IAA11) and Solyc06g066020 (IAA36) (Figure 7a)-and the GH3 member-Solyc02g092820 ( Figure  7b)-showed an increasing and then decreasing trend, but the peak of these genes' expression were at 15DPA instead of the 5DPA shown by the RNA-Seq data.

Discussion
Plant cells are surrounded by cell walls, which provide structural integrity but also spatially constrain cells and, therefore, must be modified to allow for cellular expansion [16], which is triggered by a high intracellular turgor pressure. The wall properties regulate the differential growth of the cell, resulting in a diversity of cell sizes and shapes. The acid growth theory supposes that a decrease in the pH of the cell wall matrix's solution

Discussion
Plant cells are surrounded by cell walls, which provide structural integrity but also spatially constrain cells and, therefore, must be modified to allow for cellular expansion [16], which is triggered by a high intracellular turgor pressure. The wall properties regulate the differential growth of the cell, resulting in a diversity of cell sizes and shapes. The acid growth theory supposes that a decrease in the pH of the cell wall matrix's solution alters the wall's extensibility [18][19][20]. However, different cells have displayed distinct abilities to perceive acid; mature cells are less sensitive to acidic pH and extend less than young cells [18,36]. In this study, the enlargement of pedicel AZ size was mainly observed from 0 to 15 DPA (Figure 1a), and the growth of the pedicel AZ was slow after 15 DPA, indicating that the cells of the pedicel AZ at 15 DPA might be almost mature. The cells of the pedicel AZ after 15 DPA were less sensitive to acidic pH conditions, which led to cell expansion or to an organ's growth slowing or stopping. A correlation analysis of samples at different developmental stages confirmed this result: the coefficients of the correlation between the samples at 15 and 30 DPA were very close to one, which meant the samples at these two stages were very similar to each other (Figure 1c). Thus, our results indirectly verify the supposition that mature cells are less sensitive to acidic pH than young cells are.
It was interesting to find that seven pathways-"photosynthesis", "alpha-linolenic", "glycosphingolipid biosynthesis-lacto and neolacto series", "glyco-sphingolipid biosynthesisglobo and isoglobo series", "glyoxylate and dicarboxylate metabolism", "brassinosteroid biosynthesis", and "carotenoid biosynthesis"-were not significantly changed during both 0-5 DPA and 5-15 DPA, whereas they were significantly altered during 0-15 DPA (Figure 2). This suggested that the genes involved in these pathways might show the same trend between 0 and 5 DPA and between 5 and 15 DPA, such as when they were on the rise. The pathway increased both a little bit from 0 to 5 DPA and 5 to 15 DPA, but these two increases resulted in a significant rise overall from 0 to 15 DPA. Conversely, for the downward trend, the pathways might show similar results. Two pathways-"amino sugar and nucleotide sugar metabolism" and "plant-pathogen interaction"-were changed during both 0-5 DPA and 5-15 DPA, but not in 0-15 DPA (Figure 2). This result indicated that the two pathways might show a reverse trend between 0 and 5 DPA and between 5 and 15 DPA. This meant that when the pathway first rose and then fall or when it first fell and then rose, the high peak or the low peak was at 5 DPA but showed similar levels at 0 and 15 DPA. Based on this finding, the two pathways might play important roles on the early stage of pedicel AZ development in tomatoes.
In addition to these three core components, the genes belonging to the AUX1, GH3, and SAUR families are also involved in the auxin signaling pathway (Figure 5a). Auxin is known to induce acid growth, which is regulated by the auxin-inducible SAUR proteins [24]. In this study, the number of DEGs belonging to the SAUR family involved in the auxin signaling pathway was the largest (Figure 5b). This suggests that SAUR genes might play an important role during pedicel cell development, which is consistent with the results of previous studies.
In conclusion, we analyzed gene expressions at the transcriptional level and selected some differentially expressed genes that may play an important role in pedicel AZ development. This study provides a foundation for research into plant organ abscission, as well as providing a reference for plant organ growth and cell enlargement research.