Introduction

Sexual dimorphism presents different morphological, ecological and physiological traits in gonochoristic animals and dioecious plants, despite male and female individuals sharing the same genome except for sex chromosomes (or sex-determining loci) (Mank, 2009; Barrett and Hough, 2013). Such sexual dimorphisms usually arise from differential expression of genes in both sexes, i.e., sex-biased genes (including sex-specific genes exclusively in one sex) that are located on autosomal chromosomes and sex chromosomes/or sex-determining regions (Ellegren and Parsch, 2007; Parsch and Ellegren, 2013; Charlesworth, 2018). Recently, some studies have begun to explore the strength and impact of evolutionary forces that shape different sexually dimorphic traits through sex-biased gene expression (Mank, 2009; Rowe et al., 2018). Previous studies revealed that sex-biased gene expressions have facilitated the evolution of sexual dimorphisms in many animals. However, the extent to which this occurs varies greatly among taxa, tissues, and developmental stages (Mank, 2017; Khodursky et al., 2020; Hsu et al., 2020; Lichilin et al., 2021; Toubiana et al., 2021; Djordjevic et al., 2022). Unlike most animals, the vast majority (~90%) of angiosperms are hermaphroditic plants, while only a small fraction (~5%) are dioecious plants in which individuals have exclusively male or female reproductive organs (Renner, 2014). Most dioecious plants do not possess heteromorphic sex chromosomes. Furthermore, sexual dimorphism in dioecious plants is less common and less conspicuous than that in almost animals (Barrett and Hough, 2013). Hence, the study of sex-biased gene expression is of great interest to plant evolutionary biologists, as it is essential to understand the evolution of sexual dimorphism in dioecious plants.

A common pattern that has emerged from previous studies is that sex-biased genes, particularly male-biased genes, tend to evolve rapidly in protein sequence (the ratio of non-synonymous to synonymous substitutions, d_N/d_S) compared to unbiased genes (Ellegren and Parsch, 2007; Grath and Parsch, 2016). The rapid evolution of male-biased genes was first observed in Drosophila melanogaster (Zhang et al., 2004; Zhang and Parsch, 2005) and has been supported by recent investigations in a wider range of animals (Pro-schel et al., 2006; Mank et al., 2007; Mank, 2017; Papa et al., 2017; Catalan et al., 2018; Toubiana et al., 2021). In recent years, there have been growing studies on the expression dynamics and molecular evolutionary rates of sex-biased genes in flowering plants, including Silene latifolia (Zemp et al., 2016), Salix viminalis (Darolti et al., 2018), Mercurialis annua (Cossard et al., 2019), Populus balsamifera (Sanderson et al., 2019), and Leucadendron plants (Scharmann et al., 2021). However, despite such advances, the molecular evolution pattern of sex-biased genes in plants remains inconsistent among the studied plant species (Muyle, 2019; Veltsos, 2019). In dioecious plants such as Mercurialis annua and Leucadendron, Cossard et al., (2019) and Scharmann et al., (2021) found no significant differences in evolutionary rates of proteins among female-biased, male-biased and unbiased genes, although the expression of sex-biased genes was highly different from unbiased genes in leaves. Similar patterns have also been reported in dioecious Populus balsamifera, where evolutionary rates of male-biased, female-biased and unbiased genes did not differ in reproductive tissues (Sanderson et al., 2019). However, in the dioecious plant Salix viminalis, male-biased genes have significantly lower evolutionary rates of proteins than female-biased and unbiased genes in catkin tissues (Darolti et al., 2018). To our knowledge, only five studies have investigated the expression differences and protein evolutionary rates of sex-biased genes in dioecious angiosperms. Moreover, these studies only compared gene expression in vegetative vs. vegetative tissues and vegetative vs. reproductive tissues. Therefore, more studies and taxa are needed to explore the common patterns of sequence evolution in sex-biased genes, with more focus on comparing sex-biased gene expression in reproductive vs. reproductive tissues, e.g., different floral development stages in dioecious angiosperms.

Evolutionary dynamic analyses indicate that different evolutionary forces impact the rate of sequence evolution of sex-biased genes. These forces include positive selection, which promotes the spread and adaptive fixation of beneficial alleles; sexual selection, which results from male-male competition or female choice; and relaxed purifying selection, which reduces the removal of deleterious mutations (Grath and Parsch, 2016; Mank, 2017; Dapper and Wade, 2020). For example, in animal systems, particularly in Drosophia, the elevated sequence divergence rates of male-biased genes have often been interpreted as the signature of positive selection due to adaptive evolution, suggesting that sexual selection is the primary evolutionary forces (Pro-schel et al., 2006; Assis et al., 2012). However, studies in plants have rarely reported elevated rates of sex-biased genes, except for brown alga, in which female-biased and/or male-biased genes exhibiting higher evolutionary rates than unbiased genes suggests that rapid evolution is partly driven by adaptive evolution or sexual selection (Lipinska et al., 2015; Cossard et al., 2022; Hatchett et al., 2023).

