Abstract

   Salix variegata, a typical dioecious plant with high reproductive and
   adaptive ability, has important ecological and ornamental value. To
   understand the potential mechanisms and metabolite dynamics of male and
   female flowers development, the first comparative analysis of the
   transcriptome and metabolome of S. variegata was applied. As a result,
   12,245 differentially expressed genes (DEGs) and 4,145 differently
   expressed metabolites (DEMs) were identified. Transcriptomic analysis
   showed that the male and female flowers development processes were
   related to phenylpropanoid and flavonoid biosynthesis. According to the
   metabolic profile, the main amino acids, flavonoids, phenylpropanoids,
   and their derivatives were accumulated during the development of male
   and female flowers of the S. variegata. Combined transcriptomic and
   metabolomic analyses indicated that the AUX/IAA, bHLH, MIKC, MYB, NAC,
   ERF and RLK transcription factors (TFs) and their associated key DEGs
   may mediate the metabolism of phenylpropanoids and flavonoids, which in
   turn regulate the development of male and female flowers in S.
   variegata. These results provide important insights to elucidate the
   development of male and female flowers of S. variegata at the molecular
   level. Our results will contribute to understanding the molecular and
   genetic mechanisms of male and female flower development in typical
   dioecious plants.

Supplementary Information

   The online version contains supplementary material available at
   10.1038/s41598-025-91317-0.

   Keywords: Salix variegata Franch., Flower development, Transcriptome,
   Metabolome, Coexpression network

   Subject terms: Developmental biology, Plant sciences, Systems biology

Introduction

   Plants are extremely important to the survival of human beings, as they
   can purify the air to produce oxygen, solve the problem of food sources
   for human beings, and cure human diseases, so exploring the biological
   mechanisms of plants is crucial to the survival and development of
   human beings. Previously, studies on the complex biology of plants
   mainly focused on their morphology, physiology, and biochemistry, and
   later began to analyze their intrinsic molecular mechanisms using a
   single-omics approach^[42]1,[43]2. Gradually, it has been found that it
   is difficult to fully characterize complex plant biological processes
   by a single genomics technique, and given this, a variety of combined
   genomics research techniques have been developed. By integrating
   information from different groups, the multi-omics approach can help
   people understand the regulation and causality between molecules by
   exploring the candidate key factors at a deeper level^[44]3.

   Plant growth and development is a highly complex process that involves
   the formation of various tissues and organs, color formation, and
   quality formation. For higher plants, flowering is an important part of
   their life history and is aimed at producing seeds and transmitting
   genetic information. Flower development is a prerequisite for
   successful sexual reproduction in plants, whether self-pollinated or
   cross-pollinated, and is regulated by the external environment^[45]4.
   The process of flower development in higher plants is the growth and
   development of plants to a certain stage, under the action of
   temperature, light, and other environmental factors, the hormones
   related to the regulation of flower bud differentiation continue to
   accumulate in the plant, and under the action of these hormones the
   meristematic tissue undergoes a series of morphological, structural and
   physiological changes, and its nutritive growth is aligned to no longer
   be transformed into leaf buds but rather to flower buds^[46]5. In
   recent years more and more genes and TFs related to flower development
   have been identified, such as AP2/ERF, MYB, bHLH, MADS-box, NAC, FAR1,
   TERF, Tify, and WRKY^[47]6–[48]13. These genes are responsible for the
   transformation of the plant’s stem-tip meristematic tissues to flower
   buds, then to the full bloom, final bloom, and finally to the formation
   of seeds.

   Salix variegata Franch. belonging to the Salicaceae family, is a
   typical dioecious plant, flood-resistant, drought-resistant, and
   resistant to heavy metal stress, it is an excellent soil-fixing and
   dike protection plant, but also has great ecological restoration
   potential^[49]14,[50]15. S. variegata is widely distributed in eastern
   Tibet, northern Yunnan, Guizhou, Sichuan, Chongqing, and southeastern
   Gansu in China, and often grows at riverbanks in sandy soils, pebble
   fields, and between rock crevices^[51]16. River levels subside, and
   plants in the riparian zone experience alternating annual stresses of
   flooding and drought, which also affects the growth and development of
   plants in the riparian zone^[52]17–[53]19. As a dominant species in the
   riparian zone, S. variegata is well adapted to changes in water levels
   in the riparian zone. Previous researchers have conducted many studies
   around its ecological adaptations, such as seed germination,
   physiology, and biochemistry under flooding stress, reproduction,
   morpho-physiological adaptations under heavy metal cadmium stress, soil
   microorganisms under double flooding-cadmium stress, physiological and
   biochemical characteristics, and transcriptional and metabolic traits
   of S. variegata under flooding stress^[54]15,[55]20–[56]27. Although
   many studies have been conducted on S. variegata, few have explored the
   underlying molecular mechanisms of its ecological adaptations. It has
   been found that S. variegata can rapidly recover its growth after
   flooding, and it can flower and fruit normally, its seed activity is
   also strong^[57]28. S. variegata, as a typical dioecious plant,
   generally flowers from June to December. However, there is little
   available information on the molecular mechanisms and metabolic
   dynamics of its male and female flower development. In this study,
   transcriptomic and metabolomic techniques were used to characterize the
   transcriptome and metabolome of S. variegata during male and female
   flowers development. This study had the following two main objectives:
   (1) to investigate the regulatory pathways and genes involved in the
   development of male and female flowers in S. variegata, and (2) to
   determine the dynamics of metabolites during the development of male
   and female flowers in S. variegata. This study intends to provide a
   basis for further research on sexual reproduction and breeding of S.
   variegata and even dioecious plants by elaborating, for the first time,
   the molecular genetic mechanism of male and female flowers development
   in S. variegata.

Results

