The Mechanism of the Development and Maintenance of Sexual Dimorphism in the Dioecious Mulberry Plant (Morus alba)

Citation: Shi, Y.; Ackah, M.; Amoako, F.K.; Zhao, M.; van der Puije, G.C.; Zhao, W. The Mechanism of the Development and Maintenance of Sexual Dimorphism in the Dioecious Mulberry Plant (Morus alba). Biology 2024,13, 622. https://doi.org/ 10.3390/biology13080622 Academic Editor: Wenqiang Li Received: 15 June 2024 Revised: 8 August 2024 Accepted: 13 August 2024 Published: 15 August 2024 Copyright: © 2024 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/). biology Article The Mechanism of the Development and Maintenance of Sexual Dimorphism in the Dioecious Mulberry Plant (Morus alba) Yisu Shi 1,2, Michael Ackah 1,2 , Frank Kwarteng Amoako 3, Mengdi Zhao 4,*, Grace C. van der Puije 5 and Weiguo Zhao 1,2,* 1Jiangsu Key Laboratory of Sericulture Biology and Biotechnology, School of Biotechnology, Jiangsu University of Science and Technology, Zhenjiang 212100, China; [email protected] (Y.S.); [email protected] (M.A.) 2 Key Laboratory of Silkworm and Mulberry Genetic Improvement, Ministry of Agriculture and Rural Affairs, The Sericultural Research Institute, Chinese Academy of Agricultural Sciences, Zhenjiang 212100, China 3Institute of Plant Nutrition and Soil Science, Kiel University, Hermann-Rodewald-Straße 2, 24118 Kiel, Germany; [email protected] 4Department of Materials Science and Engineering, Suzhou University of Science and Technology, Suzhou 215011, China 5Department of Crop Science, School of Agriculture, College of Agriculture and Natural Sciences, University of Cape Coast, Cape Coast P.O. Box 5007, Ghana; [email protected] *Correspondence: [email protected] (M.Z.); [email protected] (W.Z.) Simple Summary: This study highlights the importance of intersexual differentiation in dioecious plant speciation and maintenance. Using Morus alba plants, we identified 1543 differentially expressed genes (DEGs) between male and female individuals. Notably, 333 genes were exclusively expressed in male inflorescences, while 66 were unique to female inflorescences. Through comparative transcriptomics, 841 genes with significantly higher expression in males were linked to defense-related pathways, and 702 genes with female-biased expression were related to growth and development pathways. Signals of selection were detected in both maleand female-biased genes, indicating their role in evolution. These findings advance our understanding of the molecular basis of intersexual differentiation and shed light on sex evolution in M. alba plants. Abstract: Intersexual differentiation is crucial for the speciation and maintenance of dioecious plants, but the underlying mechanisms, including the genes involved, are still poorly understood. Here, we focused on a typical dioicous plant Morus alba, to explore the molecular footprints relevant to sex evolution by revealing the differentially expressed genes (DEGs) between two sexes and the testing signals of selection for these DEGs. From the results, we found a total of 1543 DEGs. Interestingly, 333 and 66 genes expression were detected only in male and female inflorescences, respectively. Using comparative transcriptomics, the expression of 841 genes were found to be significantly higher in male than in female inflorescences and were mainly enriched in defense-related pathways including the biosynthesis of phenylpropanoids, cutin, suberine and waxes. Meanwhile, the expression of 702 genes was female-biased and largely enriched in pathways related to growth and development, such as carbohydrate metabolism, auxin signaling and cellular responses. In addition, 16.7% and 17.6% signals of selection were significantly detected in femaleand male-biased genes, respectively, suggesting their non-negligible role in evolution. Our findings expanded the understanding of the molecular basis of intersexual differentiation and contribute to further research on sex evolution in dioecious plants. Keywords: sex-biased gene expression; sexual dimorphism; sexual selection; pollen competition; haploid selection; artificial selection; species diversity Biology 2024,13, 622. https://doi.org/10.3390/biology13080622 https://www.mdpi.com/journal/biology Biology 2024,13, 622 2 of 17 1. Introduction Ongoing global change poses a serious threat to the conservation of plant biodiversity and ecosystem functioning [ 1 ]. For conservation purposes, it is therefore essential to elucidate how plants respond to changing environments, which requires a comprehensive understanding of their life history [ 2 ]. Dioecious plants have many different reproductive traits that have evolved under different selective forces, and the fitness of individuals with different sexes can be influenced by different factors [ 3 – 6 ]. However, although such intersexual differences call for sex-specific conservation strategies within a single dioecious plant species [ 3 , 7 ], the details of trait differences and the underlying molecular basis and evolutionary trajectories are still poorly understood. Therefore, studies on the gene expression patterns relevant to the traits constituting sexual dimorphism are of great importance for dioecious plants. Some traits without a direct relation to gamete development but show significant sexual dimorphism are called secondary sexual characteristics [ 7 ]. These traits involved in sexual dimorphism in dioecious species mainly include flowering phenology, floral characters (structure, size, and number of flowers in an inflorescence), floral chemical composition such as nutrient content and herbivore-resistant chemicals, and physical defenses [ 7 – 9 ]. Furthermore, these traits are thought to influence fitness by affecting pollination, mating, fecundity, and seed viability [ 10 ]. Previous studies have shown that in most dioecious plants, males tend to produce more flowers with shorter flower longevity and less synchronous flowering compared to females, and their evolution is likely driven by competition for pollination efficiency and mating competition [ 11 ]. In addition, there may be stronger chemical and physical defenses of male flowers in wind-pollinated plants than insect-pollinated species, which may help to ensure sufficient pollen for pollen dispersal [ 11 – 13 ]. As for females, their flowers tend to synthesize more auxins to initiate the reproductive program and stimulate seed development, as well as nutrients that can both attract pollinators and facilitate seed growth [ 14 , 15 ]. However, the molecular basis of these differentiated reproductive traits is still unclear. Although some pathways involved in plant metabolism have been identified [ 7 , 16 , 17 ]. Fortunately, the rapid development of transcriptomics and comparative genomics has shed light on unravelling the molecular mechanisms underlying trait variation by detecting variable expression and signals of positive selection of candidate genes [5,18]. Species in the genus Morus are widely distributed and provide an important pool of species for agriculture, but many wild species are threatened [ 19 ]. Mulberry (M. alba) (2n = 28) is a dioecious wind-pollinated perennial woody plant, cultivated in Asia [ 20 , 21 ]. M. alba exhibits apparent sexual dimorphism in flower structure, number and size per inflorescence and flower earlier in males than in females [ 21 – 23 ]. This provides an excellent model for screening genes associated with sexual dimorphism. In addition, a high-quality assembly with a genome size of 346.39 Mb has been reported for this species [ 21 ], providing a profound basis for transcriptome studies. Furthermore, there are many dioecious and monoecious species in the genus Morus [ 24 ], so the presence of significant signals of positive selection in the screened genes can be easily tested by interspecific comparisons with closely related monoecious species. In this study, we used comparative transcriptome analysis to screen differentially expressed genes (DEGs) between two sexes of M. alba inflorescences collected from three male and three female trees and performed functional enrichment on selected DEGs relevant to traits of sexual dimorphism, with further evolutionary analysis to uncover their evolutionary history. We aim to answer two scientific questions: (1) whether DEGs in male and female inflorescences are involved in regulating the development of sex-specific traits, and (2) whether these genes are under positive selection. Our results will provide theoretical evidence and elucidate the molecular basis of the origin and maintenance of sexual dimorphism in plants and will aid future conservation efforts. Biology 2024,13, 622 3 of 17 2. Materials and Methods 2.1. Plant Materials The flower buds were collected from the trees of two M. alba purified breeding lines, three BaiTiao (BT), which produces only male flowers, and three Y20 which produces only female flowers [ 25 ]. These breeds are grown at the Sericulture Research Institute of the Chinese Academy of Agricultural Sciences in Zhenjiang, Jiangsu Province, China (N 32 ◦ 11 ′ 45.80 ′′ , E 119 ◦ 23 ′ 45.80 ′′ ). Maintenance conditions for the M. alba tree are as follows: minimum and maximum temperatures (9–20 °C), relative humidity (40–60%), irrigation (rainfall), day length (12 h), and light intensity (15,000 lx). Additionally, information about the appearance conditions of the trees includes age (5 year), diameter (5 cm), and height (2 m). BT and Y20 have characteristics such as larger inflorescences, a greater number of flowers, and wide adaptability and are typical cultivars of M. alba. Therefore, these two purified lines were used as the study materials to compare the differences in gene expression between male and female inflorescences of M. alba. Male and female inflorescence (1.3 m from the ground) samples were collected (that is, male inflorescence from BT purified lines and female inflorescence from the purified Y20 lines). The samples were a mixture of different inflorescences from the same plant, and six trees (3 from BT and Y20 each) were selected for sampling, representing three biological replicates. In all, a total of six catkins were sampled. Inflorescence sampling was conducted from 9:00 p.m. to 11:00 p.m. Only fully developed adult catkins (0.8–1.5 cm long and 0.3–0.5 cm in diameter) were collected, when the catkins were fully expanded, but before the flowers had opened (Figure 1). Whole catkins were cut and immediately frozen in liquid nitrogen for RNA extraction. Biology2024,13,xFORPEERREVIEW3of17   2.MaterialsandMethods 2.1.PlantMaterials TheflowerbudswerecollectedfromthetreesoftwoM.albapurifiedbreedinglines, threeBaiTiao(BT),whichproducesonlymaleflowers,andthreeY20whichproducesonly femaleflowers[25].ThesebreedsaregrownattheSericultureResearchInstituteofthe ChineseAcademyofAgriculturalSciencesinZhenjiang,JiangsuProvince,China (N32°11′45.80″,E119°23′45.80″).MaintenanceconditionsfortheM.albatreeareasfollows: minimumandmaximumtemperatures(9–20℃),relativehumidity(40–60%),irrigation (rainfall),daylength(12h),andlightintensity(15,000lx).Additionally,informationabout theappearanceconditionsofthetreesincludesage(5year),diameter(5cm),andheight (2m).BTandY20havecharacteristicssuchaslargerinflorescences,agreaternumberof flowers,andwideadaptabilityandaretypicalcultivarsofM.alba.Therefore,thesetwo purifiedlineswereusedasthestudymaterialstocomparethedifferencesingeneexpressionbetweenmaleandfemaleinflorescencesofM.alba.Maleandfemaleinflorescence (1.3mfromtheground)sampleswerecollected(thatis,maleinflorescencefromBTpurifiedlinesandfemaleinflorescencefromthepurifiedY20lines).Thesampleswereamixtureofdifferentinflorescencesfromthesameplant,andsixtrees(3fromBTandY20each) wereselectedforsampling,representingthreebiologicalreplicates.Inall,atotalofsix catkinsweresampled.Inflorescencesamplingwasconductedfrom9:00p.mto11:00p.m. Onlyfullydevelopedadultcatkins(0.8–1.5cmlongand0.3–0.5cmindiameter)werecollected,whenthecatkinswerefullyexpanded,butbeforetheflowershadopened(Figure 1).WholecatkinswerecutandimmediatelyfrozeninliquidnitrogenforRNAextraction.  Figure1.PhysicalappearanceofcollectedadultMorusalbacatkins.(a)Femalecatkinswithprotrudingpistillateflowers.(b)Malecatkinswithprotrudingstaminateflowers. 2.2.RNAExtractionandIlluminaSequencing UsingtheRNAisoPlusreagent(Takara,Shanghai,China),reversetranscriptionand totalRNAextractionwereperformedoneachsampledcatkininaccordancewiththemanufacturer’sinstructions.UsinganAgilent2100Bioanalyzer(AgilentTechnologies,Palo Alto,CA,USA)andNanoDrop1000spectrophotometer(IMPLEN,CA,USA),thequantity andpurityoftotalRNAwereassessed.Furthermore,1µgoftotalRNApersamplewas usedtoconstructanRNAsequencelibraryusingtheNEBNext ® UltraTMRNALibrary PrepKitforIllumina(NEB,USA)accordingtothemanufacturer’sspecifications.Subsequently,transcriptomesequencingwasconductedusinganIlluminaNovaseq2500platform,whichproducedpaired-endreadsof150bp. Figure 1. Physical appearance of collected adult Morus alba catkins. (a) Female catkins with protruding pistillate flowers. (b) Male catkins with protruding staminate flowers. 2.2. RNA Extraction and Illumina Sequencing Using the RNAiso Plus reagent (Takara, Shanghai, China), reverse transcription and total RNA extraction were performed on each sampled catkin in accordance with the manufacturer’s instructions. Using an Agilent 2100 Bioanalyzer (Agilent Technologies, Palo Alto, CA, USA) and NanoDrop 1000 spectrophotometer (IMPLEN, Westlake Village, CA, USA), the quantity and purity of total RNA were assessed. Furthermore, 1 µ g of total RNA per sample was used to construct an RNA sequence library using the NEBNext ® UltraTM RNA Library Prep Kit for Illumina (NEB, Ipswich, MA, USA) according to the manufacturer’s specifications. Subsequently, transcriptome sequencing was conducted using an Illumina Novaseq 2500 platform (San Diego, CA, USA), which produced paired-end reads of 150 bp. Biology 2024,13, 622 4 of 17 2.3. Quality Control of Transcriptome Data and Expression Analysis To obtain clean RNA-seq data, the raw RNA-seq reads from each sample were processed using in-house perl scripts to trim adapters and low-quality bases from the ends of the reads, with further filtering of reads containing poly-N. The clean, high-quality reads for each sample were then mapped to the M. notabilis reference genome [ 20 ] using Hisat2 v2.0.5, and the mapped reads were assembled using StringTie (v1.3.3b) for novel transcript prediction. Fragments per kilobase million (fPKM) was used to estimate gene expression levels from each gene according to the number of reads mapped to that gene region using featureCounts v1.5.0-p3. Differential expression analysis was performed to detect the significantly differentially expressed genes (DEGs) between the male and female catkins using the DESeq2 R package (v1.16.1) [ 26 ]. Furthermore, p-values were adjusted based on the false discovery rate (FDR) using the Benjamini and Hochberg approach. Genes with a fold change difference in expression (|log2(FoldChange)| > 0) or (|log2(FoldChange)| < 0) and adjusted p-values (p-adj) < 0.05 were considered as DEGs. Among these DEGs, the value of (|log2(FoldChange)| > 0) was defined as male-biased genes and (|log2(FoldChange)| < 0) as female-biased genes. The extremes were considered as genes expressed only in females and were defined as female-limited genes, and the opposite was defined as male-limited genes. Gene ontology (GO) enrichment analysis was performed on the DEGs using the ClusterProfiler R package [ 27 ]. GO terms with padj < 0.05 were considered significantly enriched by DEGs. The ClusterProfiler R package was also used to perform DEG enrichment statistics in the KEGG pathway [ 27 ]. Pop_tri_v3 (Populus trichocarpa) was used as a reference (https://plants.ensembl.org/Populus_trichocarpa/Info/Index, accessed on 15 October 2023). 2.4. RT-qPCR Verification Based on the KEGG and GO enrichment results of the DEGs, 10 female-biased/-limited expression genes enriched in growth and development pathways and 10 male-biased/- limited expression genes enriched in defense-related pathways (Table S1) were randomly selected for RT-qPCR analysis. The RNA samples of flower buds used for RT-qPCR were taken from the same individuals as the samples used for transcriptome sequencing. A total of 20 pairs of gene-specific primers were constructed based on their sequences in the reference genome utilizing Primer 6 software (Table S1), and cDNA synthesis was performed using M-MLV reverse transcriptase (RTase) (Takara, Beijing, China), with 1 µ g RNA samples as the template. The cDNA solution was diluted 5-fold, and 1 µ L of cDNA was used as the template to perform gene validation using the SYBR Green RT-PCR protocol (Roche, Indianapolis, IN, USA). The reaction system contained SYBR qPCR 2 × Taq Mix (10 µ L), forward primer (1 µ L), reverse primer (1 µ L), ddH2O (7 µ L) and cDNA (1 µ L), and the β -actin gene (as a normalization control); the primer of the β -actin gene is shown in Table S1. All reactions were performed in three biological replicates using three templates. The procedure for the RT-qPCR was as follows: 95 ◦ C for 10 min, followed by a cycle program (denaturation: 95 ◦ C for 10 s, annealing: 50 ◦ C for 10 s, and extension: 70 ◦ C for 10 s) for 45 times. The mean relative expression of each gene was normalized to the reference gene, β-actin, and calculated using the 2−∆∆Ct method [28]. 2.5. Analysis of Adaptive Evolution of Sex-Biased and Unbiased Genes To test whether the sex-biased genes expressed only in males or females and the sexbiased DEGs had undergone significant positive selection, we used the evolutionary ratio (dN/dS) for all M. alba genes using paml v4.9 [ 29 ]. For each M. alba gene, the nucleotide sequences were translated into protein sequences, and their coding sequences were searched in TBtools v11.0.2 [ 30 ]. As M. alba shares a common ancestor with M. notabilis and diverged with Prunus persica, which is a hermaphrodite [ 21 ], we compared the coding sequences from our M. alba transcriptome data with the homologous genes in the M. notabilis genome [ 20 ]. Moreover, the sequences from the P. persica genome were used as an outgroup. Based on the bit score, e-value, and local alignment length (Supplementary Materials, R script), we Biology 2024,13, 622 5 of 17 used blastn 2.12.0 to identify the optimal multidirectional matching homologs to M. alba from M. notabilis v2.0 [ 20 ] and P. persica (Peach v1.0) [ 31 ] genome annotation data. The coding and protein sequences of the three species alignments were extracted and compared using ParaAT v1.0 [ 32 ]. We estimated the dN/dS values for the M. alba and M. notabilis lineages with codeml (PAML v4.9) [ 29 ], using a model with two or more dN/dS ratios for branches, runmode as user tree, and ndata as 13,532, and other parameters were set as the default [ 33 ]. We excluded 4980 alignments from the analysis due to unrealistically high divergence (S ×dS + N ×dN > 15% of alignment) [7]. In addition, we conducted a comparison of the dN/dS values for the evolutionary lineage leading to M. alba and M. notabilis, specifically focusing on genes that exhibit sexbiased or sex-limited expression patterns, as well as genes that do not display such biases. Kruskal–Wallis tests were used to compare the evolutionary ratio (dN/dS) between the sex-biased and unbiased genes, which is also used in Sanderson et al. (2019). Initially, a non-parametric Kruskal–Wallis test was employed to assess the differences in median values of dN/dS between sex-biased genes and a randomly selected unbiased genes (Kruskal.test and dunn.test in R to perform the Kruskal–Wallis test and post hoc tests respectively). To ascertain the rate of neutral evolution for these genes, a comparative analysis was conducted on the median dS values of sex-biased genes and a randomly selected set of 3000 unbiased genes [ 7 ]. Additionally, we conducted a comparison of the 95% quantile values of dN/dS across three categories in order to assess whether genes exhibiting sex-biased or sex-limited characteristics were more prone to displaying extreme values compared to the distribution of 95% quantile values derived from 5000 bootstrap samples of unbiased genes. This analysis was performed using the boot package in the R programming language (https://CRAN. R-project.org/package=boot, accessed on 12 August 2024). To explore the pathways and functions of positively selected DEGs, KEGG and GO enrichment analysis for the positively selected DEGs were performed using the clusterProfiler R package (p-adj < 0.05). 3. Results 3.1. Sequencing Quality Control and Reference Genome Mapping Six libraries were subjected to Illumina sequencing utilizing samples from male and female flower buds. Each sample was prepared in three repetitions. A combined sum of 305,642,716 raw sequencing reads were acquired from the flower buds of BT (male) and Y20 (female). After trimming adapter and poor-quality sequences, 294,445,006 clean reads were obtained. Among these reads, 139,671,832 (95.6%) were clean reads from male flower buds and 154,773,174 (97.0%) were clean reads from female flower buds. For all the libraries, the average Q30 was more than 92% and the GC content was more than 40% (Table 1). Furthermore, more than 64.5% of the clean reads were mapped to the reference genome (M. notabilis). A total of 21,155 unigenes were obtained from the assembly, with the longest gene being 16,567 bp and the shortest having 72 bp. Table 1. Summary of quality preprocessing of RNA sequencing data. Sample ID Raw Reads Clean Reads Clean Bases Error Rate Q30 GC% Mapping Ratio MBT-1 52,365,500 50,926,984 7.64 G 0.03 93.12 44.31 64.6 MBT-2 44,385,448 41,784,712 6.27 G 0.03 93.55 42.81 58.4 MBT-3 49,325,358 46,960,136 7.04 G 0.03 92.19 40.73 72.8 FY20-1 56,900,778 54,807,454 8.22 G 0.03 92.27 44.23 55.8 FY20-2 48,518,050 47,476,836 7.12 G 0.03 93.46 44.66 66.9 FY20-3 54,147,582 52,488,884 7.87 G 0.03 92.38 43.44 68.5 MBT and FY20: Male and female mulberry (Morus alba) flower buds used in RNA sequencing. MBT-1, MBT-2, and MBT-3 and FY20-1, FY20-2, and FY20-3: three replicates of male and female mulberry (M. alba) flower buds, respectively. Biology 2024,13, 622 6 of 17 3.2. Differentially Expressed Genes between Male and Female Flower Buds Following the transcriptome assembly and annotation, differential expression analysis was performed. Surprisingly, a total of 1,543 genes were differentially expressed in male and female flower buds based on the criteria stated in Section 2.3. Out of these DEGs, 702 genes were significant in female flower buds and 841 genes were significant in male flower buds (Figure 2). A total of 19,612 genes exhibited a transcript abundance difference that did not meet the criteria (Figure 2c). These findings indicate the presence of biologically significant levels of gene expression dimorphism. Genes that were expressed exclusively in one sex type were defined as sex-limiting genes (DEGs showing some expression in one sex type but absolutely zero counts in the other). A total of 399 genes exhibited expression that was sex-limited. Out of the total, 333 genes (83.5%) were exclusively expressed in male flowers, while 66 genes (16.5%) were exclusively expressed in female flowers, suggesting that the proportion of male-biased genes was higher than that of female-biased genes. The number of genes exhibiting male bias was marginally higher compared to those exhibiting female bias. Biology2024,13,xFORPEERREVIEW6of17 3.2.DifferentiallyExpressedGenesbetweenMaleandFemaleFlowerBuds Followingthetranscriptomeassemblyandannotation,differentialexpressionanalysiswasperformed.Surprisingly,atotalof1,543genesweredifferentiallyexpressedin maleandfemaleflowerbudsbasedonthecriteriastatedinSection2.3.OutoftheseDEGs, 702genesweresignificantinfemaleflowerbudsand841genesweresignificantinmale flowerbuds(Figure2).Atotalof19,612genesexhibitedatranscriptabundancedifference thatdidnotmeetthecriteria(Figure2c).Thesefindingsindicatethepresenceofbiologicallysignificantlevelsofgeneexpressiondimorphism.Genesthatwereexpressedexclusivelyinonesextypeweredefinedassex-limitinggenes(DEGsshowingsomeexpression inonesextypebutabsolutelyzerocountsintheother).Atotalof399genesexhibited expressionthatwassex-limited.Outofthetotal,333genes(83.5%)wereexclusivelyexpressedinmaleflowers,while66genes(16.5%)wereexclusivelyexpressedinfemale flowers,suggestingthattheproportionofmale-biasedgeneswashigherthanthatoffemale-biasedgenes.Thenumberofgenesexhibitingmalebiaswasmarginallyhighercomparedtothoseexhibitingfemalebias. Figure2.ComparisonsofgeneexpressionbetweenmaleandfemaleflowerbudsofM.alba.(a) Shareddifferentiallyexpressedgenesbetweenmaleandfemaleflowerbuds.(b)Geneexpression comparisonbetweenmaleandfemaleflowerbuds.(c)Volcanoplotonsex-biasedgenesthatwere expressedsignificantlydifferentinmaleandfemaleflowerbuds. Additionally,theextentofdifferentialexpression,asmeasuredby(|log2FC|),was significantlyhigherforgeneswithmalebiascomparedtothosewithfemalebias(Figure 3a,Wilcoxonranksumtestp<0.01,R4.1.0).Themeanexpressionlevelsofmale-biased geneswerefoundtobeconsiderablylowercomparedtotheaverageexpressionlevelsof female-biasedgenes(Figure3b;Wilcoxonranksumtestp<0.01,R4.1.0),althoughfemale expressionformale-biasedgeneswassignificantlyhigherthanmaleexpressionforfemale-biasedgenes(Figure3b,Wilcoxonranksumtestp<0.01,R4.1.0).Thefindingsof thisanalysisindicatethatgeneswithabiastowardsfemalesmayoriginatefromheightenedexpressioninfemaleinflorescencesandreducedexpressioninmaleinflorescences. Figure 2. Comparisons of gene expression between male and female flower buds of M. alba. (a) Shared differentially expressed genes between male and female flower buds. (b) Gene expression comparison between male and female flower buds. (c) Volcano plot on sex-biased genes that were expressed significantly different in male and female flower buds. Additionally, the extent of differential expression, as measured by (|log2FC|), was significantly higher for genes with male bias compared to those with female bias (Figure 3a, Wilcoxon rank sum test p< 0.01, R 4.1.0). The mean expression levels of male-biased genes were found to be considerably lower compared to the average expression levels of female-biased genes (Figure 3b; Wilcoxon rank sum test p< 0.01, R 4.1.0), although female expression for male-biased genes was significantly higher than male expression for female-biased genes (Figure 3b, Wilcoxon rank sum test p< 0.01, R 4.1.0). The findings of this analysis indicate that genes with a bias towards females may originate from heightened expression in female inflorescences and reduced expression in male inflorescences. Additionally, it is seen that the bias towards males predominantly stems from the upregulation of female expression. Biology 2024,13, 622 7 of 17 Biology2024,13,xFORPEERREVIEW7of17   Additionally,itisseenthatthebiastowardsmalespredominantlystemsfromtheupregulationoffemaleexpression.  Figure3.Sex-biasedgeneexpressioninM.alba.(a)TheproportionandrangeofDEGsandunbiased genesinM.albacatkins.(b)Theupper-halfshowsthecomparisonofmaleandfemalecatkindifferentiallyexpressiongenes(DEGs)atdifferentsexbiasfoldchangethresholdsforallevaluatedcatkin male-biasedandfemale-biasedgenes.ThenumbersinthebracketsrepresentthenumberofDEGs ineachcatkincategory.Thelower-halfshowstheaveragemaleandfemalecatkinexpressionofall genesatdifferentsexbiasfoldchangethresholdsforalltheevaluatedcatkinmale-biasedandfemale-biasedgenes.SignificantdifferencesbetweenmaleandfemaleexpressionbasedonWilcoxon ranksumtestsaredenotedbyasteriskatalphalevel(p<0.0001).Thesex-biasedgenesmentioned hereincludesex-limitedgenes.BTandY20arethemaleandfemalecatkins. 3.3.FunctionalAnnotationAssociatedwithDEGs TodeterminewhethertheseDEGsareassociatedwiththedevelopmentofsex-specifictraits,GOandKEGGanalysiswasusedtoidentifypatternsofenrichmentbetween theDEGsinmaleandfemaleflowerbuds.AccordingtotheGOenrichmentresults,the mostfemale-biasedgeneswereannotatedinnucleicacidbindingtranscriptionfactoractivity(GO:0003700),transcriptionfactoractivity,andsequence-specificDNAbinding (GO:0001071)(Figure4a,b).Themostmale-biasedgeneswereannotatedforcarbohydrate Figure 3. Sex-biased gene expression in M. alba. (a) The proportion and range of DEGs and unbiased genes in M. alba catkins. (b) The upper-half shows the comparison of male and female catkin differentially expression genes (DEGs) at different sex bias fold change thresholds for all evaluated catkin male-biased and female-biased genes. The numbers in the brackets represent the number of DEGs in each catkin category. The lower-half shows the average male and female catkin expression of all genes at different sex bias fold change thresholds for all the evaluated catkin male-biased and female-biased genes. Significant differences between male and female expression based on Wilcoxon rank sum tests are denoted by asterisk at alpha level (p< 0.0001). The sex-biased genes mentioned here include sex-limited genes. BT and Y20 are the male and female catkins. 3.3. Functional Annotation Associated with DEGs To determine whether these DEGs are associated with the development of sex-specific traits, GO and KEGG analysis was used to identify patterns of enrichment between the DEGs in male and female flower buds. According to the GO enrichment results, the most female-biased genes were annotated in nucleic acid binding transcription factor activity (GO:0003700), transcription factor activity, and sequence-specific DNA binding (GO:0001071) (Figure 4a,b). The most male-biased genes were annotated for carbohydrate metabolism processes (GO:0005975). To further understand the enrichment pathways of male-/female-biased genes, these DEGs were analyzed by KEGG pathway analysis. Male-biased genes were enriched in 111 KEGG pathways and significantly enriched in the phenylpropanoid biosynthesis pathway (pop00940). Female-biased genes were as- Biology 2024,13, 622 8 of 17 signed to 101 KEGG pathways and were mainly enriched in the plant hormone signal pathway (pop04075) (Figure 4c,d). The role of plant hormones in plant sex determination has been demonstrated in many studies [ 34 , 35 ]. According to the result of KEGG enrichment analysis, 22 DEGs (10 female-biased and 12 male-biased) were enriched in plant hormone signaling, including auxin, cytokinin, salicylic acid, jasmonic acid, and abscisic acid signaling pathways (Figure 5a). Interestingly, the male-biased genes were mainly annotated in the abscisic acid and jasmonic acid signaling pathways, whereas the female-biased genes were mainly annotated in the cytokinin and auxin signaling pathways. The screened DEGs contained 88 transcription factors (TFs) that could be successfully annotated for the KEGG pathway, including 33 female-biased TFs and 55 male-biased TFs. Biology2024,13,xFORPEERREVIEW8of17 metabolismprocesses(GO:0005975).Tofurtherunderstandtheenrichmentpathwaysof male-/female-biasedgenes,theseDEGswereanalyzedbyKEGGpathwayanalysis.Malebiasedgeneswereenrichedin111KEGGpathwaysandsignificantlyenrichedinthephenylpropanoidbiosynthesispathway(pop00940).Female-biasedgeneswereassignedto 101KEGGpathwaysandweremainlyenrichedintheplanthormonesignalpathway (pop04075)(Figure4c,d).Theroleofplanthormonesinplantsexdeterminationhasbeen demonstratedinmanystudies[34,35].AccordingtotheresultofKEGGenrichmentanalysis,22DEGs(10female-biasedand12male-biased)wereenrichedinplanthormonesignaling,includingauxin,cytokinin,salicylicacid,jasmonicacid,andabscisicacidsignaling pathways(Figure5a). Interestingly,themale-biasedgenesweremainlyannotatedintheabscisicacidand jasmonicacidsignalingpathways,whereasthefemale-biasedgenesweremainlyannotatedinthecytokininandauxinsignalingpathways.ThescreenedDEGscontained88 transcriptionfactors(TFs)thatcouldbesuccessfullyannotatedfortheKEGGpathway, including33female-biasedTFsand55male-biasedTFs. Figure4.FunctionalanalysisoftheDEGsinM.albamaleandfemalecatkins.(a)GOenrichment termsofmale-biasedexpressiongenes.(b)GOenrichmenttermsoffemale-biasedexpressiongenes. Figure 4. Functional analysis of the DEGs in M. alba male and female catkins. (a) GO enrichment terms of male-biased expression genes. (b) GO enrichment terms of female-biased expression genes. (c) KEGG-enriched differential expression genes of male-biased expression genes. (d) KEGG-enriched differential expression genes of female-biased expression genes. BP; biological process. CC; cellular component. MF; molecular function. Biology 2024,13, 622 9 of 17 Biology2024,13,xFORPEERREVIEW9of17   (c)KEGG-enricheddifferentialexpressiongenesofmale-biasedexpressiongenes.(d)KEGG-enricheddifferentialexpressiongenesoffemale-biasedexpressiongenes.BP;biologicalprocess.CC; cellularcomponent.MF;molecularfunction.  Figure5.AheatmapanalysisoftheDEGsinM.albamaleandfemalecatkins.(a)AheatmapdiagramoftheDEGsbetweenmaleandfemaleflowerbudsofM.albainthephytohormonesignaling pathways:abscisicacid(i),jasmonicacid(ii),salicylicacid(iii),cytokinin(iv),andauxin(v).(b)The transcriptionfactor(TF)familiesdifferentiallyexpressedinmale-biasedgenes.(c)TheTFfamilies differentiallyexpressedinmale-biasedandfemale-biasedgenes.Thecolorscalerepresentsthelog10transformedFPKMvalue.Thesex-biasedgenesmentionedhereincludesex-limitedgenes.The male-biasedgenesareupregulatedandfemale-biasedgenesaredownregulatedbasedonthefigure legend. Amongthem,thefemale-biased/-limitedTFswereannotatedfor29KEGGpathways, whichwererelatedtogrowthanddevelopment.Theyincludemetabolicpathways (pop01100),biosynthesisofsecondarymetabolites(pop01110),starchandsucrosemetabolism(pop00500),andsoon.However,themale-biased/-limitedTFswereannotatedfor Figure 5. A heat map analysis of the DEGs in M. alba male and female catkins. (a) A heat map diagram of the DEGs between male and female flower buds of M. alba in the phytohormone signaling pathways: abscisic acid (i), jasmonic acid (ii), salicylic acid (iii), cytokinin (iv), and auxin (v). (b) The transcription factor (TF) families differentially expressed in male-biased genes. (c) The TF families differentially expressed in male-biased and female-biased genes. The color scale represents the log 10 -transformed FPKM value. The sex-biased genes mentioned here include sex-limited genes. The male-biased genes are upregulated and female-biased genes are downregulated based on the figure legend. Among them, the female-biased/-limited TFs were annotated for 29 KEGG pathways, which were related to growth and development. They include metabolic pathways (pop01100), biosynthesis of secondary metabolites (pop01110), starch and sucrose metabolism (pop00500), and so on. However, the male-biased/-limited TFs were annotated for 66 KEGG pathways, including metabolic pathways (pop01100), biosynthesis of secondary metabolites (pop01110), glycerophospholipid metabolism (pop00564), plant hormone signal transduction (pop04075), glycolysis/Gluconeogenesis (pop00010), glycerolipid metabolism (pop00561), and fatty acid metabolism (pop01212). TFs’ regulation of gene expression and the expression profiles of TF genes between male and female flower buds are shown as a heat map (Figure 5b,c). Biology 2024,13, 622 16 of 17 4. Ellegren, H.; Parsch, J. The evolution of sex-biased genes and sex-biased gene expression. Nat. Rev. Genet. 2007,8, 689–698. [CrossRef] 5. Harrison, P.W.; Wright, A.E.; Zimmer, F.; Dean, R.; Montgomery, S.H.; Pointer, M.A.; Mank, J.E. Sexual selection drives evolution and rapid turnover of male gene expression. Proc. Natl. Acad. Sci. USA 2015,112, 4393–4398. [CrossRef] 6. Zemp, N.; Tavares, R.; Muyle, A.; Charlesworth, D.; Marais, G.A.B.; Widmer, A. Evolution of sex-biased gene expression in a dioecious plant. Nat. Plants 2016,2, 16168. [CrossRef] [PubMed] 7. Sanderson, B.J.; Wang, L.; Tiffin, P.; Wu, Z.; Olson, M.S. Sex-biased gene expression in flowers, but not leaves, reveals secondary sexual dimorphism in Populus balsamifera. New Phytol. 2019,221, 527–539. [CrossRef] 8. Barrett, S.C.; Hough, J. Sexual dimorphism in flowering plants. J. Exp. Bot. 2013,64, 67–82. [CrossRef] 9. Zajitschek, F.; Bonduriansky, R.; Zajitschek, S.R.; Brooks, R.C. Sexual dimorphism in life history: Age, survival, and reproduction in male and female field crickets Teleogryllus commodus under seminatural conditions. Am. Nat. 2009,173, 792–802. [CrossRef] [PubMed] 10. Cervantes, C.; Alvarez, A.; Cuevas, E. Small but attractive: Female-biased nectar production and floral visitors in a dimorphic shrub. Plant Biol. 2018,20, 160–164. [CrossRef] 11. Tsuji, K.; Kobayashi, K.; Hasegawa, E.; Yoshimura, J. Dimorphic flowers modify the visitation order of pollinators from male to female flowers. Sci. Rep. 2020,10, 9965. [CrossRef] 12. Gershenzon, J.; Ullah, C. Plants protect themselves from herbivores by optimizing the distribution of chemical defenses. Proc. Natl. Acad. Sci. USA 2022,119, e2120277119. [CrossRef] [PubMed] 13. Wakabayashi, K.; Oguro, M.; Itagaki, T.; Sakai, S. Floral-induced and constitutive defense against florivory: A comparison of chemical traits in 12 herb species. Plant Ecol. 2018,219, 985–997. [CrossRef] 14. Gao, L.; Yu, G.; Hu, F.; Li, Z.; Li, W.; Peng, C. The patterns of male and female flowers in flowering stage may not be optimal resource allocation for fruit and seed growth. Plants 2021,10, 2819. [CrossRef] [PubMed] 15. Zhao, Y. Auxin Biosynthesis and Its Role in Plant Development. Annu. Rev. Plant Biol. 2010,61, 49–64. [CrossRef] [PubMed] 16. Liu, J.; Chen, L.-Y.; Zhou, P.; Liao, Z.; Lin, H.; Yu, Q.; Ming, R. Sex biased expression of hormone related genes at early stage of sex differentiation in papaya flowers. Hortic. Res. 2021,8, 147. [CrossRef] [PubMed] 17. Schenck, C.A.; Maeda, H.A. Tyrosine biosynthesis, metabolism, and catabolism in plants. Phytochemistry 2018,149, 82–102. [CrossRef] [PubMed] 18. Arunkumar, R.; Josephs, E.B.; Williamson, R.J.; Wright, S.I. Pollen-specific, but not sperm-specific, genes show stronger purifying selection and higher rates of positive selection than sporophytic genes in Capsella grandiflora.Mol. Biol. Evol. 2013,30, 2475–2486. [CrossRef] 19. Xuan, Y.; Ma, B.; Li, D.; Tian, Y.; Zeng, Q.; He, N. Chromosome restructuring and number change during the evolution of Morus notabilis and Morus alba.Hortic. Res. 2022,9, uhab030. [CrossRef] [PubMed] 20. He, N.; Zhang, C.; Qi, X.; Zhao, S.; Tao, Y.; Yang, G.; Lee, T.-H.; Wang, X.; Cai, Q.; Li, D.; et al. Draft genome sequence of the mulberry tree Morus notabilis.Nat. Commun. 2013,4, 2445. [CrossRef] [PubMed] 21. Jiao, F.; Luo, R.; Dai, X.; Liu, H.; Yu, G.; Han, S.; Lu, X.; Su, C.; Chen, Q.; Song, Q.; et al. Chromosome-level reference genome and population genomic analysis provide insights into the evolution and improvement of domesticated mulberry (Morus alba). Mol. Plant 2020,13, 1001–1012. [CrossRef] [PubMed] 22. Jaiswal, V.; Kumar, A. Activity and isozymes of peroxidase in Morus nigra L. during sex differentiation. Z. Pflanzenphysiol. 1981, 102, 299–302. [CrossRef] 23. Lal, M.; Jaiswal, V.S. Modification of flower sex and acid phosphatase activity by phthalimides in female plants of Morus nigra L. Plant Growth Regul. 1988,7, 29–37. [CrossRef] 24. Jolly, M.S.; Dandin, S.B.; Ravindran, S.; Kumar, R. Sexual polymorphism in the genus Morus L. Proc. Plant Sci. 1986,96, 315–320. [CrossRef] 25. Atsumi, R.; Nishihara, R.; Tarora, K.; Urasaki, N.; Matsumura, H. Identification of dominant genetic markers relevant to male sex determination in mulberry (Morus alba L.). Euphytica 2019,215, 187. [CrossRef] 26. Anders, S.; Huber, W. Differential expression analysis for sequence count data. Nat. Preced. 2010, 1. [CrossRef] 27. Yu, G.; Wang, L.-G.; Han, Y.; He, Q.-Y. clusterProfiler: An R package for comparing biological themes among gene clusters. OMICS J. Integr. Biol. 2012,16, 284–287. [CrossRef] [PubMed] 28. Livak, K.J.; Schmittgen, T.D. Analysis of relative gene expression data using real-time quantitative PCR and the 2 −∆∆CT method. Methods 2001,25, 402–408. [CrossRef] [PubMed] 29. Yang, Z. PAML 4: Phylogenetic analysis by maximum likelihood. Mol. Biol. Evol. 2007,24, 1586–1591. [CrossRef] [PubMed] 30. Chen, C.J.; Chen, H.; Zhang, Y.; Thomas, H.R.; Frank, M.H.; He, Y.H.; Xia, R. TBtools: An integrative toolkit developed for interactive analyses of big biological data. Mol. Plant 2020,13, 1194–1202. [CrossRef] [PubMed] 31. International Peach Genome Initiative; Verde, I.; Abbott, A.G.; Scalabrin, S.; Jung, S.; Shu, S.; Marroni, F.; Zhebentyayeva, T.; Dettori, M.T.; Grimwood, J.; et al. The high-quality draft genome of peach (Prunus persica) identifies unique patterns of genetic diversity, domestication and genome evolution. Nature Genet. 2013,45, 487–494. [CrossRef] [PubMed] 32. Zhang, Z.; Xiao, J.; Wu, J.; Zhang, H.; Liu, G.; Wang, X.; Dai, L. ParaAT: A parallel tool for constructing multiple protein-coding DNA alignments. Biochem. Biophys. Res. Commun. 2012,419, 779–781. [CrossRef] 33. Yang, Z.; Rannala, B. Molecular phylogenetics: Principles and practice. Nat. Rev. Genet. 2012,13, 303–314. [CrossRef] [PubMed] Biology 2024,13, 622 17 of 17 34. Devani, R.S.; Sinha, S.; Banerjee, J.; Sinha, R.K.; Bendahmane, A.; Banerjee, A.K. De novo transcriptome assembly from flower buds of dioecious, gynomonoecious and chemically masculinized female Coccinia grandis reveals genes associated with sex expression and modification. BMC Plant Biol. 2017,17, 241. [CrossRef] [PubMed] 35. Yuan, Z.; Zhang, D. Roles of jasmonate signalling in plant inflorescence and flower development. Curr. Opin. Plant Biol. 2015,27, 44–51. [CrossRef] [PubMed] 36. McKown, A.D.; Klápštˇe, J.; Guy, R.D.; Soolanayakanahally, R.Y.; La Mantia, J.; Porth, I.; Skyba, O.; Unda, F.; Douglas, C.J.; El-Kassaby, Y.A.; et al. Sexual homomorphism in dioecious trees: Extensive tests fail to detect sexual dimorphism in Populus.Sci. Rep. 2017,7, 1831. [CrossRef] [PubMed] 37. Robinson, K.M.; Delhomme, N.; Mähler, N.; Schiffthaler, B.; Önskog, J.; Albrectsen, B.R.; Ingvarsson, P.K.; Hvidsten, T.R.; Jansson, S.; Street, N.R. Populus tremula (European aspen) shows no evidence of sexual dimorphism. BMC Plant Biol. 2014,14, 276. [CrossRef] [PubMed] 38. Harkess, A.; Mercati, F.; Shan, H.; Sunseri, F.; Falavigna, A.; Leebens-Mack, J. Sex-biased gene expression in dioecious garden asparagus (Asparagus officinalis). New Phytol. 2015,207, 883–892. [CrossRef] [PubMed] 39. Bourdenx, B.; Bernard, A.; Domergue, F.; Pascal, S.; Léger, A.; Roby, D.; Pervent, M.; Vile, D.; Haslam, R.P.; Napier, J.A.; et al. Overexpression of Arabidopsis ECERIFERUM1 promotes wax very-long-chain alkane biosynthesis and influences plant response to biotic and abiotic stresses. Plant Physiol. 2011,156, 29–45. [CrossRef] [PubMed] 40. Strommer, J. The plant ADH gene family. Plant J. 2011,66, 128–142. [CrossRef] [PubMed] 41. Wang, H.; Dong, Q.; Duan, D.; Zhao, S.; Li, M.; van Nocker, S.; Ma, F.; Mao, K. Comprehensive genomic analysis of the TYROSINE AMINOTRANSFERASE (TAT) genes in apple (Malus domestica) allows the identification of MdTAT2 conferring tolerance to drought and osmotic stresses in plants. Plant Physiol. Biochem. 2018,133, 81–91. [CrossRef] Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.