
Early stage pistil degeneration in male flowers
To investigate the dynamic gene expression during the development of reproductive organs in female and male papaya flowers, we collected 17 samples from the pistils of female flowers (FP), rudimentary pistils of male flowers (MP), and stamens of male flowers (MS) across five developmental stages, as well as from the apical meristems of female (FAM) and male plants (MAM). Each sample had three biological replicates and was subjected to RNA sequencing (RNA-seq). The morphological characteristics of female flowers with a functional gynoecium and male flowers with rudimentary pistils and functional stamens were observed at five stages (Fig. 1a, Table S1). Morphological observations revealed that pistil degeneration occurs at a very early stage in male flowers, confirming the classification of male flowers as type I unisexual flowers.
Morphological characteristics of papaya floral organs and dynamic expression profiles at different developmental stages. (a) Morphological features of female and male papaya plants, and female pistils (FP), aborted pistils in male flowers (MP), and stamens in male flowers (MS) at five developmental stages. Flower lengths are 1.5–2.5 mm, 2.5–4.5 mm, 4.5–6 mm, 7–12 mm and 15–20 mm from Stage1 to Stage5, respectively. (b) Number of differentially expressed genes (DEGs) in apical meristem (AM) and stages 2–5 compared to stage 1 across female pistils (FP), aborted pistils in male flowers (MP), and stamens (MS), totaling 15 groups. (c) Venn diagram of DEGs across five DEG sets in female pistils, aborted pistils in male flowers, and stamens. FAM and MAM represent the apical meristems of female and male papaya plants, respectively. (d) Enrichment of plant hormone signal transduction pathways across 15 different combinations. For each group, the top 25 pathways with p-value < 0.05 were selected for statistical analysis
To explore the relationships among the 17 samples, we calculated their correlation coefficients. Correlation analysis of the RNA-seq datasets indicated a strong linear correlation between different floral organs during the early stages of flower development. Notably, the third (M3S) and fourth (M4S) stages of stamen development in male flowers generally showed lower correlations with other developmental stages (Fig. S1a). Conversely, a strong correlation was observed among the different stages of aborted pistil development in male flowers. Principal component analysis (PCA) of the RNA-seq datasets of the 17 samples (including three replicates for each sample) showed clear distinctions at different stages and a compact distribution of biological replicates (Fig. S1b-S1d). In the PCA of the MP, PC1 and PC2 explained 39.84% and 18.37% of the variation, respectively (Fig. S1c). The proximity of sample distributions across different groups reflected the similarity in gene expression patterns between stages (Fig. S1c), which is consistent with previous correlation analyses. Thus, the early determination of aborted pistils in male flowers causes functional loss during development, resulting in the non-specific expression of related functional genes.
GO and KEGG enrichment analysis reveals the critical role of plant hormone signaling pathways in papaya floral organ development
To understand differential gene expression, we divided the transcriptome datasets of the 17 samples into three groups based on the floral organ type: functional pistils, rudimentary pistils, and stamens. Each group used the earliest stage (the first period) as a control for differential analysis (Fig. 1b and 1c, Supplementary Data 1). The number of differentially expressed genes (DEGs) in the pistils of female flowers increased to 346 (F2P vs. F1P), 1,279 (F3P vs. F1P), 2,344 (F4P vs. F1P), 2,416 (F5P vs. F1P), and 1,518 (FAM vs. F1P). In the rudimentary pistils of male flowers, the number of DEGs was 218 (M2P vs. M1P), 1,298 (M3P vs. M1P), 3,109 (M4P vs. M1P), 3,817 (M5P vs. M1P), and 2,565 (MAM vs. M1P). The number of DEGs during different pistil stages of female flowers was similar to that of the rudimentary pistils of male flowers (Fig. 1b, Table S2). The number of DEGs increased progressively during the development of pistils in female flowers, rudimentary pistils in male flowers, and stamens in male flowers, indicating clear differences in gene expression patterns between the early and late stages of floral organ development, with only a few key genes being likely involved in regulating floral organ identity in the early stages. Notably, the number of DEGs increased dramatically during the middle to late stages of stamen development (M3S-M5S) (Fig. 1b, Table S2), suggesting that stamen development in papaya may be subjected to more extensive gene regulation. Among the female pistils, aborted pistils in male flowers, and stamens, there were 117, 38, and 127 DEGs, respectively, which were common to all five DEG groups (Fig. 1c). This indicates that these genes are differentially expressed throughout the developmental period and likely play crucial roles in the overall development of the pistils and stamens. The number of common DEGs (38) was lower in aborted pistils of male flowers than in the other two organs (117 and 127). Moreover, the DEGs between different developmental stages gradually increased, starting from the first stage, in the normal pistils of female flowers and abortive pistils of male flowers. However, the DGEs in MAM sharply decreased (MAM vs. M1S) compared to M5S vs. M1S (Fig. 1c). Additionally, in the three comparisons, only a small number of stage-specific DEGs were present in stage2 vs. stage1 comparisons (14 DEGs in F2P vs. F1P, 18 DEGs in M2P vs. M1P, and 7 DEGs in M2S vs. M1S) (Fig. 1c).
Furthermore, Gene Ontology (GO) enrichment analysis of DEGs in the pistils of female flowers (FPs vs. F1P) revealed significant enrichment in the top 25 GO terms related to responses to biotic and abiotic stimuli and hormone responses throughout pistil development (Fig. S2a-S2e). Specifically, there were 46, 158, 236, 218 DEGs in the F2P vs. F1P, F3P vs. F1P, F4P vs. F1P and F5P vs. F1P groups enriched in the GO term “response to hormone.” The number of DEGs markedly increased during the F2P developmental stages, highlighting the critical regulatory role of hormones in early pistil development. In contrast, in the later stages (F4P vs. F1P and F5P vs. F1P), 188 and 213 DEGs were significantly enriched in response to light stimulus, respectively (Fig. S2d-S2e). In the rudimentary pistils of male flowers (MPs vs. M1P), DEGs showed significant enrichment in hormone responses, which gradually increased over time (Fig. S3a-S3e). Unlike the pistils of female flowers, DEGs during the early development of rudimentary pistils in male flowers (M2P vs. M1P group) were significantly enriched in GO terms related to the transport of energy substances, including sugars and amino acids, and the response to abscisic acid (Fig. S2b and S3b). From the M2P stage, GO terms related to light stimulation or photosystem (including GO:0009416, GO:0019684, and GO:0042550) were significantly enriched (Fig. S3b-S3e). In the stamens of male flowers, DEGs during the early stages of development (M2S vs. M1S and M3S vs. M1S) were significantly enriched in GO terms related to pollen tissue structure and shoot system development, including pollen wall assembly, pollen exine formation, and shoot system development (Fig. S4b-S4c). In addition, single-organism processes were significantly enriched throughout stamen development (Fig. S4a-S4e).
Additionally, we conducted Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis on the aforementioned sets of DEGs (FPs vs. F1P, MPs vs. M1P, and MSs vs. M1S) and selected the top 30 KEGG pathways from each group for further analysis. These results suggest that the role of the plant hormone signal transduction pathway varies among female pistils, aborted pistils in male flowers, and stamens (Fig. 1d). In female pistils, the plant hormone signal transduction pathway was significantly enriched throughout development. As the developmental stages progressed, the number of DEGs in the pathway increased: 8 DEGs in F2P vs. F1P, 28 in F3P vs. F1P, 43 in F4P vs. F1P, and 39 in F5P vs. F1P, indicating its critical role in pistil development. The enrichment trend in the rudimentary pistils of male flowers was similar to that in the pistils of female flowers. However, no enrichment was observed in the M2P vs. M1P group. In the stamens of male flowers, significant enrichment was observed only in the first two comparisons (M2S vs. M1S, and M3S vs. M1S). Significant differences were observed in the enrichment of DEGs in plant hormone signal transduction pathways during the development of different floral organs, suggesting different roles of plant hormones during stamen and pistil development.
Next, we merged and removed redundancies from the DEGs enriched in plant hormone signal transduction pathways in female pistils, aborted pistils in male flowers, and stamens, resulting in three sets of DEGs specific to each tissue. A Venn diagram was constructed using these three DEG sets, which identified 8, 24, and 20 genes specific to female pistils, aborted pistils, and stamens, respectively (Fig. S1e). Based on the expression patterns and potential functions of homologous genes, we identified three highly expressed genes in the female pistil-specific enriched gene set: one involved in cytokinin regulation and two in gibberellin regulation. In the aborted-pistil-specific gene set of male flowers, five auxin-related genes were identified, all of which showed lower expression levels in aborted pistils than in normal pistils. Additionally, five genes involved in cytokinin regulation were identified in the stamen-specific gene set, all of which showed higher expression in stamens. The differential enrichment of plant hormone signal transduction pathways across female flowers, aborted pistils, and male flowers, along with the specific enrichment of hormone-regulating genes, suggests that hormone pathways are potentially involved in the regulation of sex differentiation in papaya, with certain hormones potentially playing key roles in this process.
Identification of MADS-box gene family and ABCDE subfamily members in papaya
MADS-box members are predominantly involved in floral organ development and play critical roles in floral organ formation. Using the MADS-box protein sequences of Arabidopsis thaliana and Vitis vinifera, we identified 65 MADS-box gene family members in papaya (Table S3). To further understand the evolutionary relationships among members of the MADS-box gene family, we constructed a phylogenetic tree combining the MADS-box gene families of the three aforementioned species (Fig. 2a). In papaya, MIKCC subfamily members (29 genes) and M-type subfamily members (30 genes) were almost equally abundant, whereas in Arabidopsis, MIKCC subfamily members were fewer than in the M-type members.