Transcriptome variations at different developmental stages of S. variegata
male and female flowers

   Sequencing of 18 cDNA libraries constructed from three different floral
   developmental stages of S. variegata male and female flowers yielded a
   total of 114.85 Gb of clean data, with more than 5.92 Gb of clean data
   and Q30 ≥ 94.65% for each sample (Supplementary Table [58]S1).
   Principal component analysis (PCA) results showed that the different
   groupings of samples could be differentiated from each other and that
   the male and female flowers showed similar expression patterns at the
   same floral stage of S. variegata, with PC1 (23.11%) and PC2 (11.64%)
   separating the two (Supplementary Fig. [59]S1A).

   A total of 12,245 DEGs were screened by Fold Change (FC) ≥ 1.5, False
   Discovery Rate (FDR) < 0.01 (Supplementary Table [60]S2). Among them,
   there were 2,372, 1,232, and 1,103 DEGs in F1 vs. F2, F1 vs. F3, and F2
   vs. F3, respectively. 6,981, 6,684, and 1,003 DEGs were found in M1 vs.
   M2, M1 vs. M3, and M2 vs. M3, respectively. 517, 4,982, and 2,587 DEGs
   were found in F1 vs. M1, F2 vs. M2, F3 vs. M3, respectively. Among the
   three comparison groups of female flowers, the F1 vs. F2 screened the
   most DEGs, with 1,323 up-regulated genes and 1,049 down-regulated
   genes, indicating that more genes are involved in the development of
   female flowers of S. variegata from bud to anthesis, while the F2 vs.
   F3 screened the fewest DEGs, with 1,103 DEGs. Among the three
   comparison groups of male flowers, the M1 vs. M2 screened the most
   DEGs, with 3,812 up-regulated genes and 3,169 down-regulated genes,
   suggesting that more genes are involved in the development of S.
   variegata male flowers from flower bud to anthesis, while the M2 vs. M3
   had the fewest DEGs, with 1,003 DEGs. In the comparison group of male
   and female flowers at the same flowering stage, DEGs were most abundant
   in F2 vs. M2 (4,982), containing 2,464 up-regulated and 2,518
   down-regulated genes, and least abundant in F1 vs. M1 (517), with 162
   up-regulated genes and 355 down-regulated genes. Male and female flower
   DEGs were categorized separately using Venn diagrams. There were 20
   DEGs shared among the three periods of female flowers, with 1,281, 374,
   and 650 specifically expressed DEGs in F1 vs. F2, F1 vs. F3, and F2 vs.
   F3, respectively. 383 DEGs shared among the three periods of male
   flowers, with 2,208, 1,680, and 81 specifically expressed DEGs in M1
   vs. M2, M1 vs. M3, and M2 vs. M3, respectively. 54 DEGs shared among
   the three comparison groups of male and female flowers at the same
   flowering stage, and 292, 3,646, and 1,300 specifically expressed DEGs
   in F1 vs. M1, F2 vs. M2, and F3 vs. M3, respectively (Fig. [61]1A–C).

Fig. 1.

   [62]Fig. 1
   [63]Open in a new tab

   Characterization of DEGs for three periods of S. variegata male and
   female flowers (A) Venn diagram of DEGs for three comparison groups of
   S. variegata female flowers, (F1) flower bud stage, (F2) full bloom
   stage, and (F3) final bloom stage. (B) Venn diagram for three
   comparison groups of S. variegata male flowers, (M1) flower bud stage,
   (M2) full bloom stage, (M3) final bloom stage.(C) Venn diagram for the
   three comparison groups of male and female flowers at the same
   flowering stage. (D) KEGG pathway enrichment analysis of DEGs in F1 vs.
   F2 groups. (E) KEGG pathway enrichment analysis of DEGs in the F1 vs.
   F3 group. (F) KEGG pathway enrichment analysis of DEGs in the F2 vs. F3
   group. (G) KEGG pathway enrichment analysis of DEGs in the M1 vs. M2
   group. (H) KEGG pathway enrichment analysis of DEGs in the M1 vs. M3
   group. (I) KEGG pathway enrichment analysis of DEGs in the M2 vs. M3
   group. (J) KEGG pathway enrichment analysis of DEGs in the F1 vs. M1
   group. (K) KEGG pathway enrichment analysis of DEGs in the F2 vs. M2
   group. (L) KEGG pathway enrichment analysis of DEGs in the F3 vs. M3
   group. The top 20 metabolic pathways significantly enriched for DEGs in
   different subgroups are shown. Triangles represent up-regulated genes,
   inverted triangles represent down-regulated genes, and circles
   represent genes that are both up-regulated and down-regulated, and the
   size of these three images represents the number of DEGs in the
   pathway. Colors represent significance.

   The clustering heat map was able to provide an overall picture of the
   expression of DEGs in different developmental stages of S. variegata
   male and female flowers. As can be seen from the figure, both male and
   female flowers of S. variegata showed significant differences in the
   expression of different genes and different developmental stages, and
   the biological replicates were highly similar among the groups
   (Supplementary Fig. [64]S1B). The expression of genes was categorized
   into three clusters in total, with the first cluster of genes being
   more highly expressed at the bloom stage of male flowers and least
   expressed at the bud stage of male flowers. The genes in the second
   cluster were most highly expressed at the bloom stage of female flowers
   and least highly expressed at the terminal flower stage of male
   flowers. The third cluster of genes had the highest expression at the
   bud stage of male flowers and the lowest expression at the final flower
   stage of male flowers. The expression of cluster I and cluster II genes
   tended to increase with flower development in F, and the expression of
   cluster III genes tended to decrease with female flower development.
   The expression of cluster I genes first increased and then decreased
   with the development of male flowers, and the expression of cluster III
   genes gradually decreased with the development of male flowers.

Functional analysis of DEGs

   GO and KEGG enrichment analyses were performed on the DEGs produced
   during different flower stages of S. variegata in both male and female
   flowers, and the results of GO analyses showed that the DEGs produced
   during the three developmental stages of S. variegata in both male and
   female flowers were mainly enriched in the areas of biological
   processes, cellular components, and molecular functions. Among them,
   female flowers were significantly enriched in detoxification, rhythmic
   process, transporter activity, and transcription regulator activity;
   male flowers were mainly enriched in detoxification, transporter
   activity, and transcription regulator activity (Supplementary Fig.
   [65]S1C-E). 126, 114 and 113 pathways were enriched in the F1 vs. F2,
   F1 vs. F3 and F2 vs. F3 comparison groups of female flowers,
   respectively. F1 vs. F2 was significantly enriched in pathways such as
   Plant hormone signal transduction (ko04075), Photosynthesis (ko00195),
   ABC transporters (ko02010), Carbon fixation in photosynthetic
   organisms (ko00710), Fatty acid elongation (ko00062), and Pentose
   phosphate pathway (ko00030). F1 vs. F3 was significantly enriched in
   pathways such as Plant hormone signal transduction, Pentose phosphate
   pathway, Carbon metabolism (ko01200), Photosynthesis, Carbon fixation
   in photosynthetic organisms, and ABC transporters. F2 vs. F3 was
   significantly enriched in pathways such as Plant − pathogen
   interaction (ko04626), Fatty acid elongation, Starch and sucrose
   metabolism (ko00500), Glycerolipid metabolism (ko00561), and Steroid
   biosynthesis (ko00100). 133, 135 and 108 pathways were enriched in the
   M1 vs. M2, M1 vs. M3, and M2 vs. M3 comparison groups of male flowers,
   respectively. M1 vs. M2 was significantly enriched in pathways such as
   Plant hormone signal transduction, Flavonoid biosynthesis (ko00941),
   other glycan degradation (ko00511), ABC transporters, Valine, leucine,
   and isoleucine degradation (ko00280), and Fatty acid elongation. M1 vs.
   M3 was significantly enriched in pathways such as Plant-pathogen
   interaction, Plant hormone signal transduction, MAPK signaling
   pathway-plant (ko04016), Glutathione metabolism (ko00480), Flavonoid
   biosynthesis, and beta − Alanine metabolism (ko00410). M2 vs. M3 was
   significantly enriched in pathways such as Plant hormone signal
   transduction, Cutin, suberine and wax biosynthesis (ko00073),
   Glycerolipid metabolism, ABC transporters, and Pentose and glucuronate
   interconversions (ko00040). 75, 133 and 126 pathways were enriched in
   the F1 vs. M1, F2 vs. M2, F3 vs. M3, respectively. F1 vs. M1 was
   significantly enriched in pathways such as Plant-pathogen interaction,
   Plant hormone signal transduction, Starch and sucrose metabolism, MAPK
   signaling pathway - plant, alpha-Linolenic acid metabolism (ko00592),
   and Linoleic acid metabolism (ko00591). F2 vs. M2 was significantly
   enriched in pathways such as Fatty acid elongation, Glycerolipid
   metabolism, Starch and sucrose metabolism, Carbon metabolism, and
   Photosynthesis - antenna proteins. F3 vs. M3 was significantly enriched
   in pathways such as Phenylpropanoid biosynthesis, Flavonoid
   biosynthesis, Benzoxazinoid biosynthesis (ko00402), Plant hormone
   signal transduction, and Photosynthesis - antenna proteins (ko00196)
   (Fig. [66]1D-L). Although there were differences in the pathways of DEG
   enrichment between male and female flowers at different developmental
   stages, it is still evident that DEGs are mainly involved in the
   synthesis of flavonoids and amino acids, among others, as male and
   female flowers develop.