An alternative explanation for the rapid evolution of sex-biased genes is a relaxation of purifying selection due to reduced constraints (Lahti et al., 2009; Dapper and Wade, 2020). In the model plant Arabidopsis thaliana, pollen genes were found to be evolving faster than sporophyte-specific genes due to relaxed purifying selection associated with the transition from outcrossing to selfing (Harrison et al., 2019). These trends were recently confirmed in Arabis alpina, which exhibits mating system variation across its distribution, suggesting that the efficacy of purifying selection on male gametophyte-expressed genes was significantly weaker in inbred populations (Gutierrez-Valencia et al., 2022). Together, these findings in plants reinforce the idea that both adaptive (e.g., positive selection, sexual selection) and non-adaptive (e.g., relaxed selection) evolutionary processes differentially impact the sequence evolution of sex-biased genes. Hence, investigating the potential contribution of selection forces to the emergence of specific evolutionary patterns of sex-biased genes within a focal species is of great interest.

In the family Cucurbitaceae, there are about 96 genera and 1,000 species, about 50% of species are dioecious, and 50% are monoecious (Schaefer and Renner, 2011). Phylogenetic analyses of Cucurbitaceae suggest that dioecy is the ancestral state of the family, but transitions frequently to monoecy (Zhang et al., 2006). Trichosanthes pilosa (2n = 22, Cucurbitaceae) is mainly distributed from Southwest, Southeast China to Japan, New Guinea and Western Australia. It was suggested to have originated in the late Miocene (c.13 Ma) (de Boer et al., 2012; 2015). Trichosanthes pilosa (synonym: T. ovigera) is a perennial dioecious vine that reproduces sexually and possesses a pair of heteromorphic sex chromosomes XX/XY (Ming et al., 2011). The species exhibits strong sexual dimorphism in the floral morphological structures and the life-history traits, such as racemose male flower and solitary female flower (Fig. 1), early-flowering, caducous male flower and late-flowering, long-lived female flower (Wu et al., 2011).

Fig. 1

To understand the evolution of sex-biased genes in dioecious Trichosanthes pilosa, we collected floral buds and mature flowers tissues from male and female individuals and characterized their expression profiles using Illumina RNA sequencing. Our primary objectives are to 1) compare expression divergences between males and females at two floral developmental stages; 2) explore whether there are differences in the evolutionary rates of proteins among female-biased, male-biased and unbiased genes; and if so, 3) determine the main selective forces that contribute to the differentiation of sequence evolution rates among gene categories.

Results

Transcriptome sequencing, de novo assembly and annotation

Using whole transcriptome shotgun sequencing, we sequenced 12 floral buds and mature flower tissues from females and males of dioecious T. pilosa and generated a total of nearly 276 million clean reads (Table S1). We then performed de novo assembly of transcripts from all the clean reads, followed by clustering and filtering analysis, resulting in 59,051 unigenes (Fig. S1A). To evaluate the quality of the assembled unigenes, we compared them against protein databases such as NR, KEGG, Swissport, PFAM, and GO using BLASTP and nucleotide database NT using BLASTN (Table S2). The e-value distribution of the best hits in NR database suggested that 47,241 unigenes (80%) had strong homology, with an e-value smaller than 1.0e-15 (Fig. S1B). The majority of unigenes were annotated by homologs in other species of Cucurbitaceae, such as Momordica charantia, Cucumis melo, Cucurbita pepo, Cucurbita moschata and Cucurbita maxima (Fig. S1C). Furthermore, we used BUSCO assessments based on embryophyta_odb9 database, which showed the completeness of the reference transcriptome at 78.8% (Table S2). Overall, our analyses suggested that we have generated high-quality reference transcriptomes.

Expression characteristics of sex-biased genes

We mapped the RNA-Seq reads of floral buds and mature flowers onto the reference transcriptome in dioecious T. pilosa, which resulted in approximately 75% read mappings per sample (Table S3). In floral buds, we identified 5,096 (9.50%) female-biased genes and 4,214 (7.86%) male-biased genes (Fig. 2A). In contrast, only 380 (0.70%) female-biased genes and 233 (0.43%) male-biased genes were detected in mature flowers (Fig. 2B). Using hierarchical clustering analysis, we evaluated different levels of gene expression across sexes and tissues, revealing strong sexual dimorphism (Fig. 2C). Furthermore, we observed that the number of sex-biased genes in floral buds was approximately 15 times higher than in mature flowers, indicating that sex-biased genes associated with sex differentiation and sexually dimorphic traits are predominantly expressed in floral buds. We also analyzed sex-specific genes that were exclusively expressed in floral buds and mature flowers of one sex. In floral buds, we found 253 out of 5,096 (4.96%) female-specific genes and 465 out of 4,214 (11.03%) male-specific genes. However, in mature flowers, we only identified 26 out of 380 (6.84%) female-specific genes and 52 out 233 (22.32%) male-specific genes.

Fig. 2

Tissue-biased gene expression

We compared the expression levels of transcripts in floral bud tissues and mature flower tissues within each sex to identify genes with tissue-biased expression. In male plant, the number of tissue-biased genes in mature flowers was 1,040 higher than that in floral buds (Fig. 3A and 3B). However, in female plant, the number of tissue-biased genes in mature flowers was only 536 more than that in floral buds (Fig. 3C and 3D). Our results indicated that males had a higher tissue-bias relative to females. We also identified sex-biased genes that were expressed in both types of tissues by comparing tissue-biased genes with male-biased and female-biased genes, respectively. Few of female-biased genes overlapped with tissue-biased genes, accounting for only 85 out of 5,096 (1.67%) in floral buds (Fig. 3C) and 5 out of 380 (1.32%) in mature flowers (Fig. 3D). However, a significant proportion of male-biased genes overlapped with tissue-biased genes, with 1,010 out of 4,214 (23.97%) in floral buds (Fig. 3A) and 145 out of 233 (62.23%) in mature flowers (Fig. 3B).