Phylogenetic analysis of MADS-box transcription factor gene family members and identification of ABCDE homologs in papaya. (a) Phylogenetic tree constructed using Arabidopsis (103 members), grape (54 members), and 65 MADS-box members identified in this study, with 1000 bootstrap replicates, using the neighbor-joining (NJ) method. (b) Identification of ABCDE homologous genes in papaya using Arabidopsis ABCDE genes as references, and construction of a phylogenetic tree. (c) Schematic diagram of the ABCDE flower development model and heatmap of the expression patterns of the 10 identified ABCDE subfamily MADS-box genes in papaya. In addition to RNA-seq data from the apical meristems of female (FAM) and male plants (MAM), female pistils (FP), aborted pistils in male flowers (MP), and stamens (MS) used in this study, RNA-seq data from female (F_f), male (M_f), and hermaphroditic flowers (H_f) were included to show the expression patterns of the ABCDE subfamily MADS-box genes. The TPM values were log₂-transformed before being used as input data for heatmap visualization
The phylogenetic tree was distinctly separated into five evolutionary clades: M-alpha, M-beta, M-gamma, MIKC*, and MIKCC (Fig. 2a, Table S4). A recent study suggested that the monophyletic phylogeny that groups the previously defined Type II MIKC* subfamily with the Type I M-type subfamily is significantly more robust than the previously proposed hypotheses of monophyletic MIKC* and MIKCC division [48]. Therefore, classifying MIKC* members as part of the Type I group rather than Type II is recommended. In the present study, we classified 65 MADS-box members in papaya according to this revised model, assigning them to Type I (including M-alpha, M-beta, M-gamma, and MIKC*) and Type II (containing only MIKCC). The gene structure, conserved motifs, and domains of the MADS-box genes in papaya were analyzed to further understand the MADS-box genes in papaya. MEME predicted 20 motifs, of which motifs 1 and 10 were widely distributed among the 65 MADS-box members, whereas motif 15 appeared to be specific to Type II members (Fig. S5b). Analysis of the conserved protein domains revealed that all Type II genes contained the MADS_MEF2_like domain, and the K-box domain was exclusively present in Type II genes (Fig. S5c). Gene structure analysis indicated clear differences between the M-type and MIKCC-type members. Most M-type members contain a single exon, whereas the MIKCC-type members typically exhibit longer gene structures with numerous introns (Fig. S5d).
MADS-box genes are widely involved in the morphogenesis and development of floral organs, particularly the Type II members. Following the establishment of the ABCDE flower development model, the MADS-box gene family is recognized as essential for floral organ development. Moreover, in the ABCDE model of flower development, all MADS-box genes belong to the Type II group. In our study, we identified 10 ABCDE homologous genes from 65 MADS-box gene family members in papaya, including CpAPETALA1 (CpAP1), CpFRUITFULL (CpFUL), CpAPETALA3.1 (CpAP3.1), APETALA3.2 (CpAP3.2), CpPISTILLATA (CpPI), CpPLENA (CpPLE), CpSEEDSTICK (CpSTK), CpSEPALLATA1/2 (CpSEP1/2), CpSEPALLATA3 (CpSEP3), and CpSEPALLATA4 (CpSEP4) (Fig. 2b, Table S5).
To further investigate the transcript levels of ABCDE genes in papaya floral organs, we analyzed the expression patterns of these genes at different developmental stages in female pistils, abortive pistils in male flowers, and stamens (Fig. 2c, Table S6). The A-class gene CpAP1 exhibited a broad expression pattern with high expression levels in female, male, and hermaphroditic flowers. Its expression is also elevated during the early stages of pistil development and in aborted pistils in male flowers and stamens. B-class genes, including CpAP3.1, CpAP3.2, and CpPI, were highly expressed in the stamens of male papaya flowers, with sharp upregulation during the later stages of stamen development (M4S and M5S). In addition, B-class genes were more highly expressed in male and hermaphroditic flowers than in female flowers. The expression patterns of the CD class genes CpSTK and CpPLE were notably different. CpPLE was broadly expressed throughout the development of female pistils, rudimentary pistils in male flowers, and stamens, with gradually increasing expression levels. It was also expressed at higher levels in male and hermaphroditic flowers than in female flowers. In contrast, CpSTK was almost exclusively expressed in the female pistils, which is consistent with the results of recent studies [24, 64, 65], with its expression levels sharply increasing during development. The E-class genes, CpSEP1/2, CpSEP3, and CpSEP4, were also broadly expressed throughout the development of female pistils, rudimentary pistils, and stamens of male flowers. Notably, CpSEP1/2 was highly expressed in the rudimentary stamens of male flowers, with gradual upregulation. However, SEP genes show differential expression patterns between flowers of different sexes in papaya, with higher expression observed in male and hermaphroditic flowers, particularly pronounced for CpSEP1/2 and CpSEP3.
In summary, the expression patterns of the ABCDE genes in papaya strongly support the ABCDE model of floral development, providing substantial evidence for its applicability to angiosperms, especially core eudicots.
Evolutionary insights to ABCDE subfamily MADS-box genes reveal gene contraction and expansion in papaya
To further understand the evolution of the ABCDE subfamily of MADS-box genes, we selected the ABCDE subfamily of MADS-box genes and AGL6 from gymnosperms (Ginkgo biloba, Gnetum gnemon, Picea abies, and Cycas panzhihuaensis), basal angiosperms (Amborella trichopoda and Nymphaea colorata), monocots (Oryza sativa and Zea mays), and eudicots (Arabidopsis thaliana and Carica papaya) to construct a phylogenetic tree. The phylogenetic tree clearly clustered 100 ABCDE and AGL6 proteins from 10 seed plants into five subgroups: AP1/FUL (A-class genes, 17), AP3 and PI (B-class genes, 20), AGAMOUS (CD- or C-class genes, 24), SEP (E-class genes, 24), and AGAMOUS-LIKE 6 (AGL6 genes, 15) (Fig. 3a, Table S5).