Cluster analysis of DEGs

   We obtained the expression trends of 12,245 DEGs during the development
   of male and female flowers of S. variegata by K-means clustering
   analysis. All DEGs co-clustered into 18 clusters in female flowers and
   17 clusters in male flowers (Fig. [67]2, Supplementary Fig. [68]S2).
   Similar expression trends were clustered into the same clusters in
   testing the potential role of DEGs, high expression in F3 and F2 (e.g.,
   Cluster 2, cluster 3, cluster 5, cluster 6, and cluster 13) had
   Phenylpropanoid biosynthesis (ko00940), Nitrogen metabolism (ko00910),
   Galactose metabolism (ko00052), Flavonoid biosynthesis, Steroid
   biosynthesis, and Photosynthesis. Highly expressed in F1 (e.g.,
   cluster7, cluster9, cluster10, cluster14, cluster15, and cluster17)
   were Plant hormone signal transduction, Glycolysis / Gluconeogenesis
   (ko00010), Flavonoid biosynthesis, Glycerolipid metabolism, Betalain
   biosynthesis (ko00965), Plant-pathogen interaction, and
   Phosphatidylinositol signaling system (ko04070). Highly expressed in M3
   (e.g., cluster1, cluster4, cluster8, cluster9, and cluster15) were
   Phenylpropanoid biosynthesis, Plant hormone signal transduction,
   Glycerolipid metabolism, and Plant-pathogen interaction. Highly
   expressed in M2 (cluster3, cluster5, cluster7, cluster11, cluster12,
   and cluster17) were Pentose and glucuronate interconversions, Flavonoid
   biosynthesis, Phenylpropanoid biosynthesis, Nitrogen metabolism
   (ko00910), Glucosinolate biosynthesis (ko00966), and Carotenoid
   biosynthesis (ko00906). Highly expressed in M1 (e.g., cluster2,
   cluster6, cluster10, cluster14, and cluster15) were DNA replication
   (ko03030), Plant hormone signal transduction, Starch and sucrose
   metabolism, Linoleic acid metabolism, and Glycerolipid metabolism.
   K-means clustering analysis of DEGs from different periods of male and
   female flowers yielded results similar to the functional analysis of
   DEGs, which mainly regulated the biosynthesis of flavonoids,
   phenylpropanoids, and a variety of amino acids during the development
   of male and female flowers, suggesting that the biosynthesis of these
   substances is a key pathway in the development of male and female
   flowers of S. variegata.

Fig. 2.

   [69]Fig. 2
   [70]Open in a new tab

   K-Means clustering analysis of DEGs for three developmental periods of
   S. variegata female flowers.

Transcription factor analysis

   To understand the underlying molecular mechanisms of male and female
   flower development in S. variegata, highly expressed TFs were further
   identified and analyzed. Interestingly, we found a large number of TFs
   in different flowering periods of S. variegata male and female flowers.
   Among them, RLK, ERF, MYB, bHLH, C2H2, MIKC_MADS, and NAC TFs were
   expressed in different flowering periods of male and female flowers
   (Supplementary Fig. S3). Overall, TFs were more highly expressed at all
   stages of male than female flowers and more DEGs were down-regulated
   than up-regulated at different times in both male and female flowers.
   Among the predicted gene families, RLK, C2H2, and WRKY were
   significantly up-regulated in F2 and F3 of female flowers, and ERF,
   MYB, bZIP, and TALE were significantly down-regulated in F2 and F3 of
   female flowers. RLK, NAC, and WRKY were significantly up-regulated in
   M2 and M3 of male flowers, and MYB, bHLH, C2H2, and MIKC_MADS were
   significantly down-regulated in M2 and M3 of male flower. ERF, MYB, and
   bHLH were significantly down-regulated and RLK, MIKC were significantly
   up-regulated in M2 compared to F2. Compared with F3, RLK, bHLH, ERF,
   and C2H2 were significantly down-regulated and NAC, WRKY, MIKC were
   significantly up-regulated in M3 (Supplementary Table S3).

Statistical analysis of metabolites

   A total of 5,855 metabolites were detected and quantified by LC-MS/MS,
   containing compounds such as organic acids and their derivatives,
   organic heterocyclic compounds, lipids, benzoic acids,
   phenylpropanoids, and polyketides. Among them, more substances were
   fatty acids, terpene glycans, diterpenes, linolenic acid, and its
   derivatives, etc. The results of PCA analyses showed that the
   metabolites of different periods of male and female flowers were
   differentiated. The metabolites of F3 and M3 could be clearly
   distinguished from the other periods of male and female flowers,
   suggesting that the end of the female flower and the end of the male
   flower were the most different (Supplementary Fig. [71]S4A). The
   accumulation characteristics of metabolites in S. variegata female and
   male flowers at different developmental stages could be analyzed by
   clustering heat maps, and the results showed that there were
   significant differences in metabolites between male and female flowers
   as well as at different developmental stages of the same type of flower
   (Supplementary Fig. S4B). The metabolites were categorized into four
   clusters in total, with the metabolites in cluster 1 being the most
   abundant in F2 and the least abundant in M3. Metabolites in cluster 2
   were most abundant in F3 and least abundant in M3. Metabolites in
   cluster 3 were most abundant in M3 and least abundant in M1.
   Metabolites in cluster 4 were most abundant in M1 and least abundant in
   F2. From both principal component analysis and cluster heat map
   analysis, it can be concluded that the metabolites of S. variegata
   female and male flowers changed significantly in different
   developmental stages.