Fig. 3

Elevated protein evolutionary rates of male-biased genes in floral buds

We identified 1,145 female-biased, 343 male-biased, and 2,378 unbiased one-to-one orthologous groups (OGs) from floral buds in three additional species of Trichanthes, T. anguina, T. pilosa, T. kirilowii, together with Luffa cylindrica. Additionally, we detected 45 female-biased, 13 male-biased, and 3,782 unbiased one-to-one OGs from mature flowers in all four species. To quantify the rates of protein sequences, we separately calculated ω values for each sex-biased and unbiased orthologous gene using two branch models in floral buds and mature flowers (Figs. 4 and S2).

Fig. 4

In the results of ‘two-ratio’ branch model, the median of ω values in female-biased, male-biased and unbiased genes were 0.227, 0.257 and 0.230 in floral buds, respectively (Fig. 4A and Table S4). We observed that male-biased genes had a 13.22% and 11.74% higher median than female-biased and unbiased genes in floral buds, respectively. The difference in the distribution of ω values between female-biased versus male-biased genes (P = 0.0021) and male-biased versus unbiased genes (P = 0.0051) was statistically significant in Wilcoxon rank sum tests. However, we did not find a significant difference in ω values between female-biased and unbiased genes in floral buds (Wilcoxon rank sum test, P = 0.4618). In mature flowers, the median of ω values for female-biased, male-biased, and unbiased genes were 0.269, 0.177 and 0.231, respectively (Fig. 4B and Table S4). However, we found that there was no statistically significant difference in the distribution of ω values using Wilcoxon rank sum tests for female-biased versus male-biased genes (P = 0.0556), female-biased versus unbiased genes (P = 0.0796), and male-biased versus unbiased genes (P = 0.3296).

The ‘free-ratio’ branch model yielded interesting results. In floral buds, the median ω values for female-biased, male-biased, and unbiased genes were 0.222, 0.265 and 0.226, respectively (Fig. 4C and Table S5). Male-biased genes had a significantly higher median relative to female-biased genes (19.37% higher, Wilcoxon rank sum test, P = 0.0009) and unbiased genes (17.26% higher, Wilcoxon rank sum test, P = 0.0004) in floral buds. However, there was no significant difference in ω values between female-biased and unbiased genes (Wilcoxon rank sum test, P = 0.9862). In mature flowers, the median ω values for female-biased, male-biased, and unbiased genes were 0.300, 0.148 and 0.227, respectively (Fig. 4D and Table S5). Female-biased and unbiased genes had significantly higher ω values than male-biased (Wilcoxon rank sum test, P = 0.0101, P = 0.0146, respectively). However, there was no significant difference in ω values between female-biased and unbiased genes (Wilcoxon rank sum test, P = 0.2887). Since the number and evolutionary rates of male-biased genes in mature flowers are lower than those in flora buds, we decided to focus on the latter in subsequent analyses.

Evidences and contributions of positive selection and relaxed selection for male-biased genes in floral buds

After comparing the A model and null model, we discovered that 39 out of 343 OGs (11.34%) exhibited strong evidence of having certain sites that evolved under positive selection based on foreground 2b ω value, likelihood ratio tests (LRTs, P < 0.05) and Bayes empirical Bayes (BEB) value (Fig. 5 and Table S6). As an additional approach, we utilized the aBSREL and BUSTED methods that were implemented in HyPhy v.2.5 software, and detected significant evidence of positive selection. According to our findings, 84 out of 343 OGs (24.49%) were identified to be under episodic positive selection in male-biased genes of floral buds with a site proportion of 0.17%–26.44% based on aBSREL (Table S7). In addition, 69 out of 343 OGs (20.01%) exhibited significant signs of positive selection with the site proportion of 0.28%–32.65% in male-biased genes of floral buds according to BUSTED (Table S8). Among these, a total of 32 OGs (9.30%) were identified through our tests using CodeML, aBSREL and BUSTED (Fig. 5).

Fig. 5

Relaxed selection may occur when the efficiency of natural selection is dramatically reduced or eliminated at both gene-wide and genome-wide levels. This has been proposed as an explaination for the rapid evolution of sex-biased genes (Lahti et al., 2009; Mank, 2017). Using the RELAX model, we conducted an analysis and found that 18 out of 343 OGs (5.23%) showed significant evidence of relaxed selection (K = 0.0184–0.6497) (Table S9). Interestingly, we also observed that 61 out of 343 OGs (17.73%) exhibited significant evidence of intensified positive selection (K = 2.3363–50, ω₂ ≥ 1) (Fig. 5 and Table S10), which is consistent with the results obtained from CodeML, aBSREL and BUSTED.

According to previous studies (Ellegren and Parsch, 2007; Catalan et al., 2018), genes that exhibit sex-biased expression with rapid evolutionary rates tend to display a lower codon bias usage compared to unbiased genes. In our results, we found that male-biased genes in floral buds had a significantly lower median effective number of codons (ENCs) than both female-biased and unbiased genes (Wilcoxon rank sum test, female-biased vs male-biased genes, P = 0.0001 and male-biased vs unbiased genes, P = 0.0123). This suggested that male-biased genes in floral buds exhibit a stronger codon bias usage than both female-biased and unbiased genes (Fig. S3). In summary, our analyses indicated that rapid evolutionary rates of male-biased genes in floral buds were not associated with a reduction in codon usage bias.

