Transcriptomic and Metabolomic Analyses Reveal Mechanisms of Sexual Differentiation and Dimorphism in Morus macroura

TL;DR

This study explores how male and female mulberry plants differ at the molecular level, linking these differences to traits like leaf and fruit quality.

Contribution

The integration of transcriptomic and metabolomic data reveals sex-specific gene and metabolite patterns in M. macroura.

Findings

Female plants show higher expression of genes related to ethylene signaling and flavonoid biosynthesis.

Male plants exhibit elevated gibberellin-related gene activity and glycoside accumulation.

Sex-linked metabolic networks influence nutritional and functional traits in mulberry.

Abstract

Morus macroura ‘Panzhihua No. 1’ is a dual-purpose cultivar valued for its edible leaves and fruits. Derived from an ancient mulberry tree, it is a unique local germplasm resource. During asexual propagation, M. macroura exhibits sexual variation closely associated with fruit and leaf yield. To explore the molecular mechanisms of sexual dimorphism and its effects on nutritional traits, we integrated transcriptomic and metabolomic analyses of male and female inflorescences and leaves. Key sex-biased genes were enriched in plant hormone signaling, flavonoid biosynthesis, and carbohydrate metabolism pathways. Female plants had elevated expression of ethylene-responsive transcription factor 1 (ERF1), EIN3-binding F-box proteins (EBF1/2), and Chalcone synthase (CHS) genes and higher levels of bioactive flavonoids, including isoquercitrin and kaempferol derivatives. In contrast, male plants had increased expression of gibberellin 20-oxidase (GA20ox) and DELLA genes and accumulated glycosides, which are beneficial for leaf development. These findings reveal how sex-linked metabolic networks shape mulberry tissue functional profiles, providing molecular targets for breeding.

Figures

Full text

1. Introduction

Sexual differentiation is an evolutionarily honed survival strategy optimizing plant reproduction. It ensures stable species perpetuation through efficient genetic recombination and enhances adaptive fitness via flexible resource allocation and functional specialization of reproductive structures [1,2]. Plant sexual differentiation is a highly complex process involving the dynamic interplay of sex chromosomes, sex-determining genes, hormones and environmental factors [3]. Thus, plant sex regulation exhibits distinct polymorphism and variability. In Arabidopsis, gibberellin (GA) and cytokinin (CK) act as antagonistic signals regulating female germline specification. GA promotes megaspore mother cell (MMC) formation by degrading DELLA proteins, thereby repressing CK signaling. Conversely, CK inhibits MMC formation by activating type-B Arabidopsis response regulators (ARRs). The dynamic balance between GA and CK ensures proper germline cell development [4]. In Populus deltoides, Y-chromosome-specific genes FERR-R and MSL are involved in sex determination. FERR-R is a duplicated segment of the autosomal female-promoting gene FERR. The siRNAs it produces suppress FERR via epigenetic regulation and transcriptional cleavage, thereby repressing female traits. MSL, a member of the LTR/Gypsy transposon family, generates lncRNAs that promote stamen development. These two genes synergistically maintain the male phenotype [5]. In maize, the transcription factors Grassy tiller1 (GT1) and ZmHB13 have been confirmed to regulate unisexual flower differentiation in ears and tassels by activating jasmonate (JA) biosynthesis genes ZmLOX3 and ZmOPR7, as well as relieving the inhibition of ZmJAZ repressors. Double mutants exhibit tassel feminization and additional kernel formation in ears, elucidating the full regulatory pathway involving hormones, transcription factors, and sexual phenotypes [6].

In mulberry (Morus alba L.), sexual differentiation is regulated by multiple interacting factors, with the roles of several key regulators clarified. For instance, elevated temperatures and CO_2_ concentrations significantly increase M. alba biomass and quantity of female flowers, whereas these conditions suppress male flower development [7]. Ethrel significantly increases the proportion of M. alba female spikes, while gibberellins promote male spike development [8]. Ethrel and silver nitrate treatment of male M. alba cuttings induced female inflorescence formation [9]. At the molecular level, differentially expressed genes (DEGs) between male and female mulberry flowers are mainly enriched in pathways including starch-sucrose metabolism and plant hormone signal transduction. Key transcription factors (MADS-box, MYB), ethylene- and auxin-related genes, and sex differentiation-related genes (PMADS1/2, AGAMOUS) were identified [10]. In Morus indica, DEGs between male and female flowers are mainly enriched in the MAPK signaling pathway, flavonoid biosynthesis pathway, and galactose metabolism pathway [11]. Significant differences in the expression of 1-aminocyclopropane-1-carboxylate synthase (ACS) gene in the ethylene biosynthesis pathway were observed across six developmental stages of female mulberry flower buds [12]. Additionally, studies on Morus atropurpurea have identified the ethylene-responsive transcription factor ERF as a key candidate gene associated with mulberry flowering time [13]. Studies on members of the mulberry β-galactosidase (βGAL) family have identified several hormone-responsive cis-acting elements. βGAL gene exhibit distinct expression patterns between male and female mulberry flowers, implying their potential involvement in floral bud sex differentiation [14]. To date, the regulatory mechanism of mulberry sex differentiation remains unclear, and the key regulatory genes have not been identified, requiring further investigation.