DEMs function analysis

   A total of 4,145 DEMs were screened based on the orthogonal projections
   to latent structures- discriminant analysis (OPLS-DA) model combined
   with variable importance in projection (VIP) ≥ 1, FC ≥ 1, and
   P-value ≤ 0.01 (Supplementary Table S4). Among the female flowers, the
   F2 vs. F3 comparison group had the most DEMs (1,509), of which 842
   metabolites were up-regulated and 667 metabolites were down-regulated.
   Among male flowers, the M1 vs. M3 comparison group had the highest
   number of DEMs (2,194), of which 1,324 metabolites were up-regulated
   and 870 metabolites were down-regulated. The most DEMs were accumulated
   in F3 vs. M3 and the least in F1 vs. M3 during the same flowering
   period in female and male flowers. There was a significant increase in
   the content of 1,014 DEMs and a substantial decrease in the content of
   690 DEMs in M3 compared to F3. The number and accumulation of DEMs
   increased gradually with the opening of both male and female flowers
   (Supplementary Fig. S4C-K).

   KEGG enrichment analysis between different comparison groups of female
   and male flowers showed that 87, 81, and 102 metabolic pathways were
   annotated in the three comparison groups of female flowers,
   respectively. 11, 5, and 8 of these pathways were significantly
   enriched, respectively. The three comparison groups of male flowers
   were annotated to 98, 99, and 102 metabolic pathways, respectively. 4,
   9, and 7 of these pathways were significantly enriched, respectively.
   The three comparison groups of female and male flowers of the same
   flowering stage were annotated to 71, 95, and 98 metabolic pathways.
   Nicotinate and nicotinamide metabolism (ko00760), Arginine biosynthesis
   (ko00220), Alanine, aspartate, and glutamate metabolism (ko00250), and
   Phenylalanine, tyrosine and tryptophan biosynthesis (ko00400) metabolic
   pathways were significantly enriched in F1 vs. F2. Phenylpropanoid
   biosynthesis, Betalain biosynthesis (ko00965), Glucosinolate
   biosynthesis (ko00966), Pyrimidine metabolism (ko00240), and Vitamin B6
   metabolism (ko00750) metabolic pathways were significantly enriched in
   F1 vs. F3. Phenylpropanoid biosynthesis, Carbon fixation in
   photosynthetic organisms, Alanine, aspartate and glutamate metabolism,
   Pentose and glucuronate interconversions, and Pyruvate metabolism
   (ko00620) metabolic pathways were significantly enriched in F2 vs. F3.
   Pyruvate metabolism, Aflatoxin biosynthesis (ko00254), Flavonoid
   biosynthesis, and Phosphatidylinositol signaling system (ko04070)
   metabolic pathways were significantly enriched in M1 vs. M2.
   Phenylpropanoid biosynthesis, Monoterpenoid biosynthesis, Neomycin,
   kanamycin, and gentamicin biosynthesis, Phenylalanine metabolism, and
   Tyrosine metabolism (ko00350) metabolic pathways were significantly
   enriched in M1 vs. M3. Vitamin B6 metabolism, Pentose and glucuronate
   interconversions, Nicotinate and nicotinamide metabolism, Glyoxylate
   and dicarboxylate metabolism(ko00630), and Nitrogen metabolism
   metabolic pathways were significantly enriched in M2 vs. M3.
   Phenylalanine, tyrosine and tryptophan biosynthesis, Pentose and
   glucuronate interconversions, Glyoxylate and dicarboxylate metabolism,
   and Fatty acid elongation (ko00062) metabolic pathways were
   significantly enriched in F1 vs. M1. Flavone and flavonol biosynthesis,
   Isoflavonoid biosynthesis (ko00943), Phenylpropanoid biosynthesis, and
   Flavonoid biosynthesis metabolic pathways were significantly enriched
   in F2 vs. M2. Flavone and flavonol biosynthesis, Monoterpenoid
   biosynthesis (ko00902), and Alanine, aspartate and glutamate metabolism
   (ko00250) metabolic pathways were significantly enriched in F3 vs. M3.
   Among them, Phenylpropanoid biosynthesis and Flavonoid biosynthesis
   were more enriched at the full and final bloom stages (Fig. [72]3).

Fig. 3.

   [73]Fig. 3
   [74]Open in a new tab

   KEGG metabolic pathway enrichment analysis of DEMs between different
   groupings of S. variegata male and female flowers (A) KEGG metabolic
   pathway enrichment of DEMs in F1 vs. F2. (B) KEGG metabolic pathway
   enrichment of DEMs in F1 vs. F3. (C) KEGG metabolic pathway enrichment
   of DEMs in the F2 vs. F3. (D) KEGG metabolic pathway enrichment of DEMs
   in the M1 vs. M2. (E) KEGG metabolic pathway enrichment of DEMs in the
   M1 vs. M3. (F) KEGG metabolic pathway enrichment of DEMs in the M2 vs.
   M3. (G) KEGG metabolic pathway enrichment of DEMs in the F1 vs. M1. (H)
   KEGG metabolic pathway enrichment of DEMs in the F2 vs. M2. (I) KEGG
   metabolic pathway enrichment of DEMs in the F3 vs. M3. The top 20
   metabolic pathways significantly enriched in DEGs in different
   subgroups are shown. Triangles represent up-regulated genes, inverted
   triangles represent down-regulated genes, and circles represent both
   up-and-down-regulated genes, and the size of these three images
   represents the number of DEGs in the pathway. Colors represent
   significance.

Transcriptome and metabolome in combination analysis

   Common metabolic pathways were obtained by comparing the metabolic
   pathways involved in DEGs vs. DEMs. The results of the Venn diagram
   showed that there were 79, 70, and 84 common metabolic pathways between
   DEGs and DEMs in F1 vs. F2, F1 vs. F3, and F2 vs. F3, respectively.
   There were 92, 95, and 80 common metabolic pathways between DEGs and
   DEMs in M1 vs. M2, M1 vs. M3, and M2 vs. M3, respectively. There were
   35, 90, and 89 common metabolic pathways between DEGs and DEMs in F1
   vs. M1, F2 vs. M2, F3 vs. M3 (Fig. [75]4A–I). We selected the 30
   pathways with the most significant enrichment of common pathways in
   each subgroup for enrichment analysis. The results showed that
   Phenylpropanoid biosynthesis, Pentose and glucuronate interconversions,
   and Flavonoid biosynthesis pathways were the most significant metabolic
   pathways during the development of male and female flowers of
   S.variegata (Fig. [76]4J–R).