To confirm the contributions of positive selection and relaxed selection to rapid rates of male-biased genes in floral buds, we generated three datasets of OGs by excluding different sets of genes. Specifically, we excluded 18 relaxed selective male-biased genes (5.23%), 98 positively selected male-biased genes (28.57%), and 112 male-biased genes (32.65%) under positive and relaxed selection from 343 OGs (Fig. S4). We observed that after excluding relaxed selective male-biased genes, the median (0.264) decreased by 0.34% compared to the median (0.265) of all OGs (Fig. S4A-B). However, after excluding positively selected male-biased genes, the median (0.236) reduced by 11% (Fig. S4A, C) in the results of ‘free-ratio’ branch model. This pattern was consistene in the results of ‘two-ratio’ branch model as well (Fig. S4E-G). Additionally, we analyzed female-biased and unbiased genes that underwent positive and relaxed selection in floral buds (Tables S6-S10). We identified 216 (18.86%) positively selected (Fig. S5A), and 69 (6.03%) relaxed selective female-biased genes from 1,145 OGs (Fig. S5B), respectively. Similarly, we found 436 (18.33%) positively selected (Fig. S6A), and 43 (1.81%) unbiased genes under relaxed selection from 2,378 OGs (Fig. S6B), respectively. Notably, male-biased genes have a higher proportion (10%) of positively selected genes compared to female-biased and unbiased genes. However, relaxed selective male-biased genes have a higher proportion (3.24%) than unbiased genes, but about 0.8% lower than that of female-biased genes. In summary, our analyses suggested that positive selection and relaxed selection likely drove the rapid evolutionary rates of male-biased genes compared to female-biased and unbiased genes in floral buds, but positive selection had an even greater impact.

Functional analysis of sex-biased genes in floral buds

We conducted KEGG pathway enrichment analysis on sex-biased genes in floral buds. Our results showed that 699 genes were female-biased and 358 genes were male-biased, with significant enrichment (P < 0.05) in 26 and 24 KEGG pathways, respectively (Table S11). In the floral bud stage, we observed that female-biased genes were mainly enriched in metabolic and signaling pathways, such as ribosome, Fatty acid elongation, photosynthesis and plant hormone signal transduction (Fig. S7A and Table S11). On the other hand, male-biased genes were significantly enriched in metabolic and signaling pathways, including Inositol phosphate metabolism, starch and sucrose metabolism, regulation of autophagy, plant hormone signal transduction, and Toll-like receptor signaling pathway (Fig. S7B and Table S11).

We have also found that certain male-biased genes, which are evolving under positive selection and relaxed selection (Table S12 and S13), were related to abiotic and biotic stress. For instance, mitogen-activated protein kinase kinase kinase 18 (MAPKKK18) (Zhang and Zhang 2022,), zinc finger CCCH domain-containing protein 20 (C3H20/TZF2) (Bogamuwa and Jang, 2014), and heat stress transcription factor B-3 (HSFB3) (Scharf et al., 2012) have been linked to stress. Interestingly, ten male-biased genes with rapid evolutionary rates were associated with anther and pollen development. These genes include LRR receptor-like serine/threonine protein kinase (LRR-RLK) (Cui et al., 2022), pollen receptor-like kinase 3 (PRK3) (Muschietti and Wengier, 2018), autophagy-related protein 18f (ATG18f) (Zhou et al., 2015; Li et al., 2020a), and plant homeodomain (PHD) finger protein 3 (MALE STERILITY 3) (Hou et al., 2022) in floral buds of male plants.

Discussion

The Cucurbitaceae family is an excellent model for studying the evolution of sexual systems of angiosperms, including sex-determination mechanism and sexual dimorphism. In this study, we aimed to investigate whether sex-biased genes exhibited the evidence of rapid evolutionary rates of protein sequences and identify the evolutionary forces responsible for the observed patterns in the dioecious T. pilosa. We compared the expression profiles of sex-biased gene between sexes and two tissue types and examined the signatures of rapid sequence evolution for sex-biased genes, as well as the contributions of potential evolutionary forces.

Sex-biased expression in floral buds

Numerous studies have shown that in dioecious plants, male-biased genes tend to outnumber female-biased genes. For instance, insect-pollinated dioecious plants such as Asparagus officinalis (Harkess et al., 2015) and Silene latifolia (Zemp et al., 2016), exhibit a higher proportion of male-biased genes, whereas wind-pollinated dioecious plants like Populus balsamifera (Sanderson et al., 2019) has twice as many female-biased genes than male-biased genes. These differences in these studies could be partly attributed to the impact of sexual selection on secondary sexual traits in insect-pollinated dioecious plants, as opposed to wind-pollinated ones (Delph and Herlihy, 2012; Muyle, 2019; Sanderson et al., 2019). In contrast to the above studies, our findings revealed that the number of female-biased genes in floral buds of the night-flowering, insect-pollinated dioecious plant Trichosanthes pilosa exceeded that of male-biased genes by 882. This excess of female-biased expression could be due to lower energy consumption needs and reduced chemical defense capability against herbivores in short-lived male flowers (Sanderson et al., 2019). Indeed, functional enrichment analysis in chemical defense pathways such as terpenoid backbone and diterpenoid biosynthesis, indicated that female floral buds had more expressed genes and were better equipped to defend against herbivores compared to male floral buds (Fig. S7A and Table S11). Additionally, our enrichment analysis showed that the photosynthesis, porphyrin and chlorophyll metabolism pathways were more active in female floral buds compared to male floral buds (Fig. S7A and Table S11), enabling them to acquire more resources such as carbon for fruit and seed production (Delph,1999).