Evolutionary analysis of ABCDE and AGL6 genes in 10 seed plants. (a) Statistics of A, B, CD, E, and AGL6 class genes in gymnosperms (Cycas panzhihuaensis, Ginkgo biloba, Gnetum gnemon and Picea abies), basal angiosperms (Amborella trichopoda, and Nymphaea colorata), monocots (rice and maize), and dicots (Arabidopsis and papaya). (b) Phylogenetic tree of ABCDE and AGL6 genes in 10 seed plants. Best-fit model: JTT + F + R6, constructed with 1000 bootstrap replicates using the maximum likelihood (ML) method
The A-class genes are involved in the development of sepals and petals. Our results indicate that A-class genes are specific to angiosperms, which is consistent with previous reports [49, 66, 67]. Additionally, in monocots, A-class genes seem to have undergone duplication, resulting in more paralogous genes than in dicots (Fig. 3a and 3b, Table S5). The emergence and subsequent specialization of A-class genes in angiosperms may have substantially contributed to the diversity of floral organs in these plants.
B-class genes control the development of petals and stamens. In gymnosperms, this subfamily of genes has been identified in Ginkgo biloba (two B-class genes), Cycas panzhihuaensis, and Gnetum gnemon (one B-class gene), but not in Picea abies. In angiosperms, two B-class genes were identified in Nymphaea colorata, Zea mays, and Arabidopsis thaliana, three B-class genes were obtained in Oryza sativa and Carica papaya; and Amborella trichopoda had four B-class genes, indicating a clear gene duplication event (Fig. 3a and 3b, Table S5). In the present study, B-class genes were generally more abundant in angiosperms than in gymnosperms. Notably, in papaya, despite the relatively small number of genes, two paralogous genes, CpAP3.1 and CpAP3.2, were present within the AP3 gene family. The greater abundance of B-class genes in flowering plants likely contributes to the formation of more complex floral organs, which is consistent with the conclusions of Shen et al. [49].
CD-class genes are involved in controlling stamen and pistil development. In the phylogenetic tree, the CD-class genes formed a distinct clade with 96% bootstrap support, indicating high confidence in this branch (Fig. 3b). The number of CD-class genes is higher in angiosperms than gymnosperms, suggesting that the increase in CD-class genes may be closely related to the evolution of complex floral organs in angiosperms. Notably, in papaya, we identified only two CD-class genes, CpSTK and CpPLE, which distinctly differed from other angiosperms such as rice with four CD-class genes, maize with five, and Arabidopsis with four (Fig. 3a, Table S5). In the more primitive basal angiosperms Amborella and Nymphaea, two and three CD class genes were observed, respectively, with Nymphaea having two copies of the AG gene (Table S5).
E-class genes interact with ABCD genes [68] and are involved in the development of floral organs. Phylogenetic analysis suggests that E-class genes, similar to A-class genes, are specific to angiosperms. Further analysis revealed that SEP1 and SEP2 were more closely related to SEP4 than to SEP3, which is consistent with previous research [49, 69]. Additionally, in monocots, the LOFSEP subclass of E class genes (also known as the SEP1/2/4 clade) formed a distinct branch, separating E class genes into two clades, along with those from core eudicots and basal angiosperms, with a bootstrap support of 95% (Fig. 3b). Phylogenetic analysis indicated that the LOFSEP subclass genes in monocots were more specialized than the E-class genes in the core eudicots and basal angiosperms (Fig. 3b). The emergence of these new organ forms provide protection to monocot flowers and represents an evolutionary advancement in floral organ development.
AGL6 is involved in floral development in angiosperms [38, 70] and cone formation in gymnosperms [71]. We identified AGL6 across all gymnosperm samples (Fig. 3a and 3b, Table S5). Additionally, multiple copies of AGL6 were present in Ginkgo biloba, Gnetum gnemon, and Picea abies, highlighting the importance of AGL6 in the development of gymnosperm reproductive organs. The phylogenetic tree showed a clear separation between the AGL6 subfamilies in angiosperms and gymnosperms, with the angiosperm clade exhibiting longer branches, suggesting potential evolutionary diversification of the AGL6 gene in angiosperms (Fig. 3b).
In our phylogenetic analysis, the MADS-box genes of the ABCDE subfamily and AGL6 were clearly distinguished in both gymnosperms and angiosperms. In particular, as a core eudicot, the ABCDE and AGL6 genes showed the closest phylogenetic relationships with those of Arabidopsis, which is a core eudicot. This finding is consistent with the taxonomic classifications and has been validated using molecular evolutionary analyses [9, 62]. Additionally, we observed that CD-class and E-class genes within the ABCDE subfamily, appeared to be lost in papaya, in contrast to other angiosperms, especially grasses (Poaceae) [38].
Clustering analysis of dynamic gene expression unveils flower organ-specific clusters in papaya pistil and stamen development
To explore the dynamic expression changes in floral organs across different developmental stages, we performed cluster analysis on the mRNA datasets. Cluster analysis of DEGs from F5P vs. F1P, M5P vs. M1P, and M5S vs. M1S was conducted using the Gap Statistic method to determine optimal clusters (Fig. S6a-S6c). We identified 12 clusters in the female pistils (Fig. S7, Table S7), 19 clusters in male abortive pistils (Fig. S8, Table S8) and 23 clusters in male stamens (Fig. S9, Table S9).
In female pistils, cluster 8 genes showed increased expression over time, whereas their expression remained low in male abortive pistils (Fig. 4a and 4b). Expression pattern and functional annotation analyses revealed that many genes involved in plant hormone signal transduction, especially auxin-related genes such as CpIAA32 (INDOLE-3-ACETIC ACID INDUCIBLE 32), CpGH3.1, and CpEOD1 (ENHANCER1 OF DA1), were enriched in cluster 8 (Fig. 4d, Table S10). In contrast, cluster 16 genes in male abortive pistils exhibited the opposite expression patterns (Fig. 4d and 4e), indicating their potential role in the negative regulation of female reproductive function. Several genes involved in plant hormone regulation were also identified in cluster 16. For instance, CpLOG5 (LONELY GUY 5) and CpCKX1 (CYTOKININ OXIDASE/DEHYDROGENASE 1) are directly involved in the synthesis and degradation pathways of cytokinins, while CpGASA1 (GAST1 PROTEIN HOMOLOG 1), CpSCL15 (SCARECROW-LIKE 15), and CpARF4 (AUXIN RESPONSE FACTOR 4) participate in hormone signaling pathways such as gibberellin and auxin. Additionally, many genes involved in oxidative stress response and sucrose transport were enriched in cluster 16, including CpMT2A (METALLOTHIONEIN 2A), CpPER12 (Peroxidase 12), CpLETM2 (LEUCINE ZIPPER-EF-HAND-CONTAINING TRANSMEMBRANE PROTEIN 2), CpF3H (FLAVANONE 3-HYDROXYLASE), CpAPX1 (ASCORBATE PEROXIDASE 1), CpUGT84A1, and CpSTP8 (SUGAR TRANSPORT PROTEIN 8) (Fig. 4f, Table S10).