Fig. 4.

   [77]Fig. 4
   [78]Open in a new tab

   Analysis of common pathways of transcriptome and metabolome at
   different developmental stages of S. variegata male and female flowers
   (A) Venn diagram of the pathways of DEGs and DEMs in the F1 vs. F2. (B)
   Venn diagram of the pathways of DEGs and DEMs in the F1 vs. F3. (C)
   Venn diagram of the pathways of DEGs and DEMs in the F2 vs. F3. (D)
   Venn diagram of the pathways of DEGs and DEMs in the M1 vs. M2. (E)
   Venn diagram of the pathways of DEGs and DEMs in the M1 vs. M3. (F)
   Venn diagram of the pathways of DEGs and DEMs in the M2 vs. M3. (G)
   Venn diagram of the pathways of DEGs and DEMs in the F1 vs. M1. (H)
   Venn diagram of the pathways of DEGs and DEMs in the F2 vs. M2. (I)
   Venn diagram of the pathways of DEGs and DEMs in the F3 vs. M3. (J)
   Co-enrichment pathway analysis of DEGs and DEMs in F1 vs. F2. (K)
   Co-enrichment pathway analysis of DEGs and DEMs in F1 vs. F3. (L)
   Co-enrichment pathway analysis of DEGs and DEMs in F2 vs. F3. (M)
   Co-enrichment pathway analysis of DEGs and DEMs in the M1 vs. M2. (N)
   Co-enrichment pathway analysis of DEGs and DEMs in the M1 vs. M3. (O)
   Co-enrichment pathway analysis of DEGs and DEMs in the M2 vs. M3. (P)
   Co-enrichment pathway analysis of DEGs and DEMs in the F1 vs. M1. (Q)
   Co-enrichment pathway analysis of DEGs and DEMs in the F2 vs. M2. (R)
   Co-enrichment pathway analysis of DEGs and DEMs in the F3 vs. M3.

   To elucidate the relevance of DEGs and DEMs related to
   phenylpropanoids, flavonoids, etc., we constructed the integration
   pathways for the three developmental stages of male and female flowers
   of S. variegata (Fig. [79]5A). In the biosynthesis pathway of
   phenylpropanoid and flavonoids, the expression of 4-coumarate–CoA
   ligase (4CL: Sva01G002770, Sva01G004170, Sva03G006690, Sva03G014220,
   Sva06G013710, Sva10G004990, Sva12G008250, Sva12G008260, Sva15G009870,
   Sva17G002660, Sva17G012550), Cinnamyl-alcohol dehydrogenase (CAD:
   Sva09G007840, Sva16G017940), Ferulate-5-hydroxylase (F5H: Sva05G008980,
   Sva07G001420), and Catechol-O-Methyltransferase (COMT: Sva12G000240,
   Sva12G000580, Sva14G007990, Sva14G008010, Sva16G028990, Sva16G029010)
   was F3 > F2 > F1, M3 > M2 > M1, and with the development of male and
   female flowers, the high expression of these genes promoted the
   accumulation of Ferulic acid, Sinapic acid, 5-Hydroxyconiferaldehyde,
   Sinapoyl-CoA, and Sinapyl alcohol accumulation, which ultimately
   promoted the synthesis of large amounts of Syringin. The expression of
   Anthocyanidin synthase (ANS: Sva01G009400, Sva03G008470, Sva06G008290,
   Sva06G008300, Sva16G010160), Anthocyanidin 3-O-glucosyltransferase
   (BZ1: Sva09G011020, Sva09G011040, Sva09G011050, Sva13G011910), Flavonol
   synthase (FLS: Sva03G001460, Sva10G014730, Sva19G000640), and
   Anthocyanidin 3-O-glucoside 2’’’-O-xylosyltransferase (UGT79B1:
   Sva11G006400) was highest at the full bloom stage of both male and
   female flowers, and Bifunctional dihydroflavonol 4-reductase (DFR)
   subsequently converted the two main substrates (dihydrokaempferol,
   dihydroquercetin) into Leucoyanidin and Leucopelargonidin. With the
   development of both male and female flowers of S. variegata, the
   expression of ANS, DFR, and BZ1 genes were gradually down-regulated for
   expression at subsequent developmental stages, whereas the increased
   expression of FLS promoted the synthesis of Kaempferol, which
   ultimately synthesized large amounts of
   Cyanidin-5-glucoside-3-sambubioside,
   Pelagoridin-5-glucoside-3-sambubioside, and Quercetin 3-O-rhamnoside
   7-O-glucoside. Some of the metabolites were consistent with their
   upstream gene expression, such as Sinapic acid, Syringin, and
   Cyanidin-5-glucoside-3-sambubioside. It can be seen that during the
   development of male and female flowers of S. variegata, the content of
   phenylpropanoids and flavonoids gradually increased with the
   development of flowers.

