Transcriptome and Metabolome Analyses Provide Insights into the Stomium Degeneration Mechanism in Lily

Lily (Lilium spp.) is a widely cultivated horticultural crop that has high ornamental and commercial value but also the serious problem of pollen pollution. However, mechanisms of anther dehiscence in lily remain largely unknown. In this study, the morphological characteristics of the stomium zone (SZ) from different developmental stages of ‘Siberia’ lily anthers were investigated. In addition, transcriptomic and metabolomic data were analyzed to identify the differentially expressed genes (DEGs) and secondary metabolites involved in stomium degeneration. According to morphological observations, SZ lysis occurred when flower buds were 6–8 cm in length and was completed in 9 cm. Transcriptomic analysis identified the genes involved in SZ degeneration, including those associated with hormone signal transduction, cell structure, reactive oxygen species (ROS), and transcription factors. A weighted co-expression network showed strong correlations between transcription factors. In addition, TUNEL (TdT-mediated dUTP nick-end labeling) assays showed that programmed cell death was important during anther SZ degeneration. Jasmonates might also have key roles in anther dehiscence by affecting the expression of the genes involved in pectin lysis, water transport, and cysteine protease. Collectively, the results of this study improve our understanding of anther dehiscence in lily and provide a data platform from which the molecular mechanisms of SZ degeneration can be revealed.