Expression pattern and enrichment analysis of clusters with specific expression in female pistils and aborted pistils in male flowers. (a–b) Dynamic expression trends of DEGs in cluster 8 in female pistils. (c) Schematic of papaya female flower and expression patterns of genes specifically expressed in cluster 8. (d–e) Dynamic expression trends of DEGs in cluster 16 in aborted pistils in male flowers. (f) Schematic of papaya male flower and expression patterns of genes specifically expressed in cluster 16. The TPM values were row-normalized before being used as input data for heatmap visualization. (g) GO enrichment analysis for DEGs in cluster 8 in FP and in cluster 16 in MP. (h) KEGG pathway enrichment analysis for DEGs i in cluster 8 in FP and in cluster 16 in MP. Top 25 GO terms or KEGG pathways are shown
Moreover, cluster 8 genes in female pistils revealed significant enrichment in carpel and organ morphogenesis (Fig. 4g), as well as in plant hormone signal transduction and diterpenoid biosynthesis pathways (Fig. 4h). This suggests their role in the occurrence and formation of female organs, or the execution of female-specific functions. Conversely, in cluster 16, male abortive pistil genes were enriched in phloem development and amino acid transport (Fig. 4g) as well as in secondary metabolism and starch and sucrose metabolic pathways (Fig. 4h).
Similarly, cluster 20 genes were upregulated only during the early stages of stamen development and peaked in the middle stage (M3S) in the stamens of male flowers (Fig. S10a). Moreover, many genes were involved in the formation of the pollen outer wall, glycosyl and lipid hydrolysis, and transport were found in cluster 20, including CpACOS5 (ACYL-COA SYNTHETASE 5), CpSWEET15, CpENGASE85B (ENDO-BETA-N-ACETYGLUCOSAMINIDASE 85B), CpGELP77 (GDSL-TYPE ESTERASE/LIPASE 77), CpLPTG10 (GLYCOSYLPHOSPHATIDYLINOSITOL-ANCHORED LIPID PROTEIN TRANSFER 10), and CpCYP704B1 (CYTOCHROME P450) (Fig. S10b, Table S10). These genes are closely related to cell wall remodeling, which is crucial for floral organ morphogenesis. Additionally, most genes in Arabidopsis are crucial for stamen and pollen development and are essential for maintaining male reproduction [72,73,74]. Enrichment analysis showed that the genes in cluster 20 were significantly enriched in GO terms such as anther morphogenesis and response to external stimuli (Fig. S10c and S10d), suggesting their potential role in promoting stamen morphogenesis and anther maturation.
Gene regulatory networks revealed the key roles of CD-class and E-class subfamily MADS-box genes in papaya pistil development
To preliminarily analyze the regulatory network of pistil and stamen development in papaya, we constructed gene regulatory networks (GRNs) for floral organ development. GRNs for female pistils (using DEGs of F5P versus F1P), male abortive pistils (using DEGs of M5P versus M1P), and male stamens (using DEGs of M5S versus M1S) were constructed using Mutual Rank (MR) analysis.
The GRN for female pistils contained 1,841 potential regulatory relationships, involving 1,063 genes (Table S11). Analysis of the gene expression patterns within the network revealed a gradual increase in the expression levels of CpSTK and CpSUP (SUPERMAN) in the pistils of female flowers, whereas they were almost unexpressed in the abortive pistils and stamens of male flowers (Fig. S11a and Table S12). This pistil-specific expression pattern suggests that CpSTK and CpSUP play important roles in pistil development and morphogenesis. In Arabidopsis, STK encodes a MADS-box transcription factor that is expressed within ovules [64, 65]. STK is widely recognized as a D-class gene that specifies ovule identity [46, 75] and is essential for ovule specification. Additionally, CpSTK and CpSUP both belonged to cluster 8 in the pistil of female flowers (Fig. 4a and 4c). Therefore, we extracted a subnetwork with CpSTK and CpSUP as core genes and constructed a pistil regulatory network (PRN) for pistil development (Fig. 5a). Analyzing the gene expression patterns in the PRN revealed that CpIND (INDEHISCENT) and CpZFP11 show expression patterns similar to those of the core genes (Fig. S11a, Table S12). Both ZFP11 and SUP are members of the zinc finger protein family, and studies in Arabidopsis revealed that the ERF repressor domain is essential for ZFP11 to function as a transcriptional repressor, consistent with SUP [76]. IND is the closest paralog of the HEC gene in Arabidopsis. Although IND mainly functions in the development of the fruit dehiscence zone, ind mutants exhibit minor defects in stigma elongation [77, 78], suggesting the role of CpZFP11 and CpIND in papaya pistil development and are regulated by core genes (CpSTK and CpSUP) within the network. Enrichment analysis of genes in the PRN highlighted significant enrichment in the biological process (BP) category, particularly in GO terms associated with flower organ formation, floral organ identity specification, and anatomical boundary formation (Fig. S11b, Table S17). This further suggests that the genes in the subnetwork centered on CpSTK and CpSUP are closely related to the initiation and development of floral organs and likely play important roles in pistil development.