Fig. 5.

   [80]Fig. 5
   [81]Open in a new tab

   Joint analysis of DEGs and DEMs involved in different developmental
   stages of male and female flowers of S. variegata (A) Integration
   pathways of DEGs and DEMs involved in phenylpropanoid and flavonoid
   biosynthesis. Pathways were remapped according to ko00940, ko00941,
   ko00942, ko00944, ko00945 in the KEGG database
   ([82]www.kegg.jp/kegg/kegg1.html. The Kanehisa laboratory have happily
   provided permission). Red and green colors indicate metabolites with
   relatively high levels and metabolites with relatively low levels,
   respectively. Red and blue colors indicate genes with up-regulated
   expression and genes with down-regulated expression, respectively. From
   left to right, the six squares represent the three extremes of S.
   variegata male and female flower. COMT (Sva12G000240, Sva12G000580,
   Sva14G007990, Sva14G008010, Sva16G028990, Sva16G029010). CSE,
   caffeoylshikimate esterase (Sva01G014390, Sva02G016820, Sva08G003850,
   Sva10G009810, Sva10G016780, Sva14G009920, Sva13G007050). 4CL
   (Sva01G002770, Sva01G004170, Sva03G006690, Sva03G014220, Sva06G013710,
   Sva10G004990, Sva12G008250, Sva12G008260, Sva15G009870, Sva17G002660,
   Sva17G012550); F5H (Sva05G008980, Sva07G001420). CCR, cinnamoyl-CoA
   reductase (Sva01G003320, Sva01G003330, Sva01G003340, Sva01G011550,
   Sva01G018420, Sva02G000420, Sva02G000430, Sva03G006210, Sva03G013750,
   Sva06G014290, Sva09G004650, Sva09G006230, Sva13G008720, Sva16G014550,
   Sva17G012280). CAD (Sva09G007840, Sva16G017940). UGT72E,
   coniferyl-alcohol glucosyltransferase (Sva07G002590). E2.1.1.104,
   caffeoyl-CoA O-methyltransferase (Sva09G008160, Sva10G007240,
   Sva16G018260). F3H, naringenin 3-dioxygenase (Sva05G008660,
   Sva16G030040). CYP75B1, flavonoid 3’-monooxygenase (Sva13G007050). DFR,
   (Sva02G002830, Sva03G010140, Sva08G011450). ANS (Sva01G009400,
   Sva03G008470, Sva06G008290, Sva06G008300, Sva16G010160). FLS
   (Sva03G001460, Sva10G014730, Sva19G000640). BZ1 (Sva09G011020,
   Sva09G011040, Sva09G011050, Sva13G011910). UGT79B1 (Sva11G006400). (B)
   Correlation heat map of gene expression of each module and important
   DEMs in S. variegata. (C) Network of relationships between TFs, DEGs,
   and metabolites involved in the phenylpropanoid and flavonoid pathways.
   The green, pink, and rose node sub-tables represent key genes, key TFs,
   and key metabolites related to S. variegata flower development.

   We further utilized 12,245 genes and 17 important DEMs related to
   synthesizing phenylpropanoids and flavonoids as source data for a
   weighted gene co-expression network module for analysis (WGCNA). Three
   modules- blue, turquoise, and grey- were identified in the cluster
   tree. The results of the relationship between the gene modules and the
   significant DEMs indicated that the blue and turquoise modules had a
   higher correlation with them. Genes in the blue module were
   significantly positively associated with Syringin and Sinapic_acid, and
   significantly negatively associated with
   Pelagoridin_5_glucoside_3_sambubioside, Dihydromyricetin, Sinapoyl_CoA,
   and Leucopelargonidin were significantly negatively correlated. Genes
   in the turquoise module were associated with
   Pelagoridin_5_glucoside_3_sambubioside, Dihydromyricetin, Sinapoyl_CoA,
   Dihydroquercetin, and Leucopelargonidin significantly positively and
   significantly negatively associated with Syringin and Sinapic_acid
   (Fig. [83]5B). 40 hub genes were screened based on correlation
   coefficients between candidate genes in the blue and turquoise modules,
   where correlation coefficients r > 0.8, and P-values < 0.05. Hub genes,
   high-expression TFs, and 17 significant DEMs were utilized to construct
   networks of interactions between them. 10 key genes (4CL: Sva01G002770,
   ANS: Sva06G008300, CAD: Sva09G007840, gibberellin 3beta-dioxygenase
   (GA3ox): Sva18G007280, RPM1-interacting protein 4 (RIN4): Sva16G021190,
   mitogen-activated protein kinase kinase kinase 2 (MAP3K2):
   Sva07G003520, peptide-methionine (R)-S-oxide reductase (msrB):
   Sva11G009210, glutathione peroxidase (GPX): Sva07G013260, tuliposide
   A-converting enzyme (TCEA): Sva10G009340, E3 ubiquitin-protein ligase
   RNF38/44 (RNF38_44): Sva01G007460), 7 TFs (AUX/IAA, bHLH, MIKC, MYB,
   NAC, ERF, and RLK), and 6 key metabolites (Syringin, Sinapic acid,
   Sinapoyl-CoA, Dihydroquercetin, Pelagoridin-5-glucoside-3-sambubioside,
   and Leucopelargonidin) were selected. Among these, 4CL and ANS were
   positively correlated with Syringin and negatively correlated with
   Dihydrokaempferol, Sinapic acid, and Leucopelargonidin. GA3ox, RIN4,
   MAP3K2, msrB, GPX, TCEA, RNF38_44, UX/IAA, bHLH, MIKC, MYB, NAC, ERF,
   and RLK were negatively correlated with Syringin, Sinapic acid,
   Sinapoyl-CoA, Dihydroquercetin, Pelagoridin-5-glucoside-3-sambubioside,
   and Leucopelargonidin (Fig. [84]5C). The results suggest that TFs and
   key genes may mediate the metabolism of phenylpropanoid and flavonoids,
   which further regulate flower development in S. variegata.

DEGs validation by qRT-PCR

   To check the accuracy of the differentially expressed genes, the
   expression of 19 representative DEGs was examined using quantitative
   real-time PCR (qRT-PCR) at three flowering stages of both male and
   female flowers of the S. variegata (Fig. [85]6, Supplementary Table
   S5). The DEGs were associated with phenylpropanoid and flavonoid
   biosynthesis in flower development. The relative expression levels of
   qRT-PCR and RNA-seq were logarithmized between different flowering
   stages. Correlation analysis was performed using the above logarithmic
   results, and the expression patterns were consistent with the RNA-seq
   data.

Fig. 6.

   [86]Fig. 6
   [87]Open in a new tab

   Validation of the expression of genes related to the development of
   male and female flowers in S. variegata. The expression levels of 18
   flower development-related genes were verified at three floral stages
   in male and female flowers using qRT-PCR. The relative expression
   levels of 18 DEGs were calculated using the 2^−ΔΔCt method; TIP41
   (Sva09G007630) served as an internal reference gene.