Mulberry leaves and fruits have substantial economic value. Mulberry fruits are rich in nutrients and bioactive compounds: per 100 g fresh fruit, they contain 87.68 g water, 1.44 g protein, 0.39 g total lipid, 9.8 g carbohydrate, 1.7 g total dietary fiber, and 8.1 g total sugar; minerals include 194 mg potassium, 39 mg calcium, and 38 mg phosphorus; vitamins include 36.4 mg vitamin C, 0.87 mg vitamin E, and multiple B vitamins [15]; bioactive constituents include total phenols, flavonoids, anthocyanins, cyanidin-3-glucoside, rutin, quercetin, catechin, chlorogenic acid, gallic acid, plus unique components like polysaccharides, melatonin, β-carotene, pyrrole alkaloids, and moracin C [16,17,18,19,20,21]. These components endow mulberry fruits with immunomodulatory, antitumor, hypoglycemic, and hypolipidemic effects [22]. M. alba leaves are the core resource of the sericulture industry and also serve an important role in medicine, food, feed and cosmetics [23]. Mulberry sex differentiation exerts a significant impact on leaf nutritional value and yield. Female mulberry leaves have significantly higher nitrogen (N), crude protein, and crude ash contents than male counterparts, as female plants accumulate more nitrogen-containing resources to support subsequent fruit development. In contrast, male mulberry leaves exhibit significantly higher carbon-to-nitrogen ratios (C/N), leaf area per unit, leaf thickness, and leaf mass per unit area than female leaves, reflecting their focus on vegetative growth [8,24,25]. Studies have also shown that Bombyx mori exhibits a male-biased feeding preference, and cocoon and silk quality is closely linked to sex-specific molecular interactions between silkworms and M. alba [25]. Thus, ascertaining M. alba sex carries substantial significance for practical production.

As a key member of the mulberry germplasm bank, Morus macroura was officially designated a “Newly Discovered Excellent Germplasm Resource of Sichuan Province” in 2022. ‘Panzhihua No. 1’ is a female M. macroura cultivar developed via artificial domestication and selective breeding. Its maternal plant, located in Miyi County, Panzhihua City, Sichuan Province, is the only ancient M. macroura tree in the region, aged over 650 years. ‘Panzhihua No. 1’ has strong ecological adaptability. Its tender leaves are highly palatable to silkworms, and it produces abundant high-quality fruits with a distinct milky aroma. Additionally, its fruits have high soluble solid and sugar contents, highlighting potential for dual fruit and leaf utilization [26]. Sexual variation occurred during grafting propagation and domestication, resulting in male individuals and a stable male-to-female ratio of 1:1. This study aimed to elucidate the molecular mechanisms underlying sexual differentiation in M. macroura by integrating transcriptomic and metabolomic analyses. Specifically, we aimed to: (1) characterize phenotypic differences between male and female inflorescences and leaves; (2) identify sex-biased genes and metabolites; and (3) propose candidate regulatory pathways and genes involved in sexual differentiation. Our findings provide a valuable reference for sex-specific molecular breeding in M. macroura and related M. alba species.

2. Materials and Methods

2.1. Experimental Materials

Grafted plants of M. macroura ‘Panzhihua No. 1’, cultivated at the Shu Silkworm Culture Center in Chengdu, Sichuan Province, China, were used as experimental materials. Three biological replicates were collected for each tissue type. Each replicate consisted of pooled tissues from three randomly selected female or male plants. Female inflorescences (FI), female leaves (FL), male inflorescences (MI), male leaves (ML), female stems (FS), female roots (FR), male stems (MS), and male roots (MR) were collected on April 18. Female fruits (FF) were collected on May 13, once the fruits reached maturity. Upon collection, all samples were snap-frozen in liquid nitrogen and stored at −80 °C.

2.2. Metabolite Extraction

Samples were thawed on ice, and metabolites were extracted using a 50% methanol buffer. Each 20 μL sample was mixed with 120 μL of pre-cooled 50% methanol, vortexed for 1 min, and incubated at room temperature for 10 min. The extraction mixture was stored at −20 °C overnight. After centrifugation at 4000× g for 20 min, the supernatants were transferred to fresh 96-well plates. The samples were stored at −80 °C prior to LC-MS analysis.

2.3. Metabolite Analysis by UPLC-MS/MS

All samples were analyzed via the LC-MS system according to standard instrument protocols. First, all chromatographic separations were carried out on a Thermo Scientific UltiMate 3000 HPLC (Thermo Fisher Scientific, Waltham, MA, USA). Reversed-phase separation was achieved using an ACQUITY UPLC BEH C18 column (100 mm × 2.1 mm, 1.8 μm, Waters Corp., Milford, MA, USA). The column temperature was maintained at 35 °C, with a flow rate of 0.4 mL/min. The mobile phase consisted of water with 0.1% formic acid and acetonitrile containing 0.1% formic acid. Each sample was injected at a volume of 4 µL.

To ensure the reliability and reproducibility of metabolomic data, pooled quality control (QC) samples were prepared by mixing 10 μL of supernatant from each extracted sample. These QC samples were inserted into the sample sequence, with one QC sample per 10 experimental samples, to monitor UPLC-MS/MS system stability and metabolite extraction consistency. Key QC parameters were evaluated. The relative standard deviation (RSD) of peak intensities across all QC samples was calculated using metaX software (v1.4.8, LC-Bio Technologies, Hangzhou, China), and only metabolic features with an RSD < 30% (a widely accepted threshold for untargeted plant metabolomics) were retained for subsequent annotation and statistical analysis. Mass accuracy was calibrated every 20 samples to maintain m/z deviation within ±10 ppm. Retention time drift was corrected using XCMS software (v3.12, Scripps Research, La Jolla, CA, USA) to ensure peak alignment consistency across all samples.