Construction of regulatory networks for the development of female pistils (FP), aborted pistils in male flowers (MP), and stamens (MS). (a) Construction of a regulatory network for female pistil development (PRN) with CpSTK and CpSUP as core genes. (b) Construction of a regulatory network for aborted pistil development in male flowers (APRN) with CpSEP1/2 as core genes. (c) Construction of a regulatory network for stamen development (SRN) with CpTDF1 as the core gene. Genes in each network were classified based on significant GO term enrichment results. Genes within boxed regions in the network are enriched in the same GO term
In male abortive pistils, the GRN included 6,882 potential regulatory relationships involving 2,609 genes (Table S13). Further analysis revealed that the expression of CpSEP1/2, an E-class gene, markedly increased during the development of abortive pistils in male flowers (Fig. S11c, Table S14). Flowering plants require angiosperm-specific E-class genes to facilitate the interaction of B- and C-class genes, which is necessary for the specification of third-whorl floral organs (stamens). Similarly, the development of fourth whorl floral organs (carpels and ovules) in angiosperms requires E-class and C-class genes [79]. Furthermore, CpSEP1/2 belonged to cluster 16 of the abortive pistils in male flowers (Fig. 4d and 4f). Therefore, we extracted a subnetwork using CpSEP1/2 as the core gene and constructed the abortive pistil regulatory network (APRN) (Fig. 5b). By analyzing the expression patterns of genes in the APRN, we found that CpARF4, CpABIG1, CpTPR3, CpIAA16, and CpGASA1 had expression patterns similar to those of the core genes SEP1/2 (Fig. S11c, Table S14). Subsequently, we performed an enrichment analysis of the APRN genes. In the cell component (CC) category, these genes were significantly enriched in GO terms such as cell–cell adherens junction (GO:0005913), apical junction complex (GO:0043296), and adherens junction (GO:0005912) (Fig. S11d, Table S18). This suggests that APRN plays a role in the regulation of cell polarity. In the biological processes (BP) category, these genes were significantly enriched in the single-organism processes (25/29, 86.21%, p < 0.05), responses to abscisic acid (8/29, 27.59%, p < 0.05), and photoinhibition GO terms (Table S18). “Single-organism process” refers to all physiological and developmental activities that occur within an individual organism that do not involve interactions with other organisms. The majority of APRN genes were enriched in this GO term, suggesting that they may be closely related to pistil abortion. However, the specific biological functions of the GO terms remain unclear. Additionally, the significant enrichment of genes responding to abscisic acid indicated that abscisic acid may be involved in pistil abortion.
In male stamens, GRN included 35,862 potential regulatory relationships involving 6,718 genes (Table S15). By analyzing gene expression levels, functional annotations, and previous clustering analysis results, we selected a subnetwork with CpTDF1 as the core gene to construct the stamen regulatory network (SRN) for male flower stamen development (Fig. 5c). In Arabidopsis, DEFECTIVE IN MERISTEM DEVELOPMENT AND FUNCTION 1 (TDF1) primarily regulates the differentiation and function of tapetal cells [80]. In papaya stamens, CpTDF1 is expressed almost exclusively during the early stages of stamen development (M1S-M3S), peaking at the M3S stage, and then remaining almost completely unexpressed in the later stages. The specific expression pattern of CpTDF1 highlights its important role during the early stages of stamen development. Additionally, in the CpTDF1 subnetwork, many genes exhibit similar expression patterns, and their homologs in Arabidopsis have been reported to be involved in pollen formation and maturation. (Fig. S11e, Table S16). For example, in Arabidopsis, the BHLH010 and BHLH091 transcription factors jointly regulate anther development and are essential for pollen wall formation [81]. In the papaya SRN, CpBHLH010 and CpBHLH091 exhibit high expression levels exclusively during stamen development, peaking at the M3S stage, indicating their involvement in stamen development and pollen wall formation. Within the SRN, CpAMS, a member of the bHLH family, is specifically expressed during stamen development. The AMS homolog in Arabidopsis is a key regulator of sporopollenin biosynthesis, secretion, and pollen wall formation, and plays a complex role in the early stages of pollen development [82]. Additionally, AMS is directly regulated by TDF1 and is essential for pollen development [83]. CpMS188, another transcription factor from the MYB family, is upregulated during stamen development. Its counterpart in Arabidopsis, MS188, regulates tapetum development and is vital for pollen exine synthesis [84]. Furthermore, GO enrichment analysis of genes within the SRN showed notable enrichment in stamen and anther development, particularly in the anther wall tapetum, in the biological process (BP) category (Fig. S11f, Table S19). The expression profiles and enrichment findings for the genes in the SRN suggest their potential role in regulating the identity of male flower stamens and the early development of the anther tapetum.
To further validate the expression of the ABCDE subfamily MADS-box and core genes in the papaya regulatory network, we conducted quantitative real-time PCR (qPCR). The quantitative expression data for these 12 genes were generally consistent with the trends observed in the transcriptome analysis (Fig. 6). B-class genes exhibited higher expression levels in the stamens of male flowers, highlighting their crucial role in the development and maintenance of stamen identity in papaya. The pistil-specific expression pattern of the CD class gene CpSTK indicates its crucial role in the normal development of papaya pistils. E-class genes showed broad expression in female pistils, aborted pistils, and stamens during development, consistent with their role in regulating the development of all floral organs as described in the ABCDE model. The specific high expression of CpSEP1/2 in aborted pistils suggests a possible functional specialization in papaya. Quantitative data for ABCDE gene expression reinforces the significance of the ABCDE model in understanding papaya flower development. Similarly, the quantitative expression of the core gene CpTDF1 within the SRN matched the transcriptomic findings. The targeted expression of CpTDF1 in stamens, along with the known function of its Arabidopsis counterpart, underscores the significance and evolutionary conservation of this gene during stamen development.

Expression patterns of 12 genes during the development of papaya female pistils (FP), aborted pistils in male flowers (MP), and stamens (MS). Relative quantification of 10 ABCDE subfamily MADS-box genes and CpAGL6 was calculated based on the expression level in the F1P sample, while CpTDF1 is relatively quantified using the M1S sample due to its specific expression in stamens. Using female pistils (FP), abortive pistils from male flowers (MP), stamens (MS) at five developmental stages, and leaf tissue to confirm the expression levels of 10 ABCDE subfamily MADS-box genes, CpAGL6, and CpTDF1