Introduction
Lily (Lilium spp.) is highly valued for its attractive variations in fragrance, pattern, and color, and holds a high market share in the global cut flower industry [1,2]. Control of male fertility is essential for plant reproduction and selective breeding [3]. On the other hand, pollen can also have serious negative effects. For example, many ornamental plants produce excessive amounts of pollen, including Populus, lily, and chrysanthemum, and exposure to high amounts of pollen can lead to an anaphylactic reaction in susceptible populations [4][5][6]. Moreover, large amounts of pollen are not necessary for the commercial value of ornamental plants. Therefore, sterile male or pollen-free ornamental plants can be very valuable. Although lily is an important commercial flower crop, its market popularity and landscape application are adversely affected by severe pollen pollution. Many studies have focused on anther development in model plants, but there has been little research on flower crops, and although pollen development in lily has been examined in recent years [5,7], the mechanisms of anther dehiscence remain largely unknown.
In plant evolution, selfing leads to genetic decline in offspring. Differences in stamen fertility can be an important factor in heredity and evolution. In Arabidopsis, late anther Int. J. Mol. Sci. 2021, 22, 12124 2 of 21 sterility involves anther dehiscence, pollen maturation, and filament elongation [8]. Anther dehiscence and pollen maturation occur after the anther middle layer is crushed and meiosis is completed [9]. Anther dehiscence is a multistage process that involves differentiation of specialized cells and changes in water status [3].
Plant hormones also affect dehiscence. The Arabidopsis arf17 mutant has defective endothecium thickening and male sterility, and ARF17 can directly regulate MYB108, which also has functions in anther endothecium thickening [15]. Similarly, anthers of Arabidopsis tir1 afb2 afb3 triple and tir1 afb1 afb2 afb3 quadruple mutants show early endothecium thickening [16]. At stage 10 of anther development, lignification, which leads to anther cracking in advance, occurs in the endothecium of afb1 opr3 anthers but not in the endothecium of wild-type plants [8]. Studies using dehiscent mutants show that JAs (JAs) also contribute to control anther dehiscence [3,17,18]. Jasmonate biosynthetic enzyme mutants, including the fatty acid desaturation (fad) mutant, opr3 (mutation in 12-oxophytodienoic acid reductase) mutant, delayed-dehiscence 1 (dde1) and dde2 mutant, defective in anther dehiscence 1 (dad1) mutant, and allene oxide synthase mutant [3], indicate JA involvement in anther indehiscence. The lox3 lox4 double mutant has abnormal anther maturation and defective dehiscence [17]. In addition, delayed dehiscence has also been observed in mutants defective in JA perception. The JA receptor COI 1 (coronatine insensitive 1) is involved in anther fertility [18]. As a JA signal transporter, MYB21/24 interacts with both MYC and JAZ proteins, which can repress the bHLH-MYB complex activation function. However, JA induces JAZ degradation and releases the bHLH-MYB complex to subsequently activate expression of downstream genes for anther dehiscence development [19,20].
Mechanisms regulating anther dehiscence have been studied in model plants, but the focus has largely been limited to endothecium cell lignification [3,[10][11][12][13][14][15]. Few studies have examined stomium zone (SZ) degeneration. Molecular regulation mechanisms of anther SZ degeneration remain largely unknown. In this study, SZ degeneration was investigated in the large anthers of lily using a combination of morphological, physiological, and comparative transcriptomic and metabolomic analyses. The aim of the research was to understand SZ degeneration at different stages of anther development at the morphological, physiological, and molecular levels. Transcriptomic analysis of the five developmental stages of SZ tissues was conducted to identify the differentially expressed genes (DEGs) associated with SZ lysis, and metabolomic analysis was conducted to explore differences in metabolism during SZ degeneration. The results will provide new insights into the regulation of SZ lysis and help to screen new factors to control anther dehiscence.

Occurrence of Anther Dehiscence in Lily Flower Buds 6-8 cm in Length
The Oriental hybrid lily 'Siberia' was used to investigate the relation between SZ degradation and flower development. Developmental periods of lily anthers were separated according to flower bud length. Stages of flower buds at 2-9 cm in length were designated as S2-S9 cm. In stages S2-S5 cm, SZ degeneration was not observed in anthers (Figure 1a-d). SZ lysis occurred during stages S6-S8 cm (Figure 1e-g) and was completed at stage S9 cm (Figure 1h). In addition, to further observe the morphology of developing flower buds, stamens, pistils, and perianths, the morphology of the whole flower during stages S6-S9 cm is shown in Figure 1i-t. Pigment accumulation in the SZ was observed in stages S6-S9 cm (Figure 1i-l). These results indicated that stages S6-S8 cm were the critical ones in SZ lysis.

Regulation of Stomium Zone Degradation Showing Fewer Differentially Expressed Genes
The SZ samples were collected by cutting anthers freehand, which were then immediately frozen in liquid nitrogen and stored at −80 °C for transcriptome sequencing. Based on FPKM (fragments per kb per million reads), the gene expression levels were explored in the anther SZ during five developmental stages (S5-S9 cm) using RNA-seq. More than 755.56 million clean reads were generated from 15 samples. The Q30 rates of the 15 libraries ranged from 89.24% to 94.90%. The GC rates of the 15 libraries ranged from 49.02% to 50.22%. Results are shown in Table S1.
Overlap of DEGs is shown in a Venn diagram (Figure 2b). Among the DEGs, only 189 were common to the five stages, which indicated that there were different response mechanisms involved in the regulation of anther dehiscence during SZ degradation.

Regulation of Stomium Zone Degradation Showing Fewer Differentially Expressed Genes
The SZ samples were collected by cutting anthers freehand, which were then immediately frozen in liquid nitrogen and stored at −80 • C for transcriptome sequencing. Based on FPKM (fragments per kb per million reads), the gene expression levels were explored in the anther SZ during five developmental stages (S5-S9 cm) using RNA-seq. More than 755.56 million clean reads were generated from 15 samples. The Q30 rates of the 15 libraries ranged from 89.24% to 94.90%. The GC rates of the 15 libraries ranged from 49.02% to 50.22%. Results are shown in Table S1.
To further explore the function of the 18,012 assembled genes, GO ( Figure 3a) term enrichment analysis was performed, which has three categories: biological process, cellular component, and molecular function. In the biological process category, "cellular process" was the most highly represented term, followed by "metabolic process" and "singleorganism process". In the cellular component category, the main functional terms were "cell", "cell part", and "organelle". In the molecular function category, the most enriched terms were "binding", "catalytic activity", and "transporter activity" (Figure 3a).

Jasmonate Biosynthesis and Signal Transduction Roles in Stomium Zone Degradation
To identify key genes involved in SZ degeneration in lily anthers, expression patterns of phytohormone-related DEGs (S9 cm-vs-S5 cm) were analyzed based on their function. The two largest classifications were auxin and JA signal transduction (Table S3). Auxin signal transduction includes many members, such as auxin influx carrier protein AUX1, auxin-responsive protein IAA, auxin responsive GH3 family proteins, and SAUR family proteins.

Genes Associated with Cysteine Proteinases during Stomium Zone Degradation
Genes associated with pectin degeneration, water and sucrose transport, and programmed cell death (PCD) were differentially expressed during SZ degradation (Table  S3). Pectin metabolic genes involved in cell wall degeneration were identified, including exopolygalacturonase, polygalacturonase, pectinesterase, hypotheticalprotein, pectinesterase-like, pectate lyase, pectate lyase-like, and pectinesterase inhibitor. The 25 genes associated with pectin degeneration were mostly expressed in stage S9 cm ( Figure 5). Water and sucrose transport genes included those of several aquaporins and sucrose transport proteins. Six aquaporins were primarily expressed in stage S5 cm, but sucrose transport protein genes were highly expressed in stage S8 cm. Five PCD-related genes were identified, including athepsin B-like protease, cysteine proteinase, and those of other proteases. High expression levels of PCD-related genes were primarily in stages S6 cm and S7 cm. This result suggested that the 15 protease genes significantly affected SZ cell death.  Table S3.

Genes Associated with Cysteine Proteinases during Stomium Zone Degradation
Genes associated with pectin degeneration, water and sucrose transport, and programmed cell death (PCD) were differentially expressed during SZ degradation (Table S3). Pectin metabolic genes involved in cell wall degeneration were identified, including exopolygalacturonase, polygalacturonase, pectinesterase, hypotheticalprotein, pectinesterase-like, pectate lyase, pectate lyase-like, and pectinesterase inhibitor. The 25 genes associated with pectin degeneration were mostly expressed in stage S9 cm ( Figure 5). Water and sucrose transport genes included those of several aquaporins and sucrose transport proteins. Six aquaporins were primarily expressed in stage S5 cm, but sucrose transport protein genes were highly expressed in stage S8 cm. Five PCD-related genes were identified, including athepsin B-like protease, cysteine proteinase, and those of other proteases. High expression levels of PCD-related genes were primarily in stages S6 cm and S7 cm. This result suggested that the 15 protease genes significantly affected SZ cell death.

Figure 5.
Heat maps of the differentially expressed genes associated with pectinase, water and sucrose transport, and cysteine proteinases.

Transcription Factors Roles in Anther Dehiscence
Many studies found that anther dehiscence is governed by special transcription factors (TFs), which are the master regulators controlling how anther development-related genes are turned off and on [10,13,15]. The 49 differentially expressed TFs in the group S8 cm-vs-S5 cm were examined (Figure 6a and Table S4). In the MYB TF family, four TFs were highly expressed in stage S8 cm, and one TF (TRINITY_DN15167_c0_g1_i3_1_2) was highly expressed in stage S6 cm. In the MADS TF family, TFs were highly expressed from stages S6 cm to S9 cm. In the NAC TF family, three TFs were highly expressed in stage S5 cm, two TFs were highly expressed in stage S6 cm, and one TF was highly expressed in stage S8 cm. In the TIFY TF family, three TFs were most highly expressed in stages S5 cm, S6 cm, and S8 cm. Among five bHLH TFs, three TFs were most upregulated in stage S6 cm, whereas the other two TFs had maximal FPKM values in stages S5 cm and S8 cm. Two AP2/ERF TFs had maximal FPKM values in stage S6 cm. Three C3H family members had maximal FPKM values in stages S5 cm, S6 cm, and S9 cm. Three bZIP family TFs had maximal FPKM values in stages S5 cm, S6 cm, and S8 cm. Two GRAS family TFs were highly expressed in S8 cm. Eight TFs from the HB-HD-ZIP and OFP families had significantly higher expression in stages S5 cm, S6 cm, and S8 cm. Eight other TFs in LIM, TCP, PLATZ, NF-YB, B3-ARF, GARP, C2H2, and LOB families were also identified. The highest Figure 5. Heat maps of the differentially expressed genes associated with pectinase, water and sucrose transport, and cysteine proteinases.

Transcription Factors Roles in Anther Dehiscence
Many studies found that anther dehiscence is governed by special transcription factors (TFs), which are the master regulators controlling how anther development-related genes are turned off and on [10,13,15]. The 49 differentially expressed TFs in the group S8 cmvs-S5 cm were examined (Figure 6a and Table S4). In the MYB TF family, four TFs were highly expressed in stage S8 cm, and one TF (TRINITY_DN15167_c0_g1_i3_1_2) was highly expressed in stage S6 cm. In the MADS TF family, TFs were highly expressed from stages S6 cm to S9 cm. In the NAC TF family, three TFs were highly expressed in stage S5 cm, two TFs were highly expressed in stage S6 cm, and one TF was highly expressed in stage S8 cm. In the TIFY TF family, three TFs were most highly expressed in stages S5 cm, S6 cm, and S8 cm. Among five bHLH TFs, three TFs were most upregulated in stage S6 cm, whereas the other two TFs had maximal FPKM values in stages S5 cm and S8 cm. Two AP2/ERF TFs had maximal FPKM values in stage S6 cm. Three C3H family members had maximal FPKM values in stages S5 cm, S6 cm, and S9 cm. Three bZIP family TFs had maximal FPKM values in stages S5 cm, S6 cm, and S8 cm. Two GRAS family TFs were highly expressed in S8 cm. Eight TFs from the HB-HD-ZIP and OFP families had significantly higher expression in stages S5 cm, S6 cm, and S8 cm. Eight other TFs in LIM, TCP, PLATZ, NF-YB, B3-ARF, GARP, C2H2, and LOB families were also identified. The highest FPKM values of those TFs were in stages S5 cm to S9 cm. They were initially upregulated during anther dehiscence. The results suggested TF genes regulated SZ degeneration.
Based on the expression levels, correlations between the 49 TFs were analyzed ( Figure S2). Among the TFs, there were more significant positive correlations than negative correlations. A weighted network of the TFs was constructed on the basis of their expression patterns. Forty-nine nodes and 1176 co-expressed gene pairs were identified. Because the gene pairs with a p-value < 0.05 were considered to show consistent correlations, they were then used in the weighted network construction (Figure 6b). As shown in Figure 6, there were many interconnections among the TFs, which indicated that SZ degeneration was a coordinated process. Among the TFs, ethyleneresponsive transcription factor WIN1 (TRINITY_DN10983_c0_g1_i2_2_2), transcription factor SPATULA (TRINITY_DN15114_c0_g2_i5_2_2), transcription factor MYB86 (TRIN-ITY_DN15167_c0_g1_i3_1_2), homeobox-leucine zipper protein HOX32 (TRINITY_DN25842 _c1_g1_i2_1_2), transcription repressor OFP13 (TRINITY_DN25528_c0_g1_i2_1_2), NAC domain-containing protein 83 (TRINITY_DN19575_c0_g1_i2_1_2), NAC domain-containing protein 43 (TRINITY_DN22402_c0_g1_i2_2_2), and bZIP transcription factor RISBZ5 (TRIN-ITY_DN22982_c0_g1_i3_1_2) had higher connectivity than that of the other TFs. FPKM values of those TFs were in stages S5 cm to S9 cm. They were initially upregulated during anther dehiscence. The results suggested TF genes regulated SZ degeneration. Based on the expression levels, correlations between the 49 TFs were analyzed (Figure S2). Among the TFs, there were more significant positive correlations than negative correlations. A weighted network of the TFs was constructed on the basis of their expression patterns. Forty-nine nodes and 1176 co-expressed gene pairs were identified. Because the gene pairs with a p-value < 0.05 were considered to show consistent correlations, they were then used in the weighted network construction (Figure 6b). As shown in Figure 6, there were many interconnections among the TFs, which indicated that SZ degeneration was a coordinated process. Among the TFs, ethylene-responsive transcription factor WIN1 (TRINITY_DN10983_c0_g1_i2_2_2), transcription factor SPATULA (TRIN- NAC domain-containing protein 43 (TRIN-ITY_DN22402_c0_g1_i2_2_2), and bZIP transcription factor RISBZ5 (TRIN-ITY_DN22982_c0_g1_i3_1_2) had higher connectivity than that of the other TFs. In addition, low connectivity was found for LIM domain-containing protein WLIM2b (TRINITY_DN10999_c0_g1_i1_3_1), auxin response factor 9 (TRIN-ITY_DN6130_c0_g1_i1_3_1), protein TIFY 10 b (TRINITY_DN7781_c0_g1_i1_2_2), NAC transcription factor 56 (TRINITY_DN13925_c0_g1_i5_2_2), NAC transcription factor 56 (TRINITY_DN16303_c0_g1_i1_1_2), and pathogenesis-related genes transcriptional activator PTI6 (TRINITY_DN8922_c0_g1_i1_3_1), as indicated by lower correlations with other TFs. Thus, interconnections among TFs should be considered when investigating the regulation of SZ degeneration.

Flavone and Flavanol Biosynthesis Is Enriched in the Critical Degeneration Stage of the Stomium Zone
Stages S6 cm to S8 cm were identified as the critical ones in SZ degeneration. In the metabolomic analysis of the SZ samples (S6 cm, S7 cm, and S8 cm), 4962 metabolites were identified. To identify the differentially expressed metabolites (DEMs), four biological replicates of each development stage were analyzed. In the three comparison groups (S6 cm, S7 cm, and S8 cm), the number of metabolites that were upregulated was higher than the number of metabolites that were downregulated ( Figure S2). A volcano map shows the distribution of DEMs (Figure 7a). KEGG pathway-enrichment analysis of DEMs was performed (Figure 7b) using the metabolomic data set. Comparisons of DEMs among the three developmental stages indicated that within the KEGG categories (top 20), terms related to "flavone and flavanol biosynthesis" were enriched in all three comparison groups. In addition, low connectivity was found for LIM domain-containing protein WLIM2b (TRINITY_DN10999_c0_g1_i1_3_1), auxin response factor 9 (TRINITY_DN6130_c0_g1_i1_3 _1), protein TIFY 10 b (TRINITY_DN7781_c0_g1_i1_2_2), NAC transcription factor 56 (TRINITY_DN13925_c0_g1_i5_2_2), NAC transcription factor 56 (TRINITY_DN16303_c0_g1 _i1_1_2), and pathogenesis-related genes transcriptional activator PTI6 (TRINITY_DN8922_ c0_g1_i1_3_1), as indicated by lower correlations with other TFs. Thus, interconnections among TFs should be considered when investigating the regulation of SZ degeneration.

Flavone and Flavanol Biosynthesis Is Enriched in the Critical Degeneration Stage of the Stomium Zone
Stages S6 cm to S8 cm were identified as the critical ones in SZ degeneration. In the metabolomic analysis of the SZ samples (S6 cm, S7 cm, and S8 cm), 4962 metabolites were identified. To identify the differentially expressed metabolites (DEMs), four biological replicates of each development stage were analyzed. In the three comparison groups (S6 cm, S7 cm, and S8 cm), the number of metabolites that were upregulated was higher than the number of metabolites that were downregulated ( Figure S2). A volcano map shows the distribution of DEMs (Figure 7a

Reactive Oxygen Species Mediates Stomium Zone Degeneration.
To identify the main metabolites related to SZ degeneration, a metabolites database of the enrichment pathway from S6 cm-vs-S8 cm was studied (Table S5). In figure 8, there were four primary annotations shown. In "flavonoid biosynthesis pathways", one metabolite (kaempferol) was highly expressed in stage S8 cm. The "anthocyanin biosynthesis pathway" was also enriched. Chrysanthemin and petunidin 3-glucoside were upregulated from stage S6 cm to S8 cm. Regarding phytohormones, one metabolite of the "alphalinolenic acid metabolism pathway" was downregulated from stage S6 cm to S8 cm. By contrast, metabolites in the "starch and sucrose metabolism pathway" were upregulated from stage S6 cm to S8 cm. In addition, according to the FPKM values, the ROS-scavenging genes APX1, APX4, CAT2, GPX3, and GSTF1 were upregulated during the critical stage of SZ degeneration. By contrast, CAT1, SOD1, and SOD2 were upregulated and then downregulated in stages S6 cm to S8 cm. These results suggested ROS participated in SZ cell death.

Reactive Oxygen Species Mediates Stomium Zone Degeneration
To identify the main metabolites related to SZ degeneration, a metabolites database of the enrichment pathway from S6 cm-vs-S8 cm was studied (Table S5). In Figure 8, there were four primary annotations shown. In "flavonoid biosynthesis pathways", one metabolite (kaempferol) was highly expressed in stage S8 cm. The "anthocyanin biosynthesis pathway" was also enriched. Chrysanthemin and petunidin 3-glucoside were upregulated from stage S6 cm to S8 cm. Regarding phytohormones, one metabolite of the "alphalinolenic acid metabolism pathway" was downregulated from stage S6 cm to S8 cm. By contrast, metabolites in the "starch and sucrose metabolism pathway" were upregulated from stage S6 cm to S8 cm. In addition, according to the FPKM values, the ROS-scavenging genes APX1, APX4, CAT2, GPX3, and GSTF1 were upregulated during the critical stage of SZ degeneration. By contrast, CAT1, SOD1, and SOD2 were upregulated and then downregulated in stages S6 cm to S8 cm. These results suggested ROS participated in SZ cell death.

Stomium Zone Degeneration as a Process of Programmed Cell Death
TUNEL assays were used to detect the PCD of SZ cells at different stages of anther development (S5 cm to S9 cm) (Figure 9). In the control group, a TUNEL-positive signal was not observed in the anther. In addition, a TUNEL-positive signal was not observed in the SZ at stages S5 cm and S6 cm. However, at the later S6 cm stage, positive signals were observed in SZ cells, suggesting that PCD occurred. In the S7 cm stage, additional positive signals were observed in SZ cells. At the S8 cm stage, positive signals in stomium tissues were not as strong. In the last S9 cm stage, a positive signal was barely detected in anther stomium tissues. In addition, the FPKM values of eight PCD-related genes were analyzed. In the results, vacuolar processing enzyme (Vpe) was upregulated during the critical stage of SZ degeneration, higher than the S5 cm and S9 cm expression level. Chlorophyll b reductase Nyc1 and Nol had a similar expression trend with Vpes. Respiratory burst oxidase homologue (Rboh) proteins also showed a high expression level during the critical stage of SZ degeneration. These results suggested PCD participated in SZ degeneration.

Stomium Zone Degeneration as a Process of Programmed Cell Death
TUNEL assays were used to detect the PCD of SZ cells at different stages of anther development (S5 cm to S9 cm) (Figure 9). In the control group, a TUNEL-positive signal was not observed in the anther. In addition, a TUNEL-positive signal was not observed in the SZ at stages S5 cm and S6 cm. However, at the later S6 cm stage, positive signals were observed in SZ cells, suggesting that PCD occurred. In the S7 cm stage, additional positive signals were observed in SZ cells. At the S8 cm stage, positive signals in stomium tissues were not as strong. In the last S9 cm stage, a positive signal was barely detected in anther stomium tissues. In addition, the FPKM values of eight PCD-related genes were analyzed. In the results, vacuolar processing enzyme (Vpe) was upregulated during the critical stage of SZ degeneration, higher than the S5 cm and S9 cm expression level. Chlorophyll b reductase Nyc1 and Nol had a similar expression trend with Vpes. Respiratory burst oxidase homologue (Rboh) proteins also showed a high expression level during the critical stage of SZ degeneration. These results suggested PCD participated in SZ degeneration.

Jasmonates and Regulation of the Expression of Genes Associated with Stomium Zone Degeneration
Expression of genes correlated with PCD in lily was detected. In anthers with JA treatment, three genes involved in pectin degradation (TRINITY_DN15898_c0_g1_i1_2_2, TRINITY_DN22079_c0_g2_i4_2_2, and TRINITY_DN15576_c0_g1_i1_3_1) were downregulated, and one gene associated with water and sucrose transportation (TRIN-ITY_DN13124_c0_g1_i1_2_2) was upregulated. In addition, eight genes associated with cysteine protease (TRINITY_DN15815_c0_g1_i2_2_1, TRINITY_DN22555_c0_g1_i4_2_2,  TRINITY_DN11796_c0_g1_i2_2_1,  TRINITY_DN15662_c0_g1_i2_1_2,  TRIN-ITY_DN26951_c0_g1_i1_1_2,  TRINITY_DN5007_c0_g1_i3_2_2, TRIN-ITY_DN21546_c0_g1_i1_1_2, and TRINITY_DN22130_c0_g1_i6_2_2) were downregulated with JAs treatment (Figure 10). Thus, JAs could control anther dehiscence by regulating cell death and water transportation.

Jasmonates and Regulation of the Expression of Genes Associated with Stomium Zone Degeneration
Expression of genes correlated with PCD in lily was detected. In anthers with JA treatment, three genes involved in pectin degradation (TRINITY_DN15898_c0_g1_i1_2_2, TRIN-ITY_DN22079_c0_g2_i4_2_2, and TRINITY_DN15576_c0_g1_i1_3_1) were downregulated, and one gene associated with water and sucrose transportation (TRINITY_DN13124_c0_g1_ i1_2_2) was upregulated. In addition, eight genes associated with cysteine protease (TRIN-ITY_DN15815_c0_g1_i2_2_1, TRINITY_DN22555_c0_g1_i4_2_2, TRINITY_DN11796_c0_g1_ i2_2_1, TRINITY_DN15662_c0_g1_i2_1_2, TRINITY_DN26951_c0_g1_i1_1_2, TRINITY_DN 5007_c0_g1_i3_2_2, TRINITY_DN21546_c0_g1_i1_1_2, and TRINITY_DN22130_c0_g1_i6_2_ 2) were downregulated with JAs treatment (Figure 10). Thus, JAs could control anther dehiscence by regulating cell death and water transportation.

Discussion
Studies have revealed mechanisms of anther development in different plant species [9,20,21]; however, most research has focused on anther tapetal degeneration and endothecium cell lignification. Investigations of the molecular mechanisms of anther SZ degeneration are scarce. In this study, physiological, transcriptomic, and metabolic analyses were performed in five different stages of SZ development (S5 cm, S6 cm, S7 cm, S8 cm, and S9 cm). Based on the comparative transcriptome and metabolome data (DEGs and DEMs, respectively) obtained from the different SZ developmental stages, insights into the potential molecular mechanisms underlying anther dehiscence caused by SZ degeneration are provided.
Anther opening involves forces on anther walls that facilitate rupture of the stomium [3]. Anther stomium cells undergo a PCD-related process to facilitate pollen release [22,23]. In this study, histological and morphological analyses of lily anthers showed that anther dehiscence was closely correlated with SZ degeneration (Figure 1), consistent with previous studies [9,21].
Differentially expressed genes were screened on the basis of significant differences in expression in the anther SZ at different developmental stages ( Figure 2). Overall, the number of DEGs in the group related to S5 cm and S9 cm was higher than that in different groups from S6 cm to S8 cm. These results suggest the number of DEGs in S7 cm-vs-S6 cm and S8 cm-vs-S7 cm comparisons was very small. Developmental stage or environmental change can affect the number of DEGs [24]. In addition, from stage S6 cm to S8 cm in lily anther development, four microspores are released from the callose, tapetum cells gradually degenerate, and the pollen matures [21]. Consequently, stages S6 cm-S8 cm should be considered as the key period for SZ degeneration, ultimately causing anther dehiscence in lily (Figure 1).
RNA-seq is a revolutionary tool for transcriptomics, because the high-throughput sequencing technology provides an efficient and reliable platform for molecular biology research [25]. In this study, Illumina RNA-seq technology was used to study the molecular mechanisms of anther dehiscence in lily. GO term enrichment analysis was performed (Figure 3a). In the three categories of biological process, cellular component, and molecular function, the most enriched terms were "cellular process", "cell", and "binding", respectively. In addition, the four largest KEGG classification categories were "translation", "carbohydrate metabolism", "folding, sorting and degradation", and "energy metabolism" (Figure 3b). These results are not consistent with those in the analysis of chrysanthemum anther dehiscence. In chrysanthemums, most GO-annotated unigenes were in

Discussion
Studies have revealed mechanisms of anther development in different plant species [9,20,21]; however, most research has focused on anther tapetal degeneration and endothecium cell lignification. Investigations of the molecular mechanisms of anther SZ degeneration are scarce. In this study, physiological, transcriptomic, and metabolic analyses were performed in five different stages of SZ development (S5 cm, S6 cm, S7 cm, S8 cm, and S9 cm). Based on the comparative transcriptome and metabolome data (DEGs and DEMs, respectively) obtained from the different SZ developmental stages, insights into the potential molecular mechanisms underlying anther dehiscence caused by SZ degeneration are provided.
Anther opening involves forces on anther walls that facilitate rupture of the stomium [3]. Anther stomium cells undergo a PCD-related process to facilitate pollen release [22,23]. In this study, histological and morphological analyses of lily anthers showed that anther dehiscence was closely correlated with SZ degeneration (Figure 1), consistent with previous studies [9,21].
Differentially expressed genes were screened on the basis of significant differences in expression in the anther SZ at different developmental stages ( Figure 2). Overall, the number of DEGs in the group related to S5 cm and S9 cm was higher than that in different groups from S6 cm to S8 cm. These results suggest the number of DEGs in S7 cm-vs-S6 cm and S8 cm-vs-S7 cm comparisons was very small. Developmental stage or environmental change can affect the number of DEGs [24]. In addition, from stage S6 cm to S8 cm in lily anther development, four microspores are released from the callose, tapetum cells gradually degenerate, and the pollen matures [21]. Consequently, stages S6 cm-S8 cm should be considered as the key period for SZ degeneration, ultimately causing anther dehiscence in lily (Figure 1).
RNA-seq is a revolutionary tool for transcriptomics, because the high-throughput sequencing technology provides an efficient and reliable platform for molecular biology research [25]. In this study, Illumina RNA-seq technology was used to study the molecular mechanisms of anther dehiscence in lily. GO term enrichment analysis was performed (Figure 3a). In the three categories of biological process, cellular component, and molecular function, the most enriched terms were "cellular process", "cell", and "binding", respectively. In addition, the four largest KEGG classification categories were "translation", "carbohydrate metabolism", "folding, sorting and degradation", and "energy metabolism" (Figure 3b). These results are not consistent with those in the analysis of chrysanthemum anther dehiscence. In chrysanthemums, most GO-annotated unigenes were in the categories "metabolic process", "binding", and "catalytic activity" [26,27]. Results vary greatly because they might be affected depending on the species and also the tissue sampled. However, the results were also consistent with the biological processes involved in anther dehiscence, such as cell dehydration, hormone biosynthesis, and TF binding.
Metabolites are intermediate and final products that have important roles in regulating plant growth and development. The total number of plant metabolites is as high as 200,000, which indicates the diversity of plant natural substances [28]. Because of the high diversity, plant metabolites are ideal targets to study the regulation of biosynthesis [29]. Genetic mechanisms underlying changes in plant metabolites have been studied in recent years [30][31][32]. In this study, terms related to "flavone and flavanol biosynthesis" were enriched in all three comparison groups (S8 cm-vs-S6 cm, S7 cm-vs-S6 cm, and S8 cm-vs-S7 cm) (Figure 7). In addition, "starch and sucrose metabolism", "anthocyanin biosynthesis", and "alpha-linolenic acid metabolism" were also differentially enriched pathways. The metabolites involved in those metabolic pathways might be associated with anther dehiscence [3,22,33]. Morphological diagrams of developing stamens also showed pigment accumulation in degrading stomium tissues. Anther SZ cells had nuclear deformation and chromatin condensation and marginalization via PCD. Dong et al. [22] found that overexpression of the plantacyanin gene OPX causes indehiscent anthers and that endothecium degeneration in Arabidopsis plantacyanin overexpression anthers may be caused by plantacyanin-induced precocious PCD. Similarly, in Arabidopsis, the flavonoid transporter FFT is required for anther dehiscence [9]. In the transcriptome database in this study, OXP1 and FFT were upregulated during SZ degeneration (results not shown), and according to the metabolic database, flavanol and anthocyanin accumulation also occurred during SZ degeneration ( Figure 8). Furthermore, in this study, ROS-related gene expression indicated that SZ degeneration might be a consequence of ROS accumulation, especially H 2 O 2 in specific anther cell tissues. Flavonols could prevent ROS from reaching damaging levels and restore impaired fertility [34]. Anthocyanin-deficient mutants had relatively higher ROS levels than the wild type after paraquat treatment, indicating that anthocyanins function as antioxidants to protect against ROS in plants [35]. During the early stage, calcium and diphenyleneiodonium-sensitive reactive oxygen species (ROS) production are required to induce a secondary ROS burst and JA accumulation [36]. Although the cause of accumulation of flavanol and anthocyanin is unclear, ROS-related gene expression suggest that flavanol, anthocyanin, alpha-linolenic acid, and sucrose mediate stomium degeneration by regulating ROS homeostasis in plants.
Plant hormones have significant roles in growth and development and stress resistance. Studies on Arabidopsis and crops reveal phytohormones are involved in anther dehiscence, including JAs, auxins, gibberellins, and ethylene [37]. In this study, some auxin signal transduction-related genes (AUX1, IAA, GH3, ARF, SAUR) were expressed differently in different samples. Auxin regulates early stages of stamen development [38,39]. Auxin localization, biosynthesis, transport, and responses also affect late stages of stamen development, including anther dehiscence, pollen maturation, and preanthesis filament elongation in Arabidopsis, in addition to acting on the JAs levels in the initial stages of those processes [16]. In this study, most AUX1, IAA, GH3, ARF, and SAUR genes were highly expressed in stages S5 cm and S9 cm. This result demonstrated that those genes might regulate the expression of genes involved in SZ lysis during anther development.
JAs contribute to anther dehiscence [40]. Research on the role of JAs in anther dehiscence relies on the JAs signal transduction mutant coronatine insensitive 1 (coi1) [18]. COI1 is an F-box protein that recruits JAs and forms a complex with JAZ proteins (denoted for their ZIM and Jas motifs). The COI1-ligand-JAZ ternary complexes are polyubiquitinated and subsequently degraded by the 26S proteasome, which allows JA signal transduction [18]. In other studies, bHLH factors of JAs pathways antagonize transcription activators, such as MYC2 and the WD-repeat/bHLH/MYB complex, by binding to their target sequence [19,41,42]. In this study, MYC2 and JAZs were downregulated from stage S6 cm to S8 cm (Figure 4). MYC2 and JAZ proteins might repress the transcription of JA-responsive genes as direct regulating factors of SZ degeneration. Defects in anther dehiscence have been observed in mutants of the JA synthesis pathway [15]. JA synthase DAD1 control water transport into the vascular tissues from the anther endothecium, connective tissue, and anther locules, affecting anther dehiscence [43]. In this study, JA synthase was highly expressed in advance, which could provide JA content for the critical stage of anther dehiscence.
Anther dehiscence is also associated with ethylene and gibberellin pathways [44,45]. PhETR2 in petunia regulates SZ degeneration to control anther dehiscence [41], and GA (gibberellin)-dependent anther dehiscence can be cause by defects in secondary thickening or problems of hydration [46]. In this study, several GA and ethylene-related genes, including DELLA genes, were highly expressed from stage S6 cm to S8 cm. Gibberellins regulate expression of the JAs biosynthesis gene DAD1, and DELLA can act through JAs to control the expression of anther dehiscence-related genes [47]. Thus, the complexity of phytohormones provides a powerful buffer system that is crucial for plants to regulate the complex process of anther dehiscence.
At the cellular level, anther dehiscence likely involves cell wall-degrading enzymes that break down the pectin between cells [48]. Ricinosomes can harbor KDEL-tailed cysteine proteases that are involved in the final stages of cell death [49]. One of the cysteine proteases (SlCysEP) was detected early during tomato anther development in the stomium and epidermal cells surrounding the stomium [50]. In this research, compared with genes of pectin-related enzymes and those associated with anther wall dehydration, those of cysteine proteases had higher expression in the critical stages of anther dehiscence (S6 cm to S8 cm) ( Figure 5). Around the SZ of anthers, water transport genes may serve to increase osmotic potential and create a dehydrated environment to provide the final force for anther opening [51]. These processes jointly regulate SZ cell death and anther dehiscence.
Many TFs function as regulatory components of anther dehiscence. Several members of the MYB family, including MYB21, MYB24, MYB26, MYB57, and MYB108, regulate different transcriptional pathways of anther dehiscence. Overexpression those of MYB26 and MYB108 induces anther indehiscence by endothecium secondary thickening, which are also directly regulated by ARF8 and ARF17, respectively [10,12,15]. MYB24 is primarily involved in JA-signal transduction for filament elongation and SZ degradation [52,53]. In addition, MADS-box TF AGAMOUS is also involved in anther development, including anther dehiscence [54,55], and NAC TFs are proposed as regulatory components of endothecium lignification in Arabidopsis [56,57]. These results demonstrate that there are several key TFs that directly regulate anther dehiscence, although in different ways. However, the number of TFs identified that regulate anther dehiscence in a given plant remains small. Moreover, relations among TFs have not been fully explored. In this study, transcriptome sequencing data provided a foundation to distinguish regulatory networks and new regulators involved in anther dehiscence ( Figure 6). In the transcription data in this study, JAZs and the R2R3 MYB transcription factor MYB21 were also identified. JAZs interact with MYB21 and MYB24, two R2R3-MYB TFs essential for filament elongation, anther dehiscence, and pollen maturation, to mediate JA-regulated male fertility in Arabidopsis [19]. This finding is an example of the correlations between different TFs. The weighted expression network based on correlations will provide opportunities to more fully explore the mechanisms of transcription regulation in anther dehiscence.
The fuorochrome-based TUNEL assay is used to detect apoptotic DNA fragmentation [58] and therefore PCD. Lyu et al. [58] found low temperature affected the PCD of tapetum cells and also observed a positive TUNEL signal in stomium cells at low temperatures. In this study, stronger positive TUNEL signals were detected in SZ cells during stages S7 cm and S8 cm, whereas positive signals were not observed in the same tissue during stages S5 cm (Figure 9). The results indicated that SZ cell DNA began to break down at the critical stage of anther dehiscence. Apoptotic cells remained at the edge of the broken region in the later stage. Furthermore, in this study, 11 PCD-related genes had remarkably different expression with JA treatment. CEP1 and CEP2 are KDEL-tailed cysteine endopeptidases, and CEP1 reportedly regulates tapetal PCD [59]. RD21 is also a cysteine proteinase associated with fungal defense [60]. Under heat treatment, expression levels of anther aquaporin genes (PIPs and TIPs) are significantly correlated with anther dehiscence [61]. However, there are few reports on how the expression of these genes interacts with JA. This study suggested that JA treatment affected anther dehiscence by promoting the expression of genes related to water transport, pectin lysis, and PCD. In summary, on the basis of these analyses, a new understanding of the PCD mechanisms underlying anther SZ degeneration has been obtained.

Plant Materials and Sampling
The Oriental hybrid lily 'Siberia' was cultured in a greenhouse at the Lily Resource Preserving Center, Nanjing Agricultural University (Nanjing, China). Day/night temperatures were 26 • C/18 • C, and the relative humidity was 70%. Flowers were divided into nine developmental stages according to the length of flower buds (2 cm to 9 cm corresponded to stages S2 cm to S9 cm). Stomium zone samples were collected by freehand sectioning of anthers and then immediately frozen in liquid nitrogen and stored at −80 • C. All samples were taken from 20 flower buds and pooled to generate one biological replicate. Transcriptome sequencing used samples from stage S5 cm to S9 cm (each stage included three biological replicates), and metabolome sequencing used samples from stage S6 cm to S8 cm (each stage included four biological replicates).

RNA-Seq
Total RNA was extracted using an RNAprep Pure Plant Kit DP441 (Tiangen, Beijing, China), according to the manufacturer's instructions. RNA integrity was evaluated using an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). Samples with an RNA integrity number of 7 to 10 were subjected to subsequent analysis. The cDNA libraries were constructed using a HiScript II kit 117 (Vazyme, Nanjing, China). The prepared libraries were sequenced on an Illumina HiSeq TM 2500 platform ((Illumina, San Diego, CA, USA), and 125 bp/150 bp paired-end reads were generated. The database can be found in National Center of Biotechnology Information (NCBI) database under the accession number PRJNA767274.

RNA Data Processing and Analysis
Raw reads of fastqformat were first processed using trimmomatic [62]. To obtain clean reads, original reads and low-quality reads containing poly-n were deleted. After removing the adapters and low-quality sequences, clean reads were assembled into expression sequence tag clusters (contigs) by the Trinity (v2.0.6; http://trinityrnaseq.github.io/, accessed on 2 December 2019), and the transcripts were assembled from scratch. The longest transcript was selected as a single gene for subsequent analysis according to similarity and length. Differentially expressed genes (DEGs) were identified. After annotation, the read counts and FPKM [63] of each unigene were calculated using bowtie2 [64] and express [65], respectively. The threshold for significant differential expression was a p-value < 0.05 and fold change >2 or <0.5. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were performed on DEGs using R, based on a hypergeometric distribution [66,67].

Heat Map and Weighted Network Analyses
Tbtools [68] was used to prepare heat maps of the DEGs. Pearson's correlation coefficients were used to measure the co-expression relations in the Omicshare online program (https://www.omicshare.com/tools/, accessed on 2 June 2021). Only the associations with a p-value ≤ 0.05 were selected. In addition, the Omicshare online program (https://www.omicshare.com/tools/Home/Soft/cytoscape2, accessed on 22 June 2021) was used to construct the weighted network.

Reverse-Transcription Quantitative PCR
Total RNA isolated from the SZ of stage S6 cm to S9 cm treated with 0 or 50 µM MeJA (Bioduly, Nanjing, China) for 48 h. The cDNA was prepared using a HiScript II kit 117 (Vazyme, Nanjing, China). Quantitative RT-qPCR reactions were conducted using ChamQ Universal SYBR qPCR Master Mix (Vazyme, Nanjing, China) following the manufacturer's protocol. Primers for the genes were designed using Primer express 3.0 software (Table S6). Gene expression levels were normalized using the reference gene 18S [5], and values were calculated using the 2 −∆∆Ct method [70,71]. Expression values were obtained from three biological experiments. SPSS 17.0 statistical software was used to determine statistical significance.

Conclusions
In summary, SZ degeneration is essential for anther dehiscence. On the basis of the transcriptomic and metabolomic analysis of SZ degeneration, signaling of phytohormones, including JAs, affected the transcript abundances of the genes associated with pectinase, water and sucrose transport, and cysteine proteinase, which led to stomium cell death by regulating water transport, DNA breakage, and pectin degradation ( Figure 11). Thus, a morphological process of anther SZ degeneration in lily is proposed. This work also provides data to help reveal the regulatory mechanisms of SZ degeneration in plants to improve crop breeding.
added to the sample cells, which were then incubated for 1 h at 37 °C in the dark. Samples treated without reaction mixture were the controls. Images were taken using a laser scanning confocal microscope (LSM800, Zeiss, Oberkochen, Germany).

Reverse-Transcription Quantitative PCR
Total RNA isolated from the SZ of stage S6 cm to S9 cm treated with 0 or 50 μM MeJA (Bioduly, Nanjing, China) for 48 h. The cDNA was prepared using a HiScript II kit 117 (Vazyme, Nanjing, China). Quantitative RT-qPCR reactions were conducted using ChamQ Universal SYBR qPCR Master Mix (Vazyme, Nanjing, China) following the manufacturer's protocol. Primers for the genes were designed using Primer express 3.0 software (Table S6). Gene expression levels were normalized using the reference gene 18S [5], and values were calculated using the 2 −ΔΔCt method [70,71]. Expression values were obtained from three biological experiments. SPSS 17.0 statistical software was used to determine statistical significance.

Conclusions
In summary, SZ degeneration is essential for anther dehiscence. On the basis of the transcriptomic and metabolomic analysis of SZ degeneration, signaling of phytohormones, including JAs, affected the transcript abundances of the genes associated with pectinase, water and sucrose transport, and cysteine proteinase, which led to stomium cell death by regulating water transport, DNA breakage, and pectin degradation ( Figure 11). Thus, a morphological process of anther SZ degeneration in lily is proposed. This work also provides data to help reveal the regulatory mechanisms of SZ degeneration in plants to improve crop breeding.