Raw mass spectrometry (MS) data were converted to mzXML format via ProteoWizard’s MSConvert tool (v3.0). XCMS software was utilized for peak extraction and quality control, after which adduct ions were annotated through CAMERA (v1.52). MetaX software facilitated metabolite identification according to the Metabolomics Standards Initiative (MSI) guidelines. Primary MS data were matched against public databases for initial annotation, while secondary MS data were aligned with an internal standard compound database for MS/MS validation: primary MS data were matched against public databases for initial annotation, while secondary MS data were aligned with an internal standard compound database for validation. Candidate metabolites detected in M. macroura were annotated using Human Metabolome Database (HMDB) and Kyoto Encyclopedia of Genes and Genomes (KEGG) to clarify their physicochemical characteristics and biological functions. MetaX was subsequently used to quantify metabolites and identify differential ones. Univariate analyses, including fold-change calculation and Student’s t-test, were conducted to assess differential abundance. Benjamini–Hochberg (BH) correction was applied to calculate q-values. Differentially expressed metabolites (DEMs) were defined as those meeting the following thresholds: fold-change ≥ 1.5 or ≤0.67, and FDR-adjusted p-value (q-value) ≤ 0.05.

2.4. Transcriptomic Analysis

Total RNA from samples was isolated using Trizol reagent (Thermo Fisher Scientific, Waltham, MA, USA, cat. 15596018). RNA integrity was evaluated with the Agilent Bioanalyzer 2100 system and RNA 6000 Nano LabChip Kit (Agilent Technologies, Santa Clara, CA, USA, cat. 5067-1511). Only high-quality RNA samples were selected for the library. mRNA was purified using Dynabeads Oligo(dT) (Thermo Fisher Scientific, Waltham, MA, USA), fragmented, and reverse transcribed into cDNA using SuperScript™ II Reverse Transcriptase (Invitrogen, Thermo Fisher Scientific, Waltham, MA, USA, cat. 1896649). Double-stranded cDNA was synthesized, ligated to dual-index adapters, and amplified via PCR. The cDNA libraries were sequenced on the Illumina NovaSeq 6000 platform (LC-Bio Technology Co., Ltd., Hangzhou, China) using 150 bp paired-end reads.

2.5. Gene Functional Annotation and Expression Level Analysis

