Introduction

Sexual dimorphism is the condition where sexes of the same species exhibit 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 between the two sexes, i.e., sex-biased genes (including sex-specific genes expressed 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; Grath and Parsch, 2016; Charlesworth, 2018; Tosto et al., 2023). 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; Naqvi et al., 2019; Cai et al., 2021; Mank, 2023; Murat et al., 2023; Singh and Agrawal, 2023). Previous studies revealed that sex-biased gene expressions were associated with the evolution of sexual dimorphisms in some animal species, although the extent of this bias exhibits great variation among taxa, tissues, and development stages (Mank, 2017; Hsu et al., 2020; Khodursky et al., 2020; Lichilin et al., 2021; Toubiana et al., 2021; Djordjevic et al., 2022; Yue et al., 2023). Unlike most animals, the vast majority (~90%) of flowering plants (angiosperms) are hermaphroditic, while only a small fraction (~5%) are dioecious in which individuals have exclusively male or female reproductive organs (Renner, 2014). Most dioecious plants possess homomorphic sex-chromosomes that are roughly similar in size when viewed by light microscopy (Palmer et al., 2019). Furthermore, sexual dimorphism in dioecious plants is less common and less conspicuous than in most animals (Barrett and Hough, 2013). Hence, the study of sex-biased gene expression is of great interest to plant evolutionary biologists, as it is necessary to understand the evolution of sexual dimorphism in dioecious plants (Moore and Pannell, 2011).

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 (Prö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 hermaphroditic Arabidopsis thaliana (Gossmann et al., 2014; 2016), Solanum (Moyle et al., 2021), and dioecious 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 (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 detected between male and female plants leaf tissues, 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 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 the five above-mentioned studies have investigated expression differences and protein evolutionary rates of sex-biased genes in dioecious angiosperms. Moreover, these studies only compared gene expression in vegetative versus vegetative tissues and vegetative versus reproductive tissues, limiting our understanding of sexual selection at different floral development stages. 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 versus reproductive tissues, e.g., different floral development stages in dioecious angiosperms.

Evolutionary analyses indicate that different driving 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 Drosophila, the elevated sequence divergence rates of male-biased genes have often been interpreted as the signature of adaptive evolution, suggesting that sexual selection is the primary evolutionary force (Pröschel et al., 2006; Assis et al., 2012). In brown algae, female-biased and/or male-biased genes exhibited higher evolutionary rates than unbiased genes, suggesting that rapid evolution is partly driven by adaptive evolution or sexual selection (Lipinska et al., 2015; Cossard et al., 2022; Hatchett et al., 2023). However, studies in plants have never reported elevated rates of sex-biased genes.

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 the others 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 (synonym: T. ovigera, 2n = 22, Cucurbitaceae) is mainly distributed from Southwest and Southeast China to Japan, extending to Southeast Asia, New Guinea and Western Australia. It was suggested to have originated in the late Miocene (ca. 8-6 million-year ago) (de Boer et al., 2012; 2015; Guo et al., 2020). Trichosanthes pilosa is a perennial, night-flowering, insect-pollinated dioecious vine that reproduces sexually and possesses a pair of heteromorphic sex chromosomes XX/XY (Ming et al., 2011). The male parts (e. g., anthers) of female flowers, and the female parts (e. g., pistil and ovaries) of male flowers are fully aborted. However, its male and female flowers exhibit strong sexual dimorphism in floral morphological and phenological traits, such as racemose versus solitary (Fig. 1), early-flowering versus late-flowering, and caducous versus long-lived (Wu et al., 2011).

Fig. 1

To understand the evolution of sex-biased genes in dioecious T. pilosa, we collected floral buds and flowers at anthesis from male and female individuals and characterized their expression profiles using Illumina RNA sequencing. Our primary objectives are to 1) compare expression divergence between males and females at two floral development 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 floral buds and flowers at anthesis from female and male of dioecious T. pilosa. We set up three biological replicates from three female and three male plants, including 12 samples in total (six floral buds and six flowers at anthesis). We then generated a total of nearly 276 million clean reads (Table S1). Due to the absence of a reference genome, we 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 used BUSCO assessments based on embryophyta_odb10 database, which showed the completeness of the reference transcriptome at 89.7% (Table S2). We then annotated them against protein databases including 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 species of Cucurbitaceae (61.6%, 36,375), such as Momordica charantia (16.3%, 9,625), Cucumis melo (11.9%, 7,027), Cucurbita pepo (11.9%, 7,027), Cucurbita moschata (11.5%, 6,791), Cucurbita maxima (10.1%, 5,964), and other species (38.4%, 22,676) (Fig. S1C). Overall, our assessment 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 flowers at anthesis 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 flowers at anthesis (Fig. 2B). Using hierarchical clustering analysis, we evaluated different levels of gene expression across sexes and tissues (Fig. 2C), revealing that female floral buds maybe tend to masculinize, and vice versa in the early stages of floral development. Furthermore, we observed that the number of sex-biased genes in floral buds was approximately 15 times higher than in flowers at anthesis, indicating that sex-biased genes associated with meiotic processes, 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 flowers at anthesis 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 flowers at anthesis, we only identified 26 out of 380 (6.84%) female-specific genes and 52 out of 233 (22.32%) male-specific genes. Taken together, sex bias is more prevalent in floral buds than in flowers at anthesis.

Fig. 2

Tissue-biased/stage-biased gene expression

We compared the expression levels of transcripts in floral buds and flowers at anthesis within each sex to identify genes with tissue-biased expression. In male plants, the number (M2TGs: n = 2,795) of tissue-biased genes in male flowers at anthesis (M2TGs) was 1,040 higher than that in male floral buds (M1TGs: n = 1,755, Fig. 3A and 3B). However, in female plants, the number (F2TGs: n = 660) of tissue-biased genes in female flowers at anthesis (F2TGs) was only 536 more than that in female floral buds (F1TGs: n = 124, 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 female-biased genes in floral buds (F1BGs: n = 5,096) overlapped with tissue-biased genes in female floral buds (F1TGs: n = 124) and female flowers at anthesis (F2TGs: n = 660), accounting for only 85 out of 5,096 (1.67%) (Fig. 3C). Similarly, few female-biased genes in flowers at anthesis (F2BGs: n = 380) overlapped with tissue-biased genes in female floral buds (F1TGs: n = 124) and female flowers at anthesis (F2TGs: n = 660), occupying around 5 out of 380 (1.32%) (Fig. 3D). However, a significant proportion of male-biased genes in floral buds (M1BGs: n = 4,214) overlapped with tissue-biased genes in male floral buds (M1TGs: n = 1,755) and male flowers at anthesis (M2TGs: n = 2,795), with 1,010 out of 4,214 (23.97%) (Fig. 3A). A high proportion of male-biased genes in flowers at anthesis (M2BGs: n = 233) overlapped with tissue-biased genes in male floral buds (M1TGs: n = 1,755) and male flowers at anthesis (M2TGs: n = 2,795), 145 out of 233 (62.23%) (Fig. 3B).

Fig. 3

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

We compared rates of protein evolution among male-biased, female-biased and unbiased genes in four species with phylogenetic relationships (((T. anguina, T. pilosa), T. kirilowii), Luffa cylindrica), including dioecious T. pilosa, dioecious T. kirilowii, monoecious T. anguina in Trichosanthes, together with monoecious Luffa cylindrica. To do this, we sequenced transcriptomes of T. pilosa. We also collected transcriptomes of T. kirilowii, as well as genomes of T. anguina and Luffa cylindrica (see Methods Section). We identified 1,145 female-biased, 343 male-biased, and 2,378 unbiased one-to-one orthologous groups (OGs) from floral buds. 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-ratio’ and ‘free-ratio’ branch models in juvenile and mature flowers (Figs. 4 and S2).

Fig. 4

The two-ratio branch model, where the foreground (dioecious branches) has a different ω value relative to the background (all other branches), is better supported than the fixed-ratio branch model, where all branches are constrained to have the same ω value. In the results of the ‘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 flowers at anthesis, 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, 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) possibly because of limited statistical power due to the low number of sex-biased genes in flowers at anthesis.

In free-ratios model, ω values are free to vary in each branch compared to fixed-ratio branch model and two-ratio branch model. 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 flowers at anthesis, 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 of male-biased genes and evolutionary rates of male-biased genes in flowers at anthesis are lower than those in floral buds, we decided to focus on the latter in subsequent analyses. Additionally, we also ran permutation t test, which is almost consistent with the results of Wilcoxon rank sum test. We found that only in floral buds, there were significant differences in ω values in the results of ‘free-ratio’ model (female-biased versus male-biased genes, P = 0.04282 and male-biased versus unbiased genes, P = 0.01114) and ‘two-ratio’ model (female-biased versus male-biased genes, P = 0.01992 and male-biased versus unbiased genes, P = 0.02127, respectively).

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

After comparing classical branch-site model A and branch-site model null models (see Methods Section), we discovered that 39 out of 343 OGs (11.34%) in male-biased genes of floral buds exhibited strong evidence of having certain sites that evolved under positive selection based on foreground ω value, likelihood ratio tests (LRTs, P < 0.05) and Bayes empirical Bayes (BEB) value (Fig. 5 and Table S6). As a complementary approach, we utilized the aBSREL and BUSTED methods that are implemented in HyPhy v.2.5 software, which avoids false positive results by classical branch-site models due to the presence of rate variation in background branches, 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 (e.g., the reduction of the strength of purifying selection) is reduced, leading to accumulations of deleterious mutations (Lahti et al., 2009; Dapper and Wade, 2020). This has been proposed as an explanation for the rapid evolution of sex-biased genes (Lahti et al., 2009; Mank, 2017). Using the RELAX model, we conducted an analysis and detected that 18 out of 343 OGs (5.23%) showed significant evidence of relaxed selection (K = 0.0184–0.6497) (Tables S9). We also found that most of them were members of different gene families generated by gene duplication (Table S13). Additionally, we 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 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 stronger codon bias than both female-biased and unbiased genes (Fig. S3). Similarly, given that d_N/d_S values of sex-biased genes were higher due to codon usage bias, lower d_S rates would be expected in sex-biased genes relative to unbiased genes (Ellegren & Parsch, 2007; Parvathy et al., 2022). However, we exhibited that the median of d_S values in male-biased genes were much higher than those in female-biased and unbiased genes in the results of ‘free-ratio’ (Fig. S4A, female-biased versus male-biased genes, P = 6.444e-12 and male-biased versus unbiased genes, P = 4.564e-13) and ‘two-ratio’ model (Fig. S4B, female-biased versus male-biased genes, P = 2.2e-16 and male-biased versus unbiased genes, P = 9.421e-08, respectively). In short, our analyses indicated that rapid evolutionary rates of male-biased genes in floral buds were not associated with a reduction in codon usage bias.

We also 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. S5), and 69 (6.03%) relaxed selective female-biased genes from 1,145 OGs, respectively. Similarly, we found 436 (18.33%) positively selected (Fig. S6), and 43 (1.81%) unbiased genes under relaxed selection from 2,378 OGs, 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.

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 stress and immune responses. 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. Additionally, 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., 2020), and plant homeodomain (PHD) finger protein 3 (MALE STERILITY 3) (Hou et al., 2022) in floral buds of male plants.