We identified functional enrichments in Toll-like receptor signaling, NF-kappa B signaling, and Inositol phosphate metabolism pathways in male floral buds (Fig. S7B, Tables S11 and S12). These findings suggest that male floral buds are better adaptated to biotic and abiotic environments. Moreover, the enrichment in regulation of autophagy pathways could be associated with the senescence for male floral buds (Table S14) (Liu and Bassham, 2012; Li et al., 2020). We also observed that homologous genes of two male-biased genes (Table S14) that control the raceme inflorescence development (Teo et al., 2014) were highly expressed in floral buds. This suggests that sex-biased genes, rather than sex-specific genes, are responsible for dimorphic traits of inflorescence. In short, these results indicate that changes in sex-biased genes expression play different roles on sexual dimorphism in physiological and morphological traits (Dawson and Geber, 1999).

Rapid Evolution of male-biased genes in floral buds

It has been observed that, in most animals, sex-biased genes, particularly those biased towards males, often exhibit more rapid evolutionary rates than unbiased genes (Parsch and Ellegren, 2013; Grath and Parsch, 2016; Mank, 2017; Toubiana et al., 2021). However, in dioecious angiosperms, no evidence of rapid evolution in sex-biased genes relative to unbiased genes has been found (Zemp et al., 2016; Darolti et al., 2018; Cossard et al., 2019; Sanderson et al., 2019; Scharmann et al., 2021). In contrast, our findings indicated that male-biased genes exhibited rapid evolutionary rates than both female-biased and unbiased genes in floral buds of dioecious T. pilosa. We proposed that positive selection and relaxed purifying selection may be responsible for the rapid sequence evolution of male-biased genes.

After analyzing the data, we found that around 28.57% (98 genes) of male-biased genes have undergone positive selection. Additionally, we observed that the proportion of male-biased genes under positive selection was about 10% higher than that of female-biased and unbiased genes. Furthermore, we discovered that some male-biased genes under positive selection were linked to abiotic and biotic stress responses (Table S12). Our findings are consistent with studies on Drosophila and Ectocarpus (Zhang and Parsch, 2005; Lipinska et al., 2015), suggesting that adaptive evolution is one of the primary evolutionary forces for rapid evolutionary rates. Notably, we identified several male-biased genes under positive selection that are functionally related to early flowering (phyB) (Stephenson and Bertin, 1983; Forrest, 2014; Hajdu et al., 2015) and pollen development (Skogsmyr and Lankinen, 2002; Williams and Reese, 2019) (Tables S12-S14). These findings indicate that a small fraction of male-biased genes may experience adaptive evolution due to sexual selection, driven by male-male competition.

Alternatively, relaxed constraints could contribute to the rapid evolutionary rates of sex-biased genes through three key characteristics (Dapper and Wade, 2020; Tosto et al., 2023). First, sex-biased genes are often expressed solely in reproductive tissues of one sex (e.g., sex-specific genes), particularly in the haploid phase (Immler, 2019; Beaudry et al., 2020). Sex-specific selection (e.g., relaxed purifying selection) acting on sex-specific genes could decrease the elimination of deleterious mutations (Mank, 2017), such as pollen-specific (Harrison et al., 2019) or testes-specific genes (Gershoni and Pietrokovski, 2014). However, our results did not reveal male-specific genes that underwent relaxed purifying selection. Second, sex-biased genes are often expressed in few tissues (tissue-specific genes) (Meisel, 2011; Tosto et al., 2023). These genes could evolve under positive selection or relaxed purifying selection due to low pleiotropy (Congrains et al., 2018; Whittle et al., 2021; Tosto et al., 2023). In our results, 343 male-biased genes with faster evolutionary rates relative to female-biased and unbiased genes overlapped with tissue-biased genes in floral buds (27 out of 343, 7.87%) (Fig. S8A). Furthermore, 27 out of 343 male-biased genes (that is, tissue-biased genes) in floral buds overlapped with nine out of 98 (9.18%) male-biased genes under positive selection (Fig. S8B), and one out of 18 (5.56%) male-biased genes under relaxed purifying selection (Fig. S8C). Consistent with male-biased genes in Fucus (Hatchett et al., 2023), our analyses suggested that elevated evolutionary rates may partly be linked to low pleiotropy. Finally, gene duplication has long been thought to promote functional divergences and phenotypic novelties by relaxing the constraints of purifying selection on the duplicated gene copy early in its history (Lynch and Conery, 2000; Lynch and Katju, 2004; Lahti et al., 2009). For instance, the progesterone receptor gene family in the human lineage (Marinić and Lynch, 2020) and the CYP98A9 clade in Brassicales (Liu et al., 2016) have demonstrated rapid evolution and divergent function due to relaxed purifying selection. In our results, we identified 18 out of 343 (5.25%) male-biased genes that underwent relaxed purifying selection using RELAX model (Table S13). Interestingly, the vast majority of genes under relaxed selection were members of different gene families generated by gene duplication (including whole-genome duplication), such as LOB domain-containing protein 18 (LBD18) (Zhang et al., 2020), WRKY transcription factor 72 (WRKY72) (Chen et al., 2017), and pollen receptor-like kinase 3 (PRK3) (Muschietti and Wengier, 2018).