Discussion

   S. variegata, as a perennial shrub of the Salix genus, is a typical
   dioecious plant with good flooding, drought, and heavy metal pollution
   resistance, making it an excellent ecological restoration plant^[88]28.
   It grows mainly in the riparian zone and experiences the effects of
   water level subsidence in the riparian zone all year round. The
   flowering period of S. variegata is generally June- December, which
   overlaps with the time of flooding. However, it was found that S.
   variegata has a large number of seeds and its viability is unaffected.
   The molecular regulatory mechanisms of its floral ontogeny and
   development are worthy of our attention. In the last 30 years, many
   processes in flower development have been understood in detail through
   genetic and molecular analyses, and despite significant progress, many
   questions about the regulation of flower development in angiosperms
   remain unanswered^[89]29. In angiosperms, flowering is a trait of the
   plant that adapts to its environment and is important for its survival
   and reproduction^[90]30. The occurrence and development of floral
   organs is a complex process^[91]31. The complex flowering mechanism is
   not only regulated by many highly conserved genes but also influenced
   by external environmental factors such as photoperiod, temperature, and
   humidity^[92]32. Therefore, the data of a single omics cannot fully
   analyze the process of macro-development, and the multi-omics
   integration verification can make up for the missing data of a single
   omics analysis and reduce the false positives caused by a single omics
   analysis. In this study, we used the method of transcriptome and
   metabolome joint analysis to explore the regulation mechanism of S.
   variegata male and female flower development.

   We obtained a total of 12,245 DEGs by transcriptomic analysis, 20 DEGs
   shared among the three developmental stages of female flowers, F1 vs.
   F2, F1 vs. F3 and F2 vs. F3 with 1,281, 374, and 650 specifically
   expressed DEGs, respectively, and 383 DEGs shared among the three
   periods of male flowers, M1 vs. M2, M1 vs. M3 and M2 vs. M3 had 2,208,
   1,680 and 81 specifically expressed DEGs, respectively. 54 DEGs shared
   among the three comparison groups of male and female flowers at the
   same flowering stage, and 292, 3,646, and 1,300 specifically expressed
   DEGs in F1 vs. M1, F2 vs. M2, and F3 vs. M3, respectively. The highest
   number of DEGs were screened from bud to full bloom stage in both
   female and male flowers, suggesting the need for more genes to be
   involved in the S. variegata bud-to-flowering process. TFs play an
   important role in plant growth and development, not only can they be
   directly induced by signals related to life activities, but also
   through interactions with other proteins or direct regulation of the
   expression of downstream target genes, thus participating in the
   complex network of signaling and regulation corresponding to the growth
   and development of plants^[93]33. In this study, we identified and
   analyzed the highly expressed TFs during the development of male and
   female flowers of S. variegata, and found a large number of
   differentially expressed TFs, among which AUX/IAA, bHLH, MIKC, MYB,
   NAC, ERF, and RLK gene families were more highly expressed, which are
   involved in the differentiation of male and female flowers of S.
   variegata. The MADS-box gene family is involved in various stages of
   plant development, mainly promoting plant nutrient growth, determining
   the characteristics of floral organs, regulating the timing of
   flowering, promoting the formation of pollen and embryo sacs, and
   regulating the development of seeds and fruits^[94]34,[95]35. The
   mads-box gene family is divided into two types: type I (M_Type) and
   type (MIKC), which were originally identified as floral organ signature
   genes in Arabidopsis thaliana and Antirrhinum majus^[96]36. In this
   study, MIKC genes have important regulatory roles in the development of
   male and female flowers in S. variegata, with significant up- and
   down-regulation at different developmental periods. Previous studies
   have shown that MIKCc-type genes play a more important role in plant
   floral organ development^[97]37. In contrast, we similarly found a
   portion of the MIKC gene significantly up-regulated in both male and
   female flowers of S. variegata at the bud stage, and the full bloom
   stage and the final bloom stage reduced its expression. Auxin is an
   important hormone in the development process of plant flowers, which
   can affect the development of inflorescence and flower organs^[98]38.
   Auxin signaling consists of four core components: AUX/IAA, AUX/LAX,
   transport inhibitor response 1/auxin signaling F-box protein
   (TIR1/AFB), and ARF family member proteins^[99]39. It has been shown
   that members of the MYB gene family are key players in the late stages
   of floral organ development and maturation, involved in nectar
   development, stamen development, pistil length, secondary metabolites
   such as flavonols, phenylpropanoids, and odor
   compounds^[100]40–[101]43. ERF is one of the largest families of TFs in
   plants, and it plays a key role in transcriptional regulation during
   flower senescence^[102]44. It has been shown that TFs such as NAC are
   associated with senescence in several tissues of Arabidopsis thaliana,
   including petals and angiosperms^[103]45. The NAC family is the second
   largest family of transcription factors in plants, and NAC are involved
   in the regulation of a variety of processes such as flower development
   and growth hormone signaling^[104]46. NAC have strong activity during
   flower development, which also reflects that they may be one of the
   important regulators of male and female flower development in S.
   variegata. The bHLH TF family regulates CONSTANS in Arabidopsis, which
   is essential for photoperiodic flowering. In this study, a large number
   of bHLH genes were identified, which were highly expressed during
   different developmental periods of male and female flowers of S.
   variegata and may play an essential role in the regulation of the
   development of male and female flowers of S. variegata^[105]47–[106]49.
   RLK belongs to a large subfamily of proteins that can transmit ligand
   signals through autophosphorylation and initiate signaling
   cascades^[107]50. The RLK is involved in plant flower development,
   especially in pollen germination and pollen tube
   growth^[108]51,[109]52. In this study, the RLK genes were the most
   abundant and, as in previous studies, this class of transcription
   factors was highly expressed in male flowers of S. variegata, and some
   of the RLKs were significantly up-regulated for expression during the
   bud and full bloom stage in both male and female flowers^[110]53.

   Flower development is a process of floral transition from autotrophic
   to heterotrophic, during which floral secondary metabolism
   immediately^[111]49,[112]54. As expected, the process of pathways
   involved in secondary metabolite synthesis was significantly enriched.
   At this time, metabolite accumulation is related to traits such as
   color and smell, and these changes depend on the pollinators,
   precursors, etc^[113]55. Plants are able to synthesize a range of
   phenylpropanes and flavonoids that are important in plant growth,
   development, flowering, and pigmentation^[114]56. Anthocyanins are
   associated with flower color formation, and flavonoids and
   phenylpropanoid compounds are co-pigments that stabilize anthocyanins,
   and both metabolites are synthesized via the phenylpropane-related
   pathway^[115]57. It has been demonstrated that flavonoids and
   phenylpropanoid-related genes are highly expressed during the flowering
   stage of plants^[116]58. In our study, the stigma of female flowers of
   S. variegata slowly changes from golden yellow to brown after
   flowering, and the anthers of male flowers are also golden yellow and
   gradually turn brown, which may also be caused by the gradual
   accumulation of phenylpropanoid and flavonoids and their derivatives.
   Similarly, flavonoids and phenylpropanoid-related genes were
   significantly higher in expression during the full bloom stage and the
   final bloom stage in both female and male flowers of S. variegata.
   Consistent with the changes in gene expression, the flavonoid and
   phenylpropane biosynthetic pathways were also enriched from the flower
   bud stage to the final flowering stage. To attract pollinators during
   flowering, plants often produce volatiles and pigments, including
   phenylpropanoids and flavonoids. The increased accumulation of these
   two types of substances during the development of male and female
   flowers of S. variegata may also serve to attract insect pollinators.
   These results are similar to previous studies on flower development in
   Lonicera japonica Thunb, Rosmarinus officinalis L., and Pogostemon
   cablin (Blanco) Benth., among others^[117]32,[118]58,[119]59. It has
   been shown that the synthesis of flavonoids often requires the
   production of precursors via the phenylpropane pathway, and these
   processes involve many important enzymes such as 4CL, ANS,
   F5H^[120]60,[121]61. Some studies have found that some related genes
   regulating plant hormones, such as GA signaling-related genes and
   growth hormone signaling-related genes, are negatively correlated with
   flavonoids and phenylpropanoids. In this study, we found that
   phytohormone-related signaling pathways were also enriched during the
   development of male and female flowers of S. variegata, suggesting that
   GA and growth hormone signaling may be involved in the metabolism and
   biosynthesis of flavonoids, phenylpropanoids, and their
   derivatives^[122]62. GA3ox is a target gene for GA biosynthesis, which
   can affect stamen development. Consistent with previous studies, the
   gene GA3ox was significantly up-regulated in male and female flowers
   during the development of male and female flowers in S. variegata. Many
   studies have shown that MAP3K, RIN4, regulates various physiological
   and biochemical processes throughout the plant growth cycle, including
   germ cell and embryo development, flower development^[123]63,[124]64.
   The msrB is an active enzyme in plants that is more abundant in plant
   developing organs^[125]65. TCEA, RNF38_44 are also active substances in
   plants that influences growth and development and enhances resistance
   in plants^[126]66,[127]67. The GPX plays a crucial role in inducing
   oxidative stress and promoting plant cell growth^[128]68. Functional
   clustering of these genes, which we found to be significantly
   up-regulated and mostly associated with flavonoid and phenylpropanoid
   synthesis during the developmental stages of both male and female
   flowers of S. variegata, further highlighted the key biological
   processes and molecular functions related to flavonoid production. In
   addition, KEGG enrichment analysis revealed that these genes were
   associated with flavonoid and phenylpropanoid metabolism. Through
   combined transcriptomic and metabolomic analyses, we found that TFs
   could regulate the biosynthesis of flavonoids and phenylpropanes by
   regulating their biosynthesis-related genes. The changes of DEMs and
   DEGs in these pathways were similar during the different developmental
   periods of S. variegata male and female flowers. S. variegata in both
   male and female flowers were similar in terms of their trends in
   different developmental periods. This suggests that the above genes
   play an important role in regulating the accumulation of this
   metabolite in S. variegata male and female flowers at different
   developmental periods.