Discussion

The Cucurbitaceae family, where half of the species are monoecious and half are dioecious, is an excellent model for studying the evolution of sexual systems of angiosperms, including sex-determination mechanism and sexual dimorphism (Schaefer and Renner, 2010; Boualem et al., 2015; Ma and Pannell, 2016). In this study, we aimed to investigate whether sex-biased genes exhibited evidence of rapid evolutionary rates of protein sequences and identify the evolutionary forces responsible for the observed patterns in the dioecious Trichosanthes pilosa. We compared the expression profiles of sex-biased genes 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

Several studies have shown that in dioecious plants, male-biased genes tend to outnumber female-biased genes, consistent with the patterns in most animals (Zhang et al., 2007; Djordjevic et al., 2022). 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. In contrast, the wind-pollinated dioecious plant Populus balsamifera (Sanderson et al., 2019) has twice as many female-biased genes as male-biased genes. The 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). Similar to the above study of Populus balsamifera, 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 (~21%). This excess of female-biased expression could be due to lower energy consumption needs and reduced chemical defence capability against insect herbivores in short-lived male flowers (Sanderson et al., 2019). Indeed, functional enrichment analysis in chemical pathways such as terpenoid backbone and diterpenoid biosynthesis indicated that relative to male floral buds, female floral buds had more expressed genes that were equipped to defend against herbivorous insects and pathogens, except for growth and development (Vaughan et al., 2013; Ren et al., 2022) (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, Table S11). We also found that male-biased genes with high evolutionary rates in male floral buds were associated with functions to abiotic stresses and immune responses (Tables S12 and S13), which suggest that male floral buds through rapidly evolving genes are adapted to mountain climate and the environment in Southwest China relative to female floral buds through high gene expression. In addition, the enrichment in regulation of autophagy pathways could be associated with gamete development and the senescence of male floral buds (Table S14) (Liu and Bassham, 2012; Li et al., 2020; Zhou et al., 2021). In fact, it was observed that male flowers senesced faster (Wu et al., 2011). We also found that homologous genes of two male-biased genes in floral buds (Table S14) that control the raceme inflorescence development (Teo et al., 2014) were highly expressed compared to female floral buds. Taken together, these results indicate that expression changes in sex-biased genes, rather than sex-specific genes play different roles in sexual dimorphic traits in physiology and morphology (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 experience higher 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 stress and immune 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 important driving 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 (Sandler et al., 2018; 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, we observed male-biased genes but not male-specific genes undergoing relaxed purifying selection. Second, sex-biased genes are often expressed in few tissues (tissue-biased genes) (Meisel, 2011; Tosto et al., 2023), resulting in these genes rapidly evolving under positive selection or relaxed purifying selection due to low evolutionary constraints (Congrains et al., 2018; Whittle et al., 2021; Tosto et al., 2023). In our results, 343 male-biased genes (M1-biased genes, M1BGs) with faster evolutionary rates relative to female-biased and unbiased genes overlapped with 1,755 tissue-biased genes in floral buds (M1-tissue-biased genes, M1TGs) (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). So, we obtained ten rapidly evolving tissue-biased genes which were also male-biased in male flower buds, suggesting that elevated evolutionary rates may partly be linked to low constraints, consistent with male-biased genes in Anastrepha and Fucus (Congrains et al., 2018; Hatchett et al., 2023). 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, reducing codon usage bias could theoretically accelerate evolutionary rates of sex-biased genes by decreasing synonymous substitution rates. However, our results did not support this idea due to stronger codon usage bias in male-biased genes (Fig. S3). Codon usage bias is influenced by many factors, such as levels of gene expression. Highly expressed genes have a stronger codon usage bias and could be encoded by optimal codons for more efficient translation (Frumkin et al., 2018; Parvathy et al., 2022), consistent with high levels of gene expression in males (that is, male-biased genes) in floral buds. Additionally, stronger codon usage bias may be related to higher synonymous substitution rates (Parvathy et al., 2022). Indeed, male-biased genes had significantly higher median d_S values than female-biased and unbiased genes, both in the ‘free-ratio’ analysis (Fig. S4A) and ‘two-ratio’ branch model (Fig. S4B).

Several species have been observed to exhibit rapid evolutionary rates of sequences on sex chromosomes compared to autosomes, which has been related to the evolutionary theories of fast-X or fast-Z (Meisel and Connallon, 2013; Hough et al., 2014; Wright et al., 2015; Charlesworth et al., 2018; Darolti et al., 2023). Furthermore, the quantification of gene expression by bulk RNA-seq technology, relative to single-cell transcriptome analysis, has been shown to potentially obfuscate 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). Additionally, our samples were relatively small, and may provide low power to detect differential expression and evolutionary analysis. Therefore, investigation of these interesting issues related to sex-biased gene evolution in T. pilosa can only be conducted when whole genome sequences and population datasets become available in the near future.

Methods

Plant materials and RNA isolation

Floral buds (≤ 3 mm) and flowers at anthesis were sampled from three female and three male plants (Fig. 1) from the mountainous regions of Anning (Qinglong Gorge), Yunnan Province in Southwest China. Floral buds from female and male plants were named F1 and M1, respectively. Similarly, flowers at anthesis 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.4.5 based on embryophyta_odb10 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 (Love et al., 2014). 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 used 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 coding 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 structures ((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, to avoid the effects of saturation substitution, we used separately OGs with 0 < ω < 2 and all OGs with ω > 0, plotted the distribution of ω values and compared the median of ω values in female-biased, male-biased and unbiased orthologous genes of floral buds and flowers at anthesis. All comparisons between sex-biased and unbiased genes were tested using Wilcoxon rank sum test in R software. Additionally, we also performed permutations t-tests with 100,000 permutations in the R package Deducer (Fellows, 2012).

Estimation of the strength of natural selection

The rapid evolutionary rates of sex-biased genes may be attributed to 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 that assume different ω values both among branches and across sites (Álvarez-Carretero et al., 2023), the adaptive branch-site random effects likelihood (aBSREL) model (Smith et al., 2015), the branch-site unrestricted statistical test for episodic diversification (BUSTED) model (Murrell et al., 2015), the RELAX model (Wertheim et al., 2015), and the effective number of codons (ENC) in PAML v.4.9h (Yang, 2007), HyPhy v.2.5 (Pond et al., 2020) and CodonW v.1.4.2 (http://codonw.sourceforge.net) to distinguish which evolutionary forces are driving the rapid evolutionary rates of sex-biased genes.

To determine if amino acid sites in the foreground, including the T. pilosa lineage have undergone positive selection (foreground ω >1) compared with the background for each OGs, we followed the steps of Zhang et al. (2005), and 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). The classical branch-site model assumes four site classes (0, 1, 2a, 2b), with different ω values for the foreground and background branches. In site classes 2a and 2b, the foreground branch undergoes positive selection when there is ω > 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).

To detect episodic positive selection at a proportion of sites on the foreground branch, we employed the aBSREL method (Smith et al., 2015) 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. 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 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 very grateful to Spencer C. H. Barrett (University of Toronto, Canada) for his critical reading and suggestions on the manuscript. We also thank three anonymous reviewers for their comments on improving the quality of the manuscript. We are indebted to Ting Zhang, Zhi-Yun Yang, Jiang-Li Ma, Peng-Fei Ma, Xu-Kun Wu, Shi-Yu Lv, Zhen Peng and other members of staff of Germplasm Bank of Wild Species for sampling. We also thank the iFlora HPC Center (iFlora High Performance Computing Center) of Germplasm Bank of Wild Species for computational support on data analysis.

Funding

This study was funded by the Strategic Priority Research Program of the Chinese Academy of Sciences, China (XDB31000000), the National Natural Science Foundation of China (32370233, 31570333), 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), 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 the National Wild Plant Germplasm Resource Center and the CAS Key Technology Talent Program.

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 S8

Supplementary Figure S1. The length distribution of unigenes (A) and the e-value distribution (B) and the species distribution (C) of BLAST hits for each unigene.

Supplementary Figure S2. Boxplot of d_N/d_S values (including all ω values) of female-biased, male-biased and unbiased genes in floral buds and flowers at anthesis of Trichosanthes pilosa. White dot indicates the median of dN/dS values for sex-biased and unbiased genes. Wilcoxon rank sum tests are used to test for significant differences (***P < 0.0005, **P < 0.005 and *P < 0.05). The distributions of d_N/d_S values for female-biased, male-biased and unbiased genes in floral buds (A) and flowers at anthesis (B) using ‘free-ratio’ branch model. The distributions of d_N/d_S values for female-biased, male-biased and unbiased genes in floral buds (C) and flowers at anthesis (D) using ‘two-ratio’ branch model.

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_S values of female-biased, male-biased and unbiased genes in floral buds of dioecious Trichosanthes pilosa using ‘free-ratio’ (A) and ‘two-ratio’ (B) branch model. Significant differences are represented by * in Wilcoxon rank sum tests (***P < 0.0005).

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

Supplementary Figure S7. Scatterplots of KEGG pathway of sex-biased genes in female (A) and male (B) floral buds of the dioecious Trichosanthes pilosa. The y axis represents the pathway name, and the x axis represents the ratio of genes corresponding to the pathways.

Supplementary Figure S8. The overlap between male-biased genes with faster evolutionary rates and tissue-biased genes in floral buds. M1BGs indicate male-biased genes with faster evolutionary rates in male floral buds (M1-biased genes). M1TGs indicate tissue-biased genes in male floral buds (M1-tissue-biased genes). 343 M1BGs overlapped with 1755 M1TGs (A). 27 out of 343 M1BGs&M1TGs overlapped with M1BGs_98 (B). 27 out of 343 M1BGs&M1TGs overlapped with M1BGs_18 (C).

Supplementary file 2. Supplementary Tables S1 to S14

Supplementary Table S1. Overview of sequencing reads from 12 samples of male and female plants in Trichosanthes 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 flowers at anthesis of Trichosanthes pilosa.

Supplementary Table S4. d_N, d_S and ω values of each female-biased, male-biased, unbiased orthologous genes of floral buds and flowers at anthesis for each species using ‘two-ratio’ branch model of CodeML in PAML. Two dioecious species Trichosanthes pilosa and T. kirilowii represent the foreground and two monoecious species T. anguina and Luffa cylindrica represent the background.

Supplementary Table S5. d_N, d_S and ω values of each female-biased, male-biased, unbiased orthologous genes of floral buds and flowers at anthesis 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 function in NR, KEGG, Swissport and GO databases for male-biased orthologous genes in floral buds. Two dioecious species Trichosanthes pilosa and T. kirilowii represent the foreground and two monoecious species, T. anguina and Luffa cylindrica represent the background.

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. The dioecious species Trichosanthes pilosa represents the foreground and other species represent the background.

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. The dioecious species Trichosanthes pilosa represents the foreground (unconstrained branch) and other species represent the background (constrained branch).

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. The dioecious species Trichosanthes pilosa represents the foreground (test) and other species represent the background (reference).

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. The dioecious species Trichosanthes pilosa represents the foreground (test) and other species represent the background (reference).

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

Supplementary Table S12. Functions and references associated with abiotic stress and immune 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 abiotic stress and immune 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. The reference transcriptome, orthology data, and alignments have been uploaded to ResearchGate (www.researchgate.net/publication/373194650_Trichosanthes_pilosa_datasets).