Additionally, although reducing codon usage bias could theoretically accelerate evolutionary rates of sex-biased genes by decreasing synonymous substitution rates, our results did not support this idea. In fact, we found that male-biased genes had significantly higher median d_S values than female-biased and unbiased genes, both in the ‘free-ratio’ analysis (Fig. S9A, female-biased vs male-biased genes, P = 6.444e-12 and male-biased vs unbiased genes, P = 4.564e-13) and ‘two-ratio’ branch model (Fig. S9B, female-biased vs. male-biased genes, P = 2.2e-16 and male-biased vs. unbiased genes, P = 9.421e-08, respectively), indicating higher synonymous substitution rates.

Several species have been observed to exhibit rapid evolutionary rates of sequences on sex chromosome compared to autosomes, which has been related to the evolutionary theories of fast-X or fast-Z (Meisel and Connallon, 2013; Wright et al., 2015; Charlesworth et al., 2018). Furthermore, the quantification of gene expression by bulk RNA-seq technology, relative to single-cell transcriptome analysis, has been shown to potentially obfuscating true signals in the evolution of sex-biased gene expression in complex aggregations of diverse cell types (Darolti and Mank, 2023; Tosto et al., 2023). However, investigation of these interesting issues related to sex-biased gene evolution in T. pilosa can only be conducted when whole genome sequences become available in the near future.

Methods

Plant materials and RNA isolation

Floral buds (≤ 3 mm) and mature flowers were sampled from three female and three male plants (Fig. 1) in the field of Qinglong Gorge of Anning, Yunnan Province. Floral buds from female and male plants were named F1 and M1, respectively. Similarly, mature flowers from female and male plants, there were named F2 and M2, respectively (Table S1). To exclude possible bacterial contamination, all tissues were sterilized with 75% alcohol and immediately rinsed with purified water. All samples were then snap-frozen in liquid nitrogen, and stored at −80°C. Total RNA was extracted from each sample using TRIzol reagent (Life Technologies, CA, USA) according to the manufacturer’s instructions. The quantification and qualification of RNA were assessed by the RNA Nano 6000 Assay Kit of the Bioanalyzer 2100 system (Agilent Technologies, CA, USA).

Illumina sequencing, de novo assembly and annotation

To construct the library, approximately 2 μg of total RNA was used with the Illumina NEBNext UltraTM RNA Library Prep Kit. RNA sequencing was performed on the Illumina NovaSeq 6000, generating 150 bp paired-end reads. The resulting clean reads were obtained by removing adapters, reads containing N bases and low-quality reads using Trimmomatic v.0.39 (Bolger et al., 2014). These reads were deposited in the NCBI database (PRJNA899312).

De novo assembly for clean reads from all samples was performed using Trinity v.2.10.0 (Haas et al., 2013) with min_kmer_cov: 3 and all other default parameters. To eliminate contamination, all transcripts of de novo assembly were compared to bacterial genomes downloaded from NCBI databases using BLASTN with e-value of 1.0e-05 in blast+ 2.12.0 software. We used Corset v.4.6 (Davidson and Oshlack, 2014) to obtain high quality, non-redundant consensus transcripts (unigenes). TransDecoder v.5.5.0 was run with -m 100 parameters, namely at least 100 amino acids, to predict the coding DNA and protein sequences (Haas et al., 2013).

To evaluate the accuracy and completeness of reference transcriptomes, we performed gene function annotations based on the following databases, using BLAST with a cutoff e-value of 1.0e-05: NR, NT and Swissport (Shiryev et al., 2007). We mapped the unigenes to Pfam database using InterProScan v.5.41 (Jones et al., 2014), to the GO database using Blast2GO (Conesa et al., 2005), and to the KEGG database using KEGG automatic annotation server (Moriya et al., 2007). Additionally, we estimated the completeness of reference transcriptomes using BUSCO v.5.2.2 based on embryophyta_odb9 database (Seppey et al., 2019).

Detection of sex-biased genes

Clean reads were mapped onto all unigenes using Bowtie2 (Langmead and Salzberg, 2012). Read counts were normalized to FPKM (Fragments Per Kilobase Million) value for each unigene using RSEM (Li and Dewey, 2011) in different male and female samples. Genes with zero read counts (i.e., no expression) in both two sexes and tissues were excluded. Differential expression analysis between sexes and tissue types was performed using DESeq2 R package. Unigenes with an FDR-adjusted P < 0.05 and an absolute value of log₂ ratio ≥ 1 identified by DESeq2 were considered as sex-biased genes. To perform KEGG functional enrichment, we use all KEGG annotation terms for all genes as the background and performed the analyses using KOBAS v.2.0.12 (Mao et al., 2005).

Evolutionary rate analyses

To quantify the evolutionary rates of sex-biased genes, we download published genome datasets for monoecious Trichosanthes anguina (Ma et al., 2020) and monoecious Luffa cylindrica which has a closer phylogenetic relationship with Trichosanthes (de Boer et al., 2012; Wu et al., 2020) from CuGenDB database (Zheng et al., 2019). Additionally, we also download published RNA sequencing reads of floral buds and flowers from CNCB (Accession CRA002313) and NCBI databases (Accession SRR5259239) for dioecious plant Trichosanthes kirilowii (Hu et al., 2020), and de novo assembled by previously described methods.