Conclusion

   In this study, we used the male and female flowers of S. variegata at
   three developmental stages as materials and combined the transcriptome
   and metabolites to focus on the regulation of key metabolites that may
   be involved in the DEGs during the development of the flowers of this
   species. The study of the molecular mechanisms underlying the
   development of the flowers of S. variegata, which is a typical
   dioecious plant with important ecological restoration and levee-control
   plant, will enable us to better understand the adaptive evolution of
   its adaptive evolution of the breeding system. The results of these
   studies can also provide candidate genes and the theoretical basis of
   metabolism for the further utilization of molecular biology to regulate
   the accumulation of secondary metabolites in male and female flowers of
   S. variegata. Our results will contribute to understanding of the
   molecular and genetic mechanisms of male and female flower development
   in typical dioecious plants.

Materials and methods

Materials and pre-treatment

   The samples for this study were collected from the riverside of
   Jingangbei, Beibei District, Chongqing, China (106.42092155°E,
   29.84939019°N). The plant species was identified by Professor Hongping
   Deng. It’s voucher specimen was deposited in the botanical herbarium of
   Southwest University (Holotype: SWCTU! ) under voucher number
   Sva2023081901. Based on the bloom stage of S. variegata male and female
   flowers, samples of flower buds (male and female flowers were wrapped
   in the bud, F1/M1), full bloom stage (the stigmas / filaments were
   fully unfolded, F2/M2), final bloom stage (stigma pollination was
   completed / anther powder finished) were collected in July 2023
   (Supplementary Fig. S5). At each developmental stage of the
   investigated female and male flowers, three replicates each were taken
   immediately frozen in liquid nitrogen and stored at -80 °C for
   transcriptome analysis by RNA-seq and six replicates were stored at
   -80 °C for metabolite identification (Supplementary Table S5-S6). The
   plants collected in this study and the subsequent experiments performed
   were in line with the IUCN Policy Statement on Research Involving
   Species at Risk of Extinction and the Convention on the Study Trade in
   Endangered Species of Wild Fauna and Flora. The plants were collected
   strictly following the provisions of the Regulations of the People’s
   Republic of China on the Protection of Wild Plants, and a collection
   permit was obtained.

Transcriptomic analysis

   A cDNA library from three developmental stages of male and female
   flowers of S. variegata was constructed for RNA-seq. The DNBSEQ-T24
   sequencing platform was utilized for sequencing. The S. variegata
   genome was selected as a reference (Genome data are available in the
   National Genomics Data Center with the project number PRJCA026796).
   Gene expression levels were quantified in terms of exon model per
   kilobase per million mapped fragments (FPKM) values. DEGs were screened
   using DESeq2 and adjusted p-value and fold change. DEG enrichment was
   assessed based on comparisons with Gene Ontology (GO) and KEGG
   data^[129]69. TFs were annotated using iTAK software and all TFs were
   obtained by BLAST search using Plant TFDB.

Metabolomics analysis

   Untargeted metabolites were determined in male and female flowers of S.
   variegata using liquid chromatography-tandem mass spectrometry
   (LC-MS/MS) platform (Biomarker.

   Technologies Co., Ltd.). The powder flower samples were weighed into
   50 mg, and 1000 µL of the extraction liquid (methanol: acetonitrile:
   water v/v/v = 2:2:1) containing 20 mg/L internal standard
   (2-Chloro-L-phenylalanine) was added, vortexed, and mixed for 30 s. The
   samples were ground and sonicated in an ice-water bath for 10 min and
   then allowed to stand at -20 ℃ for 1 h. The results were summarized as
   follows: (a) The extract was analyzed by LC-MS/MS. After multiple
   centrifugations, the supernatant was taken to mix the QC samples in
   preparation for online testing. LC-MS on a Waters UPLC Acquity I-Class
   PLUS (Waters Corporation, Milford, CT, USA) and a Waters UPLC Xevo
   G2-XS QTOF (Waters Corporation, Milford, CT, USA) were used to analyze
   the sample extracts. The detected peak area information was normalized
   for subsequent analysis. The reproducibility of within-group and
   quantitative control samples was examined using principal component
   analysis. The identified compounds were categorized and pathway
   information was retrieved in databases such as KEGG, HMDB, and
   lipidmaps, the multiplicity of differences was calculated and compared
   based on the grouping information, and the t-test was utilized to
   calculate the significance of the difference p-value for each compound.
   The OPLS-DA modeling was performed using the R language package ropes,
   and the reliability of the model was verified by 200 substitution
   tests. DEMs were screened by combining the multiplicity of differences,
   p-value, and VIP value of the OPLS-DA model. The screening criteria
   were FC ≥ 1, P-value ≤ 0.01, and VIP ≥ 1. A hypergeometric distribution
   test was utilized to calculate the differential metabolites with KEGG
   pathway enrichment significance.

Integrated transcriptome and metabolome analysis

   Transcriptome and metabolome integrated correlation analysis was
   performed using using MetaboAnalyst 5.0^[130]70. WGCNA analysis of
   genes associated with important DEMs was performed using the WCNA
   package in R software^[131]71. A weighted neighborhood connectivity
   matrix between different genes was constructed using power functions,
   and similar modules in the hierarchical tree were filtered using a
   dynamic tree-cutting procedure. Pearson correlation analysis of DEGs
   and DEMs was performed based on the R package. Gene and metabolite
   relationship networks were visualized by Cytoscape_v3.9.1.

Validation of transcriptomic data

   To examine the validity of the transcriptome data, 18 DEGs were
   randomly selected for quantification at three flowering stages in both
   male and female flowers of S. variegata. Total RNA was extracted using
   the TRlzol kit (Invitrogen, CA, USA) and cDNA was synthesized using
   Goldenstar RT6 cDNA Synthesis Mix (TSK314S) according to the maker’s
   instructions. The type 2 A phosphatase activator TIP41 (TIP41:
   Sva09G007630) gene was used as an internal control. All experiments
   were repeated at least 3 times.

Electronic supplementary material

   Below is the link to the electronic supplementary material.
   [132]Supplementary Material 1^ (2.5MB, xlsx)
   [133]Supplementary Material 2^ (882.7KB, docx)

Acknowledgements