To define the male and female haplotypes of the sex-determining locus in Gracilaria species, we leveraged male reference genomes along with complementary female genomic and transcriptomic data published recently in [39] (Additional file 1: Table S1) and also deposited in the Rhodoexplorer database. Our approach involved the generation of female genome assemblies for all four Gracilaria species (see Additional file 1: Table S2 for assembly metrics) and the application of various bioinformatic strategies combined with experimental validation, as outlined in the Materials and Methods section. Notably, the male genomes of G. vermiculophylla, G. chilensis, and G. gracilis as well as the female genome of G. gracilis represent continuous and high-quality assemblies, whereas the female genomes of G. caudata, G. chilensis, G. vermiculophylla, along with the male genome of G. caudata, are draft assemblies. We assume that the combination of approaches used here has provided a near-exhaustive list of male and female sex-linked genes and genomic regions in all species (Additional file 1: Table S3). It is nevertheless possible that some scaffolds, particularly those that are highly repetitive, may have been missed in the more fragmented assemblies of G. caudata male and female genomes or G. vermiculophylla and G. chilensis female genomes (Additional file 1: Table S2) [39]. Therefore, for the sex chromosome structural analysis, including gene density, repeat content, and GC content, we focused only on the most continuous assemblies of male G. vermiculophylla and G. chilensis along with both the male and female G. gracilis. All other analyses concerning sex-linked genes (gene structure, function, expression, and evolution) were conducted across all species and both sexes.
Sex chromosome architecture in Gracilariales
Similar to other organisms with haploid-diploid life cycles, such as those from green algae, bryophyte, and brown algal lineages [9], Gracilaria species have haploid sex determination with UV sex chromosomes. In these species, the V chromosome carries the male sex-determining region (V-SDR), while the corresponding U-SDR determines sex in females. To perform a fine-scale characterization of the sex chromosome in the Gracilariales, we used the male chromosome-level assembly of G. vermiculophylla [39, 40]. Gracilaria vermiculophylla has 24 chromosomes, and 95% of the assembled sequences could be placed in the 24 largest scaffolds; the male-specific V-SDR, identified using a combination of coverage analysis, k-mer profiling, and PCR validation (see Methods: Identification of non-recombining sex-specific regions), is located on scaffold 41 and contains 28 genes (Fig. 2A, Table 1). Consequently, the V sex chromosome of G. vermiculophylla has a total size of 5.14 Mbp, featuring a small (0.873 Mbp) central sex-determining region (SDR) bordered by two pseudoautosomal regions (PARs) (Fig. 2B). Similarly to G. vermiculophylla, the SDRs in the more continuous assemblies of G. gracilis male, G. gracilis female, and G. chilensis male ranged from 0.645 to 0.920 Mbp (Table 1; Additional file 1: Table S3). It appears, therefore, that the size of the sex locus in all examined Gracilaria species is relatively small and constitutes less than 2% of the total genome sequence [39, 40]. Importantly, we validated previously identified sex markers [29,30,31] using bulk-segregant analysis and hypervariable genetic markers such as Random Amplified Polymorphic DNA (RAPD), Sequence Characterized Amplified Region (SCAR), and SNPs. These rapid methods proved highly effective in detecting sex markers, even when the SDRs are small.
Sex chromosome in G. vermiculophylla. A Genome plot showing GC content, repeat density, and gene density on different chromosomes. Male sex-determining region (V-SDR) on Scaffold_41 (male V sex chromosome) is highlighted by the blue box. B Scatterplot of k-mer analysis (KQ, YGS method) and coverage of female reads over 1-kb-long windows along the Scaffold_41 (V-chromosome) in G. vermiculophylla. The region spanning from 2108 to 2981 kb, displayed no female read coverage (red), high AWK values (green) and high proportion of unmatched k-mers in the female genome (blue, from KQ and YGS analysis respectively), consistent with a non-recombining sex-determining region (V-SDR) of the V-chromosome. AWK values: amount of male-specific 15 bp k-mers (i.e., k-mers with a null ratio between male and female reads) per Kb sequence, KQ: K-mer quotient, defined as the male-to-female frequency ratio of shared k-mers, YGS: Y chromosome Genome Scan, which calculates the percentage of unmatched k-mers when mapping male k-mers to the female genome—used here analogously to identify male-specific regions in UV systems. C Structural characteristics (gene density, repeat content, GC content) of pseudoautosomal regions (PAR) compared to autosomes. Stars above the boxplots indicate significant differences (***p-value < 0.002, permutation tests, 10 k permutations)
Next, we investigated if the UV system in red algae also exhibits the structural characteristics typically associated with sex-specific, non-recombining regions, such as lower gene content and accumulation of repeats compared to autosomes [5, 14, 41, 42]. These patterns are attributed to the suppression of recombination across the SDR which can lead to genetic degeneration unless there is a strong selection on gene function to counteract this effect [43,44,45]. As expected, we found that the gene density on the sex-linked scaffolds in Gracilaria species was significantly lower (p-value < 0.05, permutation test of the difference in the mean, 10 k permutations) and the repeat content was considerably higher (p-value < 0.05, permutation test of the difference in the mean, 10 k permutations) compared to the autosomes (Table 1; Additional file 1: Table S4, S5). However, other indicators of genetic degeneration, such as lower GC content [4, 46] or shorter CDS length [14], did not show notable differences between sex-linked and non-sex-linked regions (p-value > 0.05, permutation test of the difference in the mean, 10 k permutations, Table 1; Additional file 1: Table S4). Furthermore, coding sequences of sex-linked genes in all four species did not display a significant under-representation of optimal codons, typically associated with degeneration in non-recombining regions (e.g., [47]) (Mann–Whitney rank test, p-value > 0.05, Additional file 2: Fig. S1, Additional file 1: Table S5). Therefore, although typical structural characteristics of non-recombining regions, such as low gene content and higher repeat content, are present in the Gracilaria SDR scaffolds, we found little evidence of degeneration in the genes located within these regions. In fact, the SDR genes did not present any distinguishing feature from autosomal genes (Table 1, Additional file 1: Table S5).
Previous studies in brown algal UV systems have demonstrated that pseudoautosomal regions (PARs) exhibit unique characteristics, such as the accumulation of taxonomically restricted, evolutionarily “young” genes and may act as cradles for de novo gene birth [8, 48]. We therefore examined the PARs of the G. vermiculophylla sex chromosome. Specifically, we conducted a phylostratigraphic analysis [49] to estimate the evolutionary age of each gene in the G. vermiculophylla genome. Genes were classified into five phylostrata, ranging from the most ancient (shared with cellular organisms) to the most recent or “young” (restricted to the Gracilariaceae family) (Additional file 1: Table S6, see Methods section: “Phylostratigraphy Analysis”).
We found that the average gene density in the PAR region of G. vermiculophylla was significantly higher than that of autosomes (p-value = 0.0015, permutation tests of the difference in the mean, 10 k permutations; Fig. 2C). However, we found no evidence of accumulation of evolutionary “young” genes (Additional file 2: Fig. S2A, Additional file 1: Table S7). The PAR was instead slightly depleted in young genes compared to autosomes (Χ2 = 16.514, p = 0.035), which stands in contrast to the patterns observed in brown algal UV sex chromosomes [8].
In addition, unlike other UV systems, PARs in G. vermiculophylla exhibited lower repeat content (p-value = 5 × 10−4, permutation tests of the difference in the mean, 10 k permutations) (Fig. 2C) and similar substitution rates (KS-PAR = 0.00590 vs KS-auto = 0.00678; p-value = 0.7567, permutation tests of the difference in the mean, 10 k permutations) (Additional file 1: Table S6, Additional file 2: Fig. S2B) compared to autosomes. Altogether, these characteristics may be explained by high recombination rates on the PARs [50], but currently, recombination maps are not available for any red algal species to confirm this hypothesis.
To summarize, our analyses showed that the SDRs of representative Gracilaria species are small and present typical structural characteristics of non-recombining regions, such as lower gene density and higher repeat content, but they exhibit little genetic degeneration. Contrary to the PARs of other UV systems [8], we found no evidence for distinctive features in the red algal PARs compared to autosomes [8].
Evolutionary history of the sex-determining regions in Gracilaria
To study the evolutionary history of the Gracilaria sex chromosomes, we first examined the genes present in the male and female SDRs. The male haplotype harbors between 13 (G. caudata) and 38 (G. gracilis) protein-coding genes, whereas the female haplotype contains between 12 (G. vermiculophylla) and 33 (G. caudata) genes (Additional file 1: Table S5). The sex-linked genes could be further described as “sex-limited” when present only in one haplotype of the sex locus, or “gametologs” when they have homologs in the SDRs of the opposite sex. The presence of the gametologs is consistent with the two SDR haplotypes having evolved from a common ancestral autosomal region. The proportion of gametologs among the total sex-linked genes was relatively small and varied across species, ranging from 6/46 genes in G. caudata, 8/52 genes in G. gracilis, 8/35 genes in G. chilensis, and 22/40 genes in G. vermiculophylla (Fig. 3A). In contrast, sex-limited genes may represent ancestral gametolog pairs, where one of the homologs was lost by the haplotypic counterpart. Alternatively, they may have been acquired by the U-SDR or V-SDR specifically, sometime after the sex-determining regions stopped recombining. To gain further insight into which scenario is more likely, we analyzed the conservation of sex-linked gene content across the four Gracilaria species, as well as the overall homology of these genes across the species.
Conservation of genes in the sex-determining regions (SDRs) in the four Gracilaria species. A Schematic representation of the gene content and syntenic relationships among female-linked genes (pink) and male-linked genes (blue) in the four studied species. B Venn diagram illustrating the orthogroups conservation across the four species (female conserved orthogroups in pink and male conserved orthogroups in blue). See also Additional file 1: Table S8. C Gene trees of conserved homologous gametolog pairs. Bootstrap values (from 1000 resamplings) are shown next to the branches. D Synonymous substitutions (KS) between gametolog pairs per species. Orthogroup names are indicated for homologous gametologs conserved in at least two species. E Microsynteny plot of the V-sex determining regions (V-SDRs) in the three species with continuous assemblies. Conserved SDR genes shared by all three species are marked in yellow, and genes shared by two species are marked in green. Blue box marks the newly acquired SDR genes in G. vermiculophylla
We performed an orthology search for protein-coding genes and inferred a total of 9945 orthogroups (OGs) with 101 of these orthogroups containing at least one sex-linked gene in one of the four species (Additional file 1: Table S8). Interestingly, most of the latter orthogroups (80) contained sex-linked genes in only one of the species (orthologs in other species being autosomal), and regardless of the evolutionary distance, comparisons resulted in very few shared sex-linked OGs (Fig. 3B). This result suggests that the gene content within the Gracilaria sex locus is highly dynamic (Additional file 2: Fig. S3A, Additional file 1: Table S9). The similarity between sex-limited SDR genes and their closest autosomal paralog would be consistent with gene duplication events (i.e., after the divergence of the U and V) that created either the SDR or the autosomal copy. Given that the expression levels of sex-limited genes (log2(TPM + 1)) were significantly lower than those of their autosomal paralogs (paired t-test, t = − 2.3133, df = 66, p-value = 0.02383) (Additional file 1: Table S9, Additional file 2: Fig. S3A), this suggests that the functions of these sex-limited genes may be partially compensated by their autosomal counterparts. This could indicate a potential gradual loss of SDR copies, with the autosomal paralogs arising as functional gene rescues. Alternatively, the SDR copies could be expressed at distinct developmental stages and have specialized roles not captured in our data, as our analysis focuses only on expression in sexual gametophytes. In cases where no autosomal homologs were found, it may point to a species-specific relocation of genes to the SDR or that these could be ancestral SDR genes lost entirely in all the other species (Additional file 1: Table S8).
Despite low conservation of SDR gene content, we found two orthogroups with gametolog pairs that were shared across all four species (OG0000231, OG0000351) and two with conserved, homologous male gametologs in all species (OG0000537, OG0000538), although the female gametologs have been lost in some of the Gracilaria species. Phylogenetic relationships among these conserved homologous gametologs indicate that the SDRs in Gracilaria are ancestral and stopped recombining before the speciation event, approximately 162–126 million years ago [35] (Fig. 3C). Remarkably, genes from OG0000231 shared orthology with a male sex-linked gene (Bsm1) identified in another red alga, Bostrychia moritziana, from the order Ceramiales [32, 51], which diverged from Gracilariales around 390 million years ago (MYA) [36]. This finding indicates that the sex determining regions in red algae may be even more ancient. In contrast, no conserved orthogroups were shared with SDR genes identified in Pyropia haitanensis (order Bangiales) [33], a lineage that last shared a common ancestor with Gracilaria more than 817 million years ago [36]. Additionally, we found six male-limited genes that were consistently sex-linked in all Gracilaria species (Fig. 3B, Additional file 1: Table S8). These genes may represent the ancestral V-sex chromosome. Notably, no conserved female-limited genes were found on the U chromosome.
Finally, we examined the rates of divergence between gametologs to uncover the volutionary strata’ of sex chromosomes. By analyzing the synonymous substitution rates (KS) between gametologs within species, we can estimate the relative timing of successive recombination suppression between U and V sex chromosomes, providing insight into the history and evolutionary dynamics of the SDR regions. Gene-by-gene analysis showed, as expected, that the KS is highest and saturated for the gametologs associated with the ancestral SDR in all four species (OG000231, OG000351, OG0000538, OG0000537) (Fig. 3D, Additional file 1: Table S10). Based on gene trees and KS values for two additional gametolog pairs in G. caudata (OG000670) and G. gracilis (OG0001184), it appears that these genes became sex-linked at later stages, after which female gametologs were lost, and some male genes relocated outside the V-SDR (Fig. 3D). Lastly, the SDR region in G. vermiculophylla included a cluster of nine consecutive gametologs (based on the physical position in the male sex chromosome assembly) with lower KS values (< 1), located in an inverted region in relation to the G. chilensis SDR (Fig. 3D, E). Three of these gametologs share homology with male sex-linked genes in G. gracilis, but the phylogenetic analysis suggests independent acquisition, as the genes follow the species tree, with male and female copies grouping together in G. vermiculophylla (Additional file 2: Fig. S4, Additional file 1: Table S8). The remaining gametologs are either not sex-linked in any other Gracilaria species or lack orthologs altogether, suggesting a recent U- and V-SDR expansion in G. vermiculophylla which is further supported by synteny analysis (Fig. 3E).
Taken together, our analyses suggest that the sex chromosome system is at least 100 million years old and ancestral to the Gracilaria species studied here. The fact that one ancestral gametolog, shared by all four species, is also sex-linked in the distantly related red alga Bostrychia moritziana, suggests that this sex chromosome system may have already originated over 390 MYA [32, 36, 51]. Although we did not identify clear evolutionary strata on the sex chromosomes, we observed a recent SDR expansion in G. vermiculophylla involving a series of inversions and/or gene shuffling events.
Sex differences in gene expression and function
We focused on the functions and expression patterns of sex-linked genes, particularly the conserved ones, assuming that these genes play a pivotal role in the cascade of differentiation and development in males and females. Gene ontology analysis revealed that genes in sex-linked scaffolds across all four species were predominantly associated with DNA binding, phosphorylation, and energy metabolism Additional file 1: Table S5, Additional file 2: Fig. S5). Particularly, conserved male-specific genes included putative key regulators, such as a homeobox transcription factor (OG0003259), a glycosyltransferase (OG0000538), and a serine/threonine-protein kinase (OG0005139). The conserved gametolog pairs were found to be involved in signaling pathways and transcription regulation (e.g., OG0000231: Zinc finger and ankyrin-repeat protein, OG0000351: ATP-NAD kinase, PpnK-type) (Additional file 1: Table S5). The significant correlation in transcript abundances between gametolog pairs (log2(TPM + 1)) in males and females (Pearson’s r = 0.541, p-value = 0.0248, Additional file 2: Fig. S3B) suggests that these genes have been retained on the U- and V-SDRs due to their crucial roles during the haploid phase of the life cycle. These roles may be related to sex determination or reproduction but also involve broader gametophytic developmental processes. Therefore, we examined differential expression levels of autosomal genes between male and female gametophytes to identify sex-biased genes that could be regulated by the SDRs.
To explore the extent of sex-biased gene expression, we analyzed transcriptomic data from male and female gametophytes. Gracilaria male and female gametophytes are morphologically difficult to distinguish prior to the development of reproductive structures, where sex becomes easily identifiable after the development of spermatangial sori (male reproductive structures) or the carposporophyte (i.e., cystocarp) (Fig. 4A). However, in some cases, male gametophytes have been reported to be smaller, as seen in G. gracilis and G. chilensis, and may exhibit physiological differences, such as variations in growth speed [52,53,54]. Our analysis included gametophytes representing two developmental stages: younger gametophytes that had just started forming reproductive structures (G. chilensis and G. caudata) and fully grown mature gametophytes (G. vermiculophylla and G. gracilis). For G. chilensis and G. caudata, gametophytes were cultured under controlled laboratory conditions and sex determination was conducted through microscopic observation of male sori or female trichogynes (receptive structures associated with the female gametangia). These gametophytes were young, and because they were grown in isolation, females were fertile but unfertilized and did not develop cystocarps. In contrast, gametophytes of G. vermiculophylla and G. gracilis were collected directly from the field. These thalli were fully mature, and females were fertilized, bearing cystocarps (which were removed from the tissue prior to RNA extraction).
Sex-biased expression in Gracilaria. A Gametophytes of Gracilaria showing male gametophytes with the presence of spermatangial sori and female gametophytes with carposporophytes (cystocarps). G. caudata and G. gracilis gametophytes are shown in their natural habitat. Male and female thalli of G. chilensis and G. vermiculophylla are shown side by side, highlighting how morphologically similar the two sexes appear to the naked eye. B Proportion of sex-biased genes in Gracilaria species: green, male-biased genes; purple, female-biased genes; grey, unbiased genes. C Venn diagrams representing the number of shared orthogroups (OGs) containing male-biased genes (top) and female-biased genes (bottom) in all four studied species. D Expression levels (log2(TPM)) of sex-biased and unbiased genes in male and female gametophytes of Gracilaria: green, male-biased genes; purple, female-biased genes; grey, unbiased genes. Significant differences are indicated with asterisk, pairwise Wilcoxon test with Holm correction, ***p-value < 0.0001, **p-value < 0.001, *p-value < 0.05). Photo credit: E.M. Plastino/F. Marchi (G. caudata), C. Destombe (G. gracilis), M.L. Guillemin (G. chilensis), S.A. Krueger-Hadfield (G. vermiculophylla)
We observed significant variability in the proportion of the transcriptome that experiences sex-biased expression patterns depending on the developmental stage of the gametophytes (Fig. 4B, Additional file 1: Table S11). In G. caudata and G. chilensis, only a small number of genes were sex-biased—336 (3.9% of total predicted genes) and 611 (7.8% of total predicted genes), respectively—whereas G. gracilis and G. vermiculophylla showed moderate levels of biased expression, with 2213 (23.5% of total predicted genes) and 1475 (21.8% of total predicted genes), respectively. Most sex-biased genes were found in orthogroups unique to each species, indicating that sex-biased expression is highly species-specific. Only four orthogroups exhibited conserved male bias, while two orthogroups showed conserved female bias (Fig. 4C). However, we observed a greater number of conserved orthogroups between species pairs representing similar developmental stages: young gametophytes in G. chilensis and G. caudata versus fully mature and fertilized gametophytes in G. vermiculophylla and G. gracilis. Functional analysis revealed that male-biased genes of both G. chilensis and G. caudata were significantly enriched in Gene Ontology (GO) terms related to cell division, DNA replication, and cytoskeletal organization (GO:0044774, GO:0000727, GO:0007020, GO:1,905,775, GO:0006267, GO:1,902,975, GO:0006268) (Additional file 1: Table S12). These GO terms are characteristic of a highly proliferative tissue, suggesting they could be related to spermatial production in the young male gametophytes. In contrast, female-biased genes of G. vermiculophylla and G. gracilis female gametophytes were enriched in GO terms involved in response to external stress and defense-related processes (GO:0098869, GO:0045454, GO:0042744) (Additional file 1: Table S12). This implies a protective role against environmental challenges, such as intense light exposure, pathogen attacks, or other stressors that trigger reactive oxygen species (ROS) production.
Although the specific identities of sex-biased genes varied between species, a pattern known as turnover of sex-biased gene expression—a common pattern emerged across all Gracilaria species when we examined the expression levels (log2(TPM + 1)) of sex-biased and unbiased genes in males and females. In female gametophytes, female-biased genes (FBGs) were expressed at significantly higher levels, while male-biased genes (MBGs) were expressed at significantly lower levels compared to unbiased genes. In contrast, male gametophytes exhibited more consistent expression levels across male-biased, female-biased, and unbiased genes, suggesting that sex-biased expression primarily results from the up-regulation or down-regulation of specific genes in females (Fig. 4D).
Finally, we examined the evolutionary forces shaping sex-biased genes, focusing on their chromosomal distribution and patterns of coding sequence evolution, measured as non-synonymous to synonymous substitution rates (KN/KS). Sex-biased genes are often expected to evolve more rapidly than unbiased genes due to positive selection or relaxed evolutionary constraints [55,56,57,58]. Additionally, sex-biased genes might be enriched in PARs if their partial linkage to the SDR provides reproductive or survival advantages, facilitating the spread of advantageous alleles [59]. However, our analysis revealed no significant enrichment of sex-biased genes in the PARs of the G. vermiculophylla sex chromosome (Additional file 1: Table S13, chi-square test, p-value > 0.05). Furthermore, the evolutionary rates of male- and female-biased genes (measured as KN/KS) did not differ from those of unbiased genes, indicating that these genes are likely evolving under purifying selection (Additional file 1: Table S14, Additional file 2: Fig. S6).
In sum, although Gracilaria SDRs contain several conserved genes enriched in functions related to signaling and transcription regulation, the downstream networks of sex-biased gene expression are highly variable in this group. Sex-biased genes experience high turnover rates, with little conservation of sex-biased expression among orthologous genes, indicating that such expression is highly species-specific. However, sex bias consistently appears to stem from gene expression regulation in females, suggesting that regulatory elements may play a crucial role in sexual differentiation.