We identified one-to-one OGs using OrthoFinder v.2.3.3 with default parameters from T. anguina, T. pilosa, T. kirilowii, and Luffa cylindrica (Emms and Kelly, 2019). Then, we employed TranslatorX with -c 1 -p M -g -b5 n parameters (i.e., the multiple alignment and the trimming using Muscle and GBlocks, respectively), translated nucleotide sequences and back-translated amino acid alignments into nucleotide alignments to ensure codon-to-codon alignment (Abascal et al., 2010). The remaining gapless alignments (≥100 bp in length) were retained.

To investigate the evolutionary rates of code sequences, we estimated nonsynonymous substitution (d_N), synonymous substitution (d_S) rates, as well as protein substitution rates (d_N/d_S, ω), using two branch models from CodeML package in PAML v.4.9h with the F3x4 codon frequencies (CodonFreq = 2) (Yang, 2007). According to the phylogenetic relationships of Trichosanthes (de Boer et al., 2012; Guo et al., 2020), we set up tree structure ((T. anguina, T. pilosa), T. kirilowii, L. cylindrica) in the control file of CodeML. First, we employed a ‘two-ratio’ branch model (model = 2, Nssites = 0) that assumes the foreground (two dioecious species) has a different ω value from the background (two monoecious species) to estimate and compare the divergences of the foreground. Second, to reduce the potential bias of ω value due to the conflation of two dioecious species, we also implemented a ‘free-ratio’ branch model (model = 1, Nssites = 0), which assumes an independent ω ratio for each branch. Finally, we excluded all OGs with ω > 2, plotted the distribution of ω values and compared the median of ω values in female-biased, male-biased and unbiased orthologous genes of floral buds and mature flowers. All comparisons between sex-biased and unbiased genes were tested using Wilcoxon rank sum test in R software.

Estimation of the strength of natural selection

The rapid evolutionary rates of sex-biased genes may be attributed by positive selection, relaxed selection and lower codon usage bias (Catalan et al., 2018; Dapper and Wade, 2020). Therefore, we conducted separate analyses using classical branch-site models, the adaptive branch-site random effects likelihood (aBSREL) model, the branch-site unrestricted statistical test for episodic diversification (BUSTED) model, the RELAX model, and the effective number of codons (ENC) in PAML v.4.9h, HyPhy v.2.5 and CodonW v.1.4.2 to distinguish which evolutionary forces are driving the rapid evolutionary rates of sex-biased genes (Pond et al., 2020).

To determine if amino acid sites in the foreground, including the T. pilosa lineage have undergone positive selection (foreground 2b ω >1) compared with the background for each OGs, we used branch-site model A (model = 2, Nssite = 2, fix_omega = 0, omega = 1.5) and branch-site model null (model = 2, Nssite = 2, fix_omega = 1, omega = 1). We examined the significance of likelihood ratio tests (LRTs, P < 0.05) to identify positively selected sites between model A and model null by comparing LRTs to the Chi-square distribution with two degrees of freedom. We adjusted the LRTs P value for multiple comparisons using Benjamini and Hochberg’s (FDR) algorithm. When the P value was significant, we used Bayes Empirical Bayes (BEB) estimates to identify sites with a high posterior probability (pp ≥ 0.95) of being under positive selection (Yang et al., 2005; Zhang et al., 2005).

To detect episodic positive selection at a proportion of sites on the foreground branch, we employed the aBSREL method in the HyPhy v.2.5 packages to compare the fully adaptive model (ω > 1) to the null model that allows no positive selection rate classes by LRTs, which is an improved algorithm of branch-site models in PAML (Smith et al., 2015). For relatively small datasets, such as those with fewer than 10 taxa, the aBSREL method may not have enough power to detect positive selection. Therefore, we also ran the BUSTED method to identify gene-wide evidence of episodic positive selection at least one site on at least one branch (Murrell et al., 2015). We set T. pilosa as the foreground and assessed the statistical significance (P < 0.05) using LRTs with the Holm-Bonferroni correction.

To test the relaxation of selective strength, we utilized the RELAX model in the HyPhy v.2.5 software (Wertheim et al., 2015; Schrader et al., 2021). The RELAX model estimates three ω parameters (ω₀ ≤ ω₁ ≤ 1 ≤ ω₂), and determines the proportion of sites in the test (foreground) and reference (background) branches using a branch-site model. The first two ω classifications indicate that sites have undergone purifying selection, and the third classification indicates that sites have been under positive selection. Additionally, the model introduces a selection intensity parameter (K value) to compare a null model (K = 1) with an alternative model, thereby assessing the the strength of natural selection. When K > 1, it suggests intensified natural selection, when K < 1, indicates relaxed natural selection in the test branch relative to the reference branch. We quantified the statistical confidence of K value (P < 0.05) using LRTs and the Holm-Bonferroni correction.

To investigate codon usage bias, which refers to the differences in the frequency of occurrence of synonymous codons in coding DNA, we employed CodonW v.1.4.2. This program considers the ENC values from 20 to 61 as a measure of the departure of the genetic codes for a given gene (Wright, 1990), with lower ENC values represent stronger codon usage bias (Hambuch and Parsch, 2005). We performed a Wilcoxon rank sum test to determine if there were deviations in ENC values among female-biased, male-biased, and unbiased genes in floral buds.

Acknowledgements