Raw sequencing reads were filtered to remove adapters, poly-A/G sequences, low-quality bases (Q < 20), and reads with more than 5% unknown nucleotides using Cutadapt (v1.9). Data quality was verified via FastQC (v0.11.9). Clean reads from all samples were aligned to the Morus reference genome (available at NCBI) using HISAT2 (v2.2.1, https://www.ncbi.nlm.nih.gov/genome/?term=morus+notabilis (accessed on 10 May 2024)). The detailed sequencing and mapping statistics for each library are provided in Supplementary Table S1. Transcript assembly and quantification were conducted using StringTie (v2.1.6), and gene expression levels were calculated as FPKM (fragments per kilobase of transcript per million mapped). Gene annotation was performed using BLASTX (v2.16.0) searches against the NCBI NR (Non-Redundant) and COG (Clusters of Orthologous Groups) databases. Functional pathways were annotated using the KEGG Automatic Annotation Server (KAAS). Differential expression analysis between female and male samples (FI vs. MI; FL vs. ML) was performed using DESeq2 software (v1.42.0) and by edgeR, with thresholds set as |log2FC| ≥ 1 and FDR-adjusted p-value < 0.05. Genes that meet the above criteria are selected as differentially expressed genes (DEGs).

2.6. Weighted Gene Co-Expression Network Analysis (WGCNA)

WGCNA was conducted using the WGCNA package (v1.71) in R (v4.1.3) [27] to investigate co-expression networks associated with sexual differentiation. Twelve samples (FI, FL, MI, and ML, with three biological replicates each) were included in the analysis. A soft-thresholding power of β = 20 was chosen to achieve a scale-free topology network (scale-free R^2^ > 0.85). Modules with at least 30 genes were identified using the dynamic tree-cutting algorithm. Highly similar modules were merged with a mergeCutHeight of 0.25, resulting in 13 distinct modules. Module-tissue correlations were calculated to identify modules associated with each tissue type. KEGG enrichment analysis of tissue-associated modules was performed using clusterProfiler (v4.6.0), with a significance threshold of adjusted p < 0.05. Pearson correlation coefficients were calculated between module eigengenes and key differential metabolites to assess gene-metabolite associations.

2.7. qRT-PCR Analyses of DEGs

Total RNA was extracted from M. macroura samples at four time points using the RNAprep Pure Plant Kit (Tiangen Biotech, Beijing, China), following the manufacturer’s instructions. RNA integrity, purity, and concentration were then assessed. Expression patterns of 15 differentially expressed genes (DEGs) involved in sex differentiation, identified through RNA-seq, were analyzed. Primers for qRT-PCR were designed using PrimerQuest (https://sg.idtdna.com/Primerquest/Home/Index (accessed on 15 April 2025)) and are listed in Supplementary Table S2. First-strand cDNA was synthesized from total RNA using the FastKing RT Kit (Tiangen Biotech, Beijing, China) with gDNA removal, according to the manufacturer’s protocol. qPCR was performed under the following conditions: initial denaturation at 96 °C for 20 s; 40 cycles of 96 °C for 10 s and 60 °C for 20 s; final extension at 60 °C for 30 s. Each 10 μL qPCR reaction contained 1 μL of cDNA template, 5 μL of qPCR Master Mix, 0.2 μL each of forward and reverse primers, and 3.2 μL of nuclease-free water. The M. macroura actin gene served as the internal reference. Relative gene expression levels were calculated using the 2^−ΔΔCT^ method [28]. All qRT-PCR assays were performed in triplicate.

3. Results

3.1. Phenotypes of M. macroura

Male M. macroura plants derived from the same maternal source exhibit significant phenotypic differences from female counterparts in both inflorescence and leaf morphology. The female inflorescence (FI, Figure 1A) comprises densely clustered unisexual florets arranged in a compact, cylindrical spike. Each floret bears a prominent pistil with a bifurcated stigma, adapted for efficient pollen capture. After pollination, the stigma withers, and the ovary wall and perianth enlarge to form a fleshy syncarp. Developing fruits change color from green to red, eventually maturing into dark purple. The male inflorescence (MI, Figure 1B) is protandrous, with florets maturing earlier than those in the female inflorescence. The flowers are arranged in pendulous, catkin-like spikes composed of numerous small staminate florets. Each floret contains four basifixed anthers that undergo longitudinal dehiscence upon maturation, releasing abundant wind-dispersed pollen. During anthesis, the floral axis elongates significantly, enhancing pollen dispersal efficiency. Female leaves (FL, Figure 1C) are broadly ovate to orbicular, with entire margins and a glabrous, glossy adaxial surface. The lamina displays a cordate base and acute apex, typical characteristics of M. macroura leaves. Male leaves (ML, Figure 1D) are deeply lobed, featuring 3–5 palmate lobes. The adaxial surface is pubescent with trichomes concentrated along the veins, whereas the abaxial surface bears a denser indumentum.

3.2. Metabolomic Profiling and Differential Metabolite Analysis

Quality control (QC) samples were analyzed to assess the reliability of the metabolomic data. Principal component analysis (PCA) of all samples, including QC samples, showed tight clustering of the QCs (Supplementary Figure S1), indicating excellent stability and reproducibility of the LC-MS system throughout the analysis. Metabolic features with a relative standard deviation (RSD) > 30% in QC samples were filtered out, ensuring that only high-quality data were used in downstream analysis. Stringent quality control led to the confident detection of 25,345 metabolic features in the inflorescences and leaves of M. macroura using UPLC-MS/MS. Metabolite identification was carried out at multiple confidence levels. High-confidence identification (Level 1/Level 2) was achieved for 893 features by matching MS/MS fragmentation patterns against commercial, public spectral databases (e.g., HMDB, MassBank), and an in-house standard library. A further 14,759 features were putatively annotated at Level 2 by accurate mass matching (mass error < 10 ppm) to the HMDB and KEGG databases. The remaining 9693 features were detected but could not be annotated and were classified as unknown (Level 3).

Cluster heatmap analysis was conducted to compare overall metabolic profiles between FI and MI, and between FL and ML (Figure 2A). Cluster analysis showed high reproducibility among samples within each group. As shown in Figure 2B, metabolites in M. macroura inflorescences and leaves were mainly 272 lipids and lipid-like molecules, 133 phenylpropanoids and polyketides, and 105 organic acids and their derivatives. These findings are consistent with previously reported nutrient compositions of mulberry leaves [29,30]. Principal component analysis (PCA) was conducted to examine overall metabolic differences between male and female M. macroura inflorescences and leaves, as well as within-group variation (Figure 2C). In the FI vs. MI comparison, the first two principal components (PC1 + PC2) explained 67.34% of the variance, and 52.25% in the FL vs. ML comparison. The PCA score plots showed clear separation between groups, indicating significant metabolic differences between FI and MI and between FL and ML.

Based on fold-change (FC) ≥1.5 or ≤0.67 and p-value ≤ 0.05, 7786 DEMs were identified in the FI vs. MI comparison, including 3849 upregulated and 3937 downregulated metabolites (Figure 2D). These DEMs were mainly classified into glycerophospholipids, carboxylic acids and their derivatives, flavonoids, and benzene and substituted derivatives. Salicylic acid, 3,4-dihydroxybenzaldehyde, benzoic acid, isoquercitrin, and jasmonic acid exhibited significantly higher accumulation in FI, whereas luteolin, isorhamnetin, quercetin, 2-acyl-sn-glycero-3-phosphoserine, and quercetin-3-(6″-malonyl)-glucoside were predominantly enriched in MI. In the FL vs. ML comparison, 4763 DEMs were identified, including 2471 upregulated and 2292 downregulated metabolites. Most DEMs belonged to flavonoids, benzene and substituted derivatives, and organic oxides. Among these, salicylic acid, kaempferol-3-O-glucoside, robinin, and quercetin 3,4′-di-O-glucoside were markedly upregulated in FL, whereas quercetin-3-(6″-malonyl)-glucoside, protohypericin, L-arginine, and 1,5-dicaffeoylquinic acid displayed higher abundance in ML (Figure 2E). Detailed information on these DEMs is provided in Supplementary Table S3.

3.3. KEGG Analysis of Differential Metabolites

Pathway enrichment analysis of DEMs in both groups was conducted using the KEGG database. In FI vs. MI, 87 DEMs were enriched across 179 pathways, while in FL vs. ML, 39 DEMs were mapped to 116 pathways (Supplementary Table S4). DEMs in the FI vs. MI comparison were significantly enriched in Environmental Information Processing, Metabolism, and Organismal Systems categories. The metabolic pathway contained the highest number of enriched DEMs, totaling 61. Other enriched pathways included protein digestion and absorption, plant secondary metabolite biosynthesis, ABC transporters, glycerophospholipid and glycerolipid metabolism, flavonoid biosynthesis, and plant hormone signal transduction. In FL vs. ML, DEMs were primarily enriched in Metabolism, Environmental Information Processing, and Cellular Processes. Enriched pathways included glycerophospholipid metabolism, flavonoid biosynthesis and degradation, plant hormone signal transduction, and starch and sucrose metabolism. The metabolic pathway remained the most enriched, with 29 DEMs, followed by plant secondary metabolite biosynthesis (22 DEMs) and flavonoid biosynthesis (7 DEMs) (Figure 2F).

3.4. Transcriptomic Analysis and Identification of DEGs

To ensure the reliability of downstream differential expression analysis, we evaluated the quality of the RNA-seq data. A total of 12 cDNA libraries were sequenced, generating an average of 39.38 million clean reads per sample after filtering, with Q30 percentages ranging from 96.73% to 97.71%. The clean reads were aligned to the Morus notabilis reference genome, achieving an average mapping rate of 73.1% (ranging from 69.37% to 75.14%), with an average unique mapping rate of 70.5% (Supplementary Table S1). These metrics indicate high-quality sequencing data suitable for subsequent differential expression analysis. Differentially expressed genes (DEGs) were identified in the comparison groups using FDR-adjusted p-value < 0.05 and |log2FC| ≥ 1 as threshold criteria. In FI vs. MI, 5042 DEGs were identified, including 1696 upregulated in FI and 3346 downregulated. In FL vs. ML, 2513 DEGs were identified, with 1142 upregulated in FL and 1371 downregulated (Figure 3A). A Venn diagram (Figure 3B) revealed 1196 co-expressed DEGs shared between the two comparison groups. Detailed information on these DEGs is available in Supplementary Table S5. These findings suggest that DEGs identified in both FI vs. MI and FL vs. ML play key regulatory roles in sexual dimorphism and sex determination in M. macroura.

3.5. Functional Enrichment Analysis of DEGs

To elucidate the biological functions of DEGs in M. macroura, GO enrichment analysis was conducted for functional categorization across the comparison groups. A total of 4737 and 2320 DEGs were identified in the FI vs. MI and FL vs. ML comparisons, respectively. These DEGs were assigned to 3062 and 2139 GO terms, respectively. Detailed information is provided in Supplementary Table S6.

In the FI vs. MI comparison, 99 GO terms showed significant enrichment, with the top 50 visualized in Figure 3C. Enrichment was primarily observed in the biological process (BP; 32 terms) and molecular function (MF; 13 terms) categories. The biological process subcategory contained 1653 functional subclasses, comprising 9919 annotations. Among them, defense response (211), flavonoid biosynthetic process (85), response to wounding (68), and signal transductionsignal transduction (105) subcategories were most significantly enriched. The MF subcategory included 1061 functional subclasses and 8118 annotations. Significantly enriched subclasses included protein serine/threonine kinase activity (256), naringenin-chalcone synthase activity (14), and DNA-binding transcription factor activity (335). The cellular component (CC) subcategory comprised 348 functional subclasses with 7864 annotations. The most enriched functional subclasses in the CC category were the extracellular region (341), plasma membrane (751), and membrane (705).

In the FL vs. ML comparison, 79 GO terms were significantly enriched, with the top 50 terms displayed in Figure 3D. Enrichment was primarily observed in biological process (BP; 30 terms) and molecular function (MF; 16 terms) categories. The biological process subcategory included 1181 functional subclasses with 4937 annotations. Significantly enriched subclasses included defense response (151), protein phosphorylation (148), response to wounding (52), and signal transduction (73). The MF subcategory included 1061 functional subclasses, comprising 8118 annotations. Significantly enriched subclasses included protein serine/threonine kinase activity (174), ADP binding (67), and kinase activity (116). The CC subcategory comprised 218 functional subclasses with 3799 annotations. The most enriched functional subclasses in the CC category were plasma membrane (433), plasmodesma (121), and the extrinsic component of the plasma membrane (12).

To further elucidate the functional roles of DEGs between floral and foliar tissues, KEGG enrichment analysis was performed. In the FI vs. MI, 1669 DEGs were assigned to 130 KEGG pathways (detailed in Supplementary Table S7), with the top 20 significantly enriched pathways visualized in Figure 3E. As shown in Figure 3F, these pathways were predominantly categorized into Cellular Processes, Organismal Systems, Metabolism, and Environmental Information Processing; most of them fell into the metabolism category. Key enriched pathways included flavonoid biosynthesis, circadian rhythm-plant, MAPK signaling pathway-plant, plant hormone signal transduction, glycolysis/gluconeogenesis, and starch and sucrose metabolism. The pathways with a high number of DEGs were plant-pathogen interaction (116), plant hormone signal transduction (94), MAPK signaling pathway-plant (75), and sucrose metabolism (73). In the FL vs. ML, 868 DEGs were mapped to 128 KEGG pathways (detailed in Supplementary Table S7), with the top 20 enriched pathways visualized in Figure 3G. The functional categorization of these pathways remained consistent with that of FI vs. MI, and the majority still fell into metabolic categories, as depicted in Figure 3H. Enriched pathways included circadian rhythm-plant, flavonoid biosynthesis, plant hormone signal transduction, glycolysis/gluconeogenesis, and plant-pathogen interaction, among others. Pathways with the highest DEG counts were plant-pathogen interaction (69), plant hormone signal transduction (61), flavonoid biosynthesis (37), and MAPK signaling pathway-plant (34).

3.6. Tissue-Specific Co-Expression Modules Associated with Sexual Differentiation

WGCNA identified 13 distinct modules ranging in size from 280 to 4860 genes based on a soft-thresholding power of β = 20 (Supplementary Figure S2). Module-condition correlation analysis revealed modules specifically associated with each tissue type (Supplementary Figure S3): MEgreen (584 genes) with FI (r = 0.93, p = 2 × 10^−5^), MEpink (280 genes) with FL (r = 0.95, p = 2 × 10^−6^), MEblue (2444 genes) with MI (r = 0.96, p = 6 × 10^−7^), and MEturquoise (4860 genes) with ML (r = 0.78, p = 0.003).

Correlation analysis between module eigengenes and key metabolites revealed distinct metabolic associations for each tissue-specific module (Figure 4A). The FI-associated MEgreen module was positively correlated with salicylic acid (r = 0.597, p = 0.041) in FI samples. The FL-associated MEpink module showed strong positive correlations with isoquercitrin (r = 0.967, p = 3.0 × 10^−7^), trehalose (r = 0.944, p = 4.0 × 10^−6^), and cyanidin (r = 0.622, p = 0.031), and negative correlations with jasmonic acid (r = −0.846, p = 0.0005) and L-arginine (r = −0.736, p = 0.006). The MI-associated MEblue module was positively correlated with luteolin (r = 0.972, p = 1.4 × 10^−7^), quercetin (r = 0.950, p = 2.2 × 10^−6^), and jasmonic acid (r = 0.651, p = 0.022), and negatively correlated with trehalose (r = −0.762, p = 0.004). The ML-associated MEturquoise module exhibited a negative correlation with jasmonic acid (r = −0.676, p = 0.016) in ML samples.

To further explore the functional basis underlying these metabolic associations, KEGG enrichment analysis was performed on each tissue-associated module, revealing biologically relevant pathways. The FI-associated MEgreen module was significantly enriched in flavonoid biosynthesis (map00941, p = 2.44 × 10^−15^), plant hormone signal transduction (map04075, p = 6.39 × 10^−5^), and circadian rhythm (map04712, p = 1.66 × 10^−11^) (Figure 4B). The FL-associated MEpink module showed enrichment in starch and sucrose metabolism (map00500, p = 7.00 × 10^−8^) and plant hormone signal transduction (map04075, p = 0.036) (Figure 4C). The MI-associated MEblue module was significantly enriched in starch and sucrose metabolism (map00500, p = 1.86 × 10^−15^), galactose metabolism (map00052, p = 3.99 × 10^−15^), plant hormone signal transduction (map04075, p = 1.12 × 10^−17^), and flavonoid biosynthesis (map00941, p = 2.82 × 10^−8^) (Figure 4D). The ML-associated MEturquoise module exhibited strong enrichment in photosynthesis (map00195, p = 3.69 × 10^−28^), carbon fixation (map00710, p = 5.06 × 10^−15^), and starch and sucrose metabolism (map00500, p = 3.49 × 10^−25^) (Figure 4E). These pathway enrichments are consistent with the observed metabolite correlations, further validating the functional relevance of these modules in sex- and tissue-specific metabolic regulation.

3.7. Validation of DEGs by qRT-PCR

To verify the transcriptome sequencing results, 15 DEGs (Supplementary Figure S4) were selected from four pathways for qRT-PCR validation, including plant hormone signal transduction (EBF1, EBF2, ERF1_1, ERF1_2, DELLA, GA2ox, MYC2), flavonoid biosynthesis (CHI, CHS1, CHS2, DFR), starch and sucrose metabolism (GLC1, GLC2), and galactose metabolism (βGAL, αGAL). The results showed significant expression differences between male and female plants (p < 0.05), consistent with RNA sequencing trends. This confirmed the reliability of Illumina sequencing data. Furthermore, key gene expression levels across tissues were determined to explore their potential roles in sex differentiation, with detailed analyses in the Discussion.

4. Discussion

Dioecy is rare in angiosperms, accounting for only 6% of species [31]. In contrast, M. alba exhibits remarkable sexual plasticity, with diverse sex phenotypes including female, male, monoecious, and male-biased individuals [32]. Sex expression in M. alba is co-regulated by multiple factors. Numerous studies have confirmed that environmental stressors, along with increasing ambient concentrations of endocrine-disrupting chemicals, increasingly trigger sex reversal in M. alba, challenging its inherently stable sex phenotypes [10,33]. In M. macroura, sexual differentiation not only determines plant reproductive fate by forming distinct sex phenotypes but also affects the nutritional value of its inflorescences and leaves, making it a valuable dual-purpose cultivar. Transcriptomic and metabolomic analyses of two comparison groups, inflorescences (FI vs. MI) and leaves (FL vs. ML), revealed overlapping enrichment in pathways of plant hormone signal transduction, flavonoid biosynthesis, carbohydrate metabolism, circadian rhythm, and plant-pathogen interaction. These results indicate that sexual differentiation is closely correlated with and deeply integrated into metabolic regulation. Flavonoid-mediated redox balance, hormone-regulated organogenesis, and carbon allocation to reproductive organs may act synergistically to shape sex-specific traits. Similar mechanisms have also been observed in Morus indica [11].

4.1. Plant Hormone Signal Transduction Is a Critical Pathway in Mediating Sexual Differentiation in M. macroura

4.1.1. Auxin–Ethylene Signaling Is Associated with Female Flower Development

Auxin (AUX) and ethylene (ETH) are established regulators of female floral fate in many plant species. In Cucumis sativus, AUX concentration gradients determine the developmental fate of pistil and stamen primordia [34,35]. AUX homeostasis governs gynoecium development in Carica papaya, with an optimal AUX concentration and normal polar AUX transport critical for functional gynoecium development [36]. Auxin response factor (ARF) is closely associated with floral primordium initiation and induces indole-3-acetic acid amido synthetase (GH3) gene expression [37,38,39]. In M. macroura, six ARF and four GH3 genes were downregulated in female tissues, while most auxin-responsive (IAA) genes were upregulated in FI and downregulated in FL (Figure 5A). This expression pattern is consistent with the conserved Aux/IAA-ARF regulatory module in Dendrobium officinale, which mediates organ-specific sexual development via auxin signaling [40].

ETH is known to participate in floral organ morphogenesis in M. alba and has been suggested to play a regulatory role in floral development [32]. In M. macroura, the ethylene signaling pathway showed enhanced activity in female tissues based on our transcriptomic data (Figure 5B), and qRT-PCR validation revealed sex- and tissue-specific expression of ethylene-related genes: ERF1 (ethylene-responsive transcription factor 1) and EBF1/2 (EIN3-binding F-box proteins). As negative regulators of the ethylene signaling pathway, EBF1 and EBF2 expression patterns were validated by qRT-PCR (Figure 5C): EBF1 was significantly more highly expressed in MI than in other male tissues; EBF2 was also significantly more highly expressed in MI and MS than in other tissues, with the lowest level in MR. EBF1/2-encoded proteins degrade EIN3 and attenuate ethylene signal transduction [41,42,43]; their high expression in MI may therefore be associated with limiting ethylene levels, potentially contributing to male floral development. However, in female tissues, the observed expression patterns of EBF1/2 are reminiscent of a feedback regulatory mechanism that could fine-tune ethylene response activity rather than completely suppress it. EBF1 was relatively highly expressed in FI, FS, and FF—significantly higher than in other female tissues, while EBF2 was most highly expressed in FI. This pattern parallels the EIN3-EBF negative feedback loop described in Arabidopsis, where ethylene signaling remains active despite EBF upregulation, owing to functional redundancy among EIN3-like (EIL) proteins [44]. As a downstream transcription factor of the ethylene signaling pathway, ERF1 was highly expressed in FR and reproductive tissues, suggesting that ethylene signaling plays a role in female floral development.

4.1.2. The Gibberellin–DELLA Module Is Linked to Male Trait Formation

Gibberellins (GA) antagonize ETH and promote male reproductive organ development by repressing ethylene biosynthesis and destabilizing DELLA proteins [45,46,47]. In M. macroura, most DELLA genes were upregulated in male tissues (Figure 6A), which may reflect a compensatory regulatory mechanism whereby DELLA gene transcriptional levels are elevated to counteract GA-mediated DELLA protein degradation and maintain the homeostasis of functionally active DELLA proteins. DELLA proteins have been proposed to switch function from growth repressors to activators and could potentially drive male trait development by interacting with male-specific transcription factors. GA20ox (gibberellin 20-oxidase) and GA3ox (gibberellin 3-oxidase), key genes in the GA biosynthetic pathway, were also upregulated in male inflorescences of M. macroura, which is consistent with their reported roles in stamen development and pollen maturation [48,49]. Tissue-specific expression of GA20ox and DELLA genes supports their putative roles in floral organ patterning. In male plants, GA20ox was highly expressed in roots and inflorescences, while DELLA expression showed a decreasing trend from roots to inflorescences; this reduction, likely due to GA-mediated degradation, may permit floral development. In female plants, stably expressed DELLA genes may sustain ethylene signaling and thus promote pistil development (Figure 6C).

4.1.3. Potential Role of Jasmonic Acid in Regulating Flowering and Metabolism in Female Tissues

Previous studies have shown that jasmonic acid (JA) modulates transcription factor expression via the JAZ-MYC module, participating in the regulation of floral development and nutrient metabolism [50,51]. In FI of M. macroura, three jasmonate ZIM-domain protein (JAZ) genes were downregulated and two were upregulated, while all four myelocytomatosis (MYC) genes were significantly upregulated (Figure 6B). JAZ1 is relatively highly expressed in FI and MI, with significantly higher levels in MI than other organs. MYC2 is most highly expressed in roots of both sexes, and significantly more highly expressed in FI than MI. This positively correlates with JA accumulation, suggesting active JA signaling in these tissues (Figure 6D). The observed expression patterns suggest that JA pathway activation may be associated with the coordinated regulation of floral development and the synthesis of secondary metabolites, including flavonoids and phenolic compounds [52,53,54]. In addition to AUX, ETH, GA and JA, we identified 50 ABA-related, 15 SA-related and 30 BR-related DEGs, most of which were upregulated in FI and FL (Supplementary Figure S5). These results highlight complex hormonal crosstalk in sexual dimorphism regulation, where antagonistic and synergistic hormonal actions orchestrate organ identity and development. Detailed information regarding the DEGs of plant hormone signal transduction is presented in Supplementary Table S8.

4.2. Flavonoid Biosynthesis Pathway Correlates with Sex-Specific Traits

Flavonoids are core metabolites sustaining plant reproductive fitness and nutritional quality [55,56,57]. Chalcone synthase (CHS), chalcone isomerase (CHI) and dihydroflavonol reductase (DFR), key structural genes of the flavonoid biosynthetic pathway, were significantly upregulated in female tissues, correlating with the accumulation of isoquercitrin and kaempferol derivatives in FI. These metabolites could potentially contribute to the high soluble solids content and superior functional value of FF. In contrast, MI were enriched in luteolin and quercetin derivatives, which may enhance pollen viability and improve leaf palatability for silkworms (Figure 7A), with details provided in Supplementary Table S9. Their expression patterns were consistent with those reported for Broussonetia papyrifera [51]. qRT-PCR validated the sex-specific expression patterns of CHI, CHS1 and CHS2 across different tissues. CHI, CHS1, and CHS2 were relatively highly expressed in FR and MR. CHS1 and CHS2 showed low expression in leaves and stems of both sexes, with higher levels in MI than in FI. In contrast, CHI was more highly expressed in FI than in MI (Figure 7B). These results suggest that CHS1 and CHS2 may mediate the synthesis of male-specific flavonoids, while CHI underpins female reproductive function.

4.3. Carbohydrate Metabolism Correlates with Sex-Specific Organ Development

Carbohydrate metabolism is a pivotal pathway providing energy and precursor metabolites for hormone biosynthesis, organ development and nutrient storage [58,59]. In M. macroura, 15 glucoendo-1,3-beta-glucosidase (GLC) genes were differentially expressed, with 4 upregulated in FI and 8 in FL (Figure 8A). These genes could be involved in starch degradation and sucrose accumulation in female tissues. In the galactose metabolic pathway, 11 β-galactosidase (βGAL) genes showed contrasting expression patterns, with 3 upregulated in FI and 11 in ML. This result was consistent with reports on Camellia sinensis [60]. βGALs expression patterns in M. macroura suggest potential roles in FF development or cell wall structural remodeling in ML. Four α-galactosidase (αGAL) genes also showed sex-specific expression, with 1 upregulated in female plants and 3 in male plants (Figure 8B). This expression pattern is similar to that reported for αGAL genes in Vernicia fordii [61]. Their expression in male plants may contribute to raffinose oligosaccharide breakdown and help maintain plant sugar homeostasis. Tissue-specific expression of GLC1/GLC2 and βGAL1/βGAL2 provides further insights into their potential biological functions: GLC1 was most highly expressed in roots of both female and male M. macroura plants, showing similar trends between sexes and generally higher expression in males. GLC2 expression in females increased first and then decreased from FI to FR, peaking in FS; in males, it increased progressively from MI to MR, with significantly lower levels in MI and ML than MS and MR (Figure 8C). βGAL1 showed identical expression patterns in both sexes, with peak expression in inflorescences (FI, MI) and significantly higher expression in MI than FI. βGAL2 was significantly more highly expressed in FI and FS compared with MI and MS. αGAL1 expression was comparable across organs of both sexes, with slightly higher expression in MI and MR than in FI and FR (Figure 8D). Taken together, these expression patterns suggest that these genes may be involved in sugar allocation and could contribute to organ-specific physiological functions across plant organs.

4.4. Glycolysis/Gluconeogenesis Is a Central Metabolic Hub

The glycolysis and gluconeogenesis pathway is known to be pivotal for adenosine triphosphate (ATP) production and phytohormone precursor synthesis [62,63]; its intermediate metabolites, including phosphoenolpyruvate (PEP) and pyruvate, have been shown to directly participate in GA and AUX biosynthesis [39,64]. A total of 94 DEGs were identified in this pathway, and key enzyme-encoding genes including hexokinase (HK), phosphoglucomutase (PGM) and glucose-6-phosphate isomerase (GPI) all exhibited differential expression patterns: PGM and GPI genes were downregulated in FI and upregulated in FL, whereas the expression trends of HK genes varied across different tissues (Figure 8E), with details provided in Supplementary Table S10. These core enzymes may participate in linking glycolysis with starch-sucrose and galactose metabolism via shared intermediate metabolites, including glucose-6-phosphate (G6P), glucose-1-phosphate (G1P), and fructose-6-phosphate (F6P). HK initiates the glycolytic cascade by catalyzing glucose conversion to G6P [65,66]; phosphoglucomutase mediates G6P involvement in starch synthesis or degradation [67]; in galactose metabolism, the Leloir pathway converts galactose to G1P via GALK and GALT, and channels G1P into central carbon metabolism [68,69]. This coordinated metabolic network allows plants to flexibly regulate carbon allocation and balance energy production and reproductive development demands. Although research on glycolysis and gluconeogenesis pathway involvement in plant sexual differentiation remains limited, our results suggest that this pathway may form a key metabolic foundation for hormonal regulation and organ-specific growth in M. macroura.

Based on our integrated transcriptomic and metabolomic analyses and previous reports on gene functions, we have summarized the key genes associated with sexual differentiation in M. macroura (Figure 9). The schematic integrates the sex-biased genes identified in this study with their associated metabolic pathways. In female plants, genes involved in auxin signaling (IAA), ethylene signaling (EBF1/2, ERF1), and flavonoid biosynthesis (CHI, DFR, CHS) were predominantly upregulated. In male plants, upregulated genes included those related to auxin response (ARF, GH3), gibberellin signaling and biosynthesis (DELLA, GA20ox, GA3ox), jasmonate signaling (JAZ), and carbohydrate metabolism (βGAL, αGAL, GLC, PGM, GPI). The interaction among these pathways, particularly the crosstalk between phytohormone signaling and metabolic processes, appears to contribute to the sexual dimorphism observed in M. macroura.

5. Conclusions

Sexual differentiation in M. macroura is a complex biological process. This study investigated the regulatory mechanisms governing sexual differentiation in M. macroura. Integrated transcriptomic and metabolomic analyses identified that pathways such as phytohormone signaling, flavonoid biosynthesis, and carbohydrate metabolism are potentially involved in sexual regulation, with dynamic interactions likely coordinating their functions. Female M. macroura plants exhibited increased ethylene signaling (ERF1, EBF1/2) and flavonoid-related gene expression (CHS, CHI, DFR), which were correlated with stress-adaptive metabolites such as isoquercitrin and salicylic acid. Conversely, male plants showed upregulated expression of genes involved in gibberellin biosynthesis (GA20ox, GA3ox) and carbohydrate metabolism (αGALs, βGALs, GLCs, PGMs, GPIs, HKs), suggesting their potential involvement in the accumulation of glycolipids and kaempferol derivatives that may support pollen development. These findings reveal a potential ternary regulatory network involving plant hormone signal transduction, flavonoid biosynthesis, and carbohydrate metabolism. Furthermore, the large number of metabolites and differentially expressed genes identified in this study serve as a valuable reservoir of candidate genes and pathways for future functional investigations. Their definitive roles in M. macroura sexual differentiation await confirmation through targeted experimental approaches such as genetic transformation or gene editing. This work not only provides critical insights but also offers novel avenues for breeding strategies in this and related dioecious plant species.