We are indebated to Ting Zhang, Zhi-Yun Yang, Jiang-Li Ma, Peng-Fei Ma, Xu-Kun Wu, Shi-Yu Lv and other members of staff of Germplasm Bank of Wild Species for sampling. We are very grateful to Spencer C. H. Barrett (University of Toronto, Canada) for his critical reading and suggestions on the manuscript. We also thank to the iFlora HPC Center (iFlora High Performance Computing Center) of Germplasm Bank of Wild Species for computationoal support on data analysis. This study was funded by the Strategic Priority Research Program of the Chinese Academy of Sciences, China (XDB31000000), the Key R and D Program of Yunnan Province, China (202103AC100003), the Key Basic Research Program of Yunnan Province, China (202101BC070003), the Science and Technology Basic Resources Investigation Program of China (2019FY100900), the National Natural Science Foundation of China (31570333), and the open research project of “Cross-Cooperative Team” of the CAS’ Germplasm Bank of Wild Species, Kunming Institute of Botany. The study was also supported by National Wild Plant Germplasm Resource Center.

Author contributions

HTL, DZL and LZ conceived the project the study. HTL and LZ designed and performed the experiment. LZ, JH and HTL collected the samples. LZ analyzed the data and wrote the manuscript. HTL, DZL, LZ, and WZ revised the manuscript. All authors approved the final manuscript.

Additional files

Supplementary files

Supplementary file 1. Supplementary Figures S1 to S9

Supplementary Figure S1. The length distribution of unigenes and the results of alignment by BLAST.

Supplementary Figure S2. Boxplot of d_N/d_S values of all female-biased, male-biased and unbiased genes (including ω ≥ 2) in floral buds and mature flowers of the dioecious T. pilosa. Significant differences are represented by * in Wilcoxon rank sum tests (***P < 0.0005, **P < 0.005 and *P < 0.05).

Supplementary Figure S3. Violin plots of ENCs values of female-biased, male-biased and unbiased genes in floral buds. Significant differences using Wilcoxon rank sum tests are represented by * (***P < 0.0005 and *P < 0.05).

Supplementary Figure S4. Boxplot of d_N/d_S values of male-biased genes, including all 343 OGs, excluding 18 OGs under relaxed selection, excluding 98 OGs under positive selection, and excluding 112 OGs under both positive and relaxed selection, and boxplot of d_N/d_S values of all female-biased and unbiased genes in floral buds.

Supplementary Figure S5. Venn diagrams of female-biased genes under positive selection and relaxed selection in floral buds. Overlaps of female-biased genes were detected to be under positive selection using aBSREL, BUSTED, CodeML and RELAX, and to be under relaxed selection using RELAX.

Supplementary Figure S6. Venn diagrams of unbiased genes under positive selection and relaxed selection in floral buds. Overlaps of unbiased genes were identified to be under positive selection using aBSREL, BUSTED, CodeML and RELAX, and to be under relaxed selection using RELAX.

Supplementary Figure S7. Scatterplots of KEGG pathway of sex-biased in floral buds of the dioecious T. pilosa.

Supplementary Figure S8. The overlap between male-biased genes with faster evolutionary rates and tissue-biased genes in floral buds.

Supplementary Figure S9. Boxplot of d_S values of female-biased, male-biased and unbiased genes in floral buds of dioecious T. pilosa. Significant differences are represented by * in Wilcoxon rank sum tests (***P < 0.0005).

Supplementary file 2. Supplementary Tables S1 to S14

Supplementary Table S1. Overview of sequencing reads from 12 samples of male and female plants in T. pilosa.

Supplementary Table S2. Numbers of unigenes annotated in public databases.

Supplementary Table S3. The mapping rate of reads for each sample in floral buds and mature flowers of the dioecious T. pilosa.

Supplementary Table S4. d_N, d_S and ω values of each female-biased, male-biased, unbiased orthologous genes of floral buds and mature flowers for each species using ‘two-ratio’ branch model of CodeML in PAML.

Supplementary Table S5. d_N, d_S and ω values of each female-biased, male-biased, unbiased orthologous genes of floral buds and mature flowers for each species using ‘free-ratio’ branch model of CodeML in PAML.

Supplementary Table S6. Genes under positive selection identified by branch-site model of CodeML in PAML and functions in NR, KEGG, Swissport and GO databases for male-biased orthologous genes in floral buds.

Supplementary Table S7. Genes under episodic positive selection tested by aBSREL model in HyPhy and functions in NR, KEGG, Swissport and GO databases for male-biased orthologous genes in floral buds.

Supplementary Table S8. Genes under episodic positive selection found by BUSTED model in HyPhy and functions in NR, KEGG, Swissport and GO databases for male-biased orthologous genes in floral buds.

Supplementary Table S9. Genes under relaxed selection detected by RELAX model in HyPhy and functions in NR, KEGG, Swissport and GO databases for male-biased orthologous genes in floral buds.

Supplementary Table S10. Genes under intensified positive selection identified by RELAX model in HyPhy and functions in NR, KEGG, Swissport and GO databases for male-biased orthologous genes in floral buds.

Supplementary Table S11. KEGG pathway enrichment analysis of female-biased and male-biased genes in floral buds of the dioecious T. pilosa.

Supplementary Table S12. Functions and references associated with biotic and abiotic stress responses, organ developments of male-biased genes under significant positive selection (P < 0.05) in floral buds.

Supplementary Table S13. Functions and references associated with biotic and abiotic stress responses, organ developments of male-biased genes under significant relaxed selection (P < 0.05) in floral buds.

Supplementary Table S14. The expressions and functions of some male-biased genes associated with senescence, raceme inflorescence development and early flowering in floral buds.

Data Availability

All RNA-Sequencing clean reads have been deposited in the databases of the National Center for Biotechnology Information (NCBI) under BioProject ID PRJNA899312.