Selection drives convergent gene expression changes
during transitions to co-sexuality in haploid sexual
systems

Guillaume Cossard, Olivier Godfroy, Zofia Nehr, Corinne Cruaud, J. Mark
Cock, Agnieszka Lipinska, Susana Coelho

► To cite this version:

Guillaume Cossard, Olivier Godfroy, Zofia Nehr, Corinne Cruaud, J. Mark Cock, et al. Selection drives convergent gene expression changes during transitions to co-sexuality in haploid sexual systems. Nature Ecology & Evolution, 2022, 6, pp.579-589. 10.1038/s41559-022-01692-4 . hal-03619076

HAL Id: hal-03619076

https://cnrs.hal.science/hal-03619076

Submitted on 24 Mar 2022

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
1 SELECTION DRIVES CONVERGENT GENE EXPRESSION CHANGES

2 DURING TRANSITIONS TO CO-SEXUALITY IN HAPLOID SEXUAL

3 SYSTEMS

4 Guillaume G. Cossard<sup>1,3</sup>, Olivier Godfroy<sup>1</sup>, Zofia Nehr<sup>1</sup>, Corinne Cruaud<sup>2</sup> J. Mark Cock<sup>1</sup>, Agnieszka

5 Lipinska<sup>1,3</sup>, Susana M. Coelho<sup>1,3*</sup>

6 <sup>1</sup>Sorbonne Université, UPMC Univ Paris 06, CNRS, Integrative Biology of Marine Models, Station Biologique de

7 Roscoff, CS 90074, F-29688, Roscoff, France <sup>2</sup>Genoscope, Institut de Biologie François Jacob, CEA, Université Paris-

8 Saclay, Evry, France <sup>3</sup>Max Plank Institute for Biology Tübingen, Max-Planck-Ring 5, 72076, Tübingen, Germany

9 *susana.coelho@tuebingen.mpg.de

10 ABSTRACT

11 Co-sexuality has evolved repeatedly from unisexual (dioicous) ancestors across a wide range of taxa.

12 However, the molecular changes underpinning this important transition remain unknown, particularly in

13 organisms with haploid sexual systems such as bryophytes, red algae and brown algae. Here we explore

14 four independent events of emergence of co-sexuality from unisexua lancestors in brown algal clades to

15 examine the nature, evolution and degree of convergence of gene expression changes that accompany

16 the breakdown of dioicy. The amounts of male versus female phenotypic differences in dioicous species

17 were not correlated with the extent of sex-biased gene expression, in stark contrast to what is observed

18 in animals. Although sex-biased genes exhibited a high turnover rate during brown alga diversification,

19 some of their predicted functions were conserved across species. Transitions to co-sexuality consistently

20 involved adaptive gene expression shifts and rapid sequence evolution, particularly for male-biased

21 genes. Gene expression in co-sexual species was more similar to that in females rather than males of

22 related dioicous species, suggesting that co-sexuality may have arisen from ancestral females. Finally,

23 extensive convergent gene expression changes, driven by selection, were associated with the transition to

24 co-sexuality. Together, our observations provide insights on how co-sexual systems arise from ancestral,

25 haploid UV sexual systems.

26

27

28
29 INTRODUCTION

30 Eukaryotic organisms exhibit a wide diversity of sexual systems, ranging from separate sexes (referred to
31 as gonochorism in animals and dioecy in plants) to co-sexuality (combined sexes) and several theories
32 have been developed to explain what conditions favour which strategy$^{1-7}$. The evolution of this diversity
33 often involved transitions between sexual systems. For example, separate sexes have evolved from co-
34 sexual ancestors independently many times in several eukaryotic lineages, and the fundamental
35 mechanisms and evolutionary drivers of this important transition have been intensively studied in many
36 organisms (reviewed in $^{2,8}$). Frequently, organisms with separate sexes display marked sexual dimorphism
37 in a range of morphological, behavioral and physiological traits. Females and males are nevertheless
38 genetically similar with the exception of the sex-specific regions of their sex chromosomes. While sex-
39 chromosomes necessarily play a role in the expression differences between sexes, most of sex-biased
40 gene expression involves autosomal genes $^{9-11}$. Differences in autosomal gene expression patterns
41 between sexes may be associated with different physiological processes directly linked to the production
42 of male or female gametes (primary sexual dimorphism) or to the consequences of sexual selection
43 and/or sexual specialization (secondary sexual dimorphism) that may occur once separate sexes have
44 evolved $^{12}$.

45 While the emergence of separate sexes from co-sexual ancestors and the evolution of sexual dimorphism
46 have been thoroughly investigated $^{11,13-15}$, less attention has been devoted to the opposite transition, i.e.
47 from separate sexes to co-sexuality. Transitions to co-sexuality have occurred frequently during
48 eukaryotic evolution and are relatively common in animals (e.g.$^{13,16-20}$). In flowering plants this transition
49 was believed to be rare but recent studies are increasingly providing evidence that dioecy-to-monoecy
50 transitions may have occurred frequently $^{21,22}$. Evolutionary models intending to decipher the causes of
51 such transitions invoke the sex-allocation theory 5 and the deterministic fate of genetic modifiers causing
52 the acquisition of an opposite-sex function $^{23,24}$. However, empirical knowledge on the proximate
53 mechanisms and forces driving the shift from separate sexes to co-sexuality remains largely elusive.

54 Transitions from separate sexes to co-sexuality are also prevalent in eukaryotic lineages other than
55 animals and flowering plants, and in particular those that express sex during the haploid stage of their life
56 cycles. In organisms such as bryophytes, liverworts, green, red and brown algae, male and female sexes
57 are expressed during the haploid (gametophyte) stage 25. The terms 'dioicy' (i.e., separate sexes during
58 the haploid phase of the life cycle, as opposed to 'dioecy' where separate sexes occur in the diploid
59 phase) and monoicy (i.e., co-sexuality during the haploid phase of the life cycle, as opposed to 'monoecy'
60 where co-sexuality occurs in the diploid phase) are used to describe the sexual systems of these
61 organisms 26. Genetic sex determination in dioicous organisms occurs during meiosis (and not at
62 fertilisation as in XY and ZW systems) 27, depending on whether spores inherit a U or V sex-chromosome
63 26,28. Spores receiving a V chromosome will develop into a male multicellular individual (male
64 gametophyte) and the spores inheriting a U chromosome will grow into females (female gametophytes).
65 Organisms with haploid sex determination may also produce male and female sexual structures in the
66 same (co-sexual) individual (monoicy) 29,30. Despite the prevalence of haploid sexual systems among
67 eukaryotes, the gene expression changes and evolutionary forces underlying transitions from dioicy to
68 monoicy have remained largely unknown.
In this context, the brown algae represent a particularly attractive group for studies of the evolution of sexual systems and breakdown of dioicy. The brown algae are a complex multicellular lineage that is part of the stramenopile (or heterokont) supergroup, which also includes diatoms and oomycetes, and they have diverged from the Archaeplastida lineage at the time of the eukaryotic crown radiation ³¹. Most brown algae have a haplo-diplontic life cycle, with a haploid gametophyte generation alternating with a diploid sporophyte generation. In these brown algae, sexuality is expressed in the haploid generation, with male and female gametes either produced by the same haploid individual (monoicy) or on two separate haploid individuals (dioicy). Dioicy is the prevalent reproductive system ²⁹,³². This situation contrasts markedly with that described for flowering plants, where only about 6% of extant species have separate sexes, and is more similar to that of bryophytes and liverworts ³⁰. Dioicous brown algae may exhibit a broad range of levels of sexual dimorphism, both at the level of the gametophytes but also with respect to the difference between male and female gametes size ²⁹,³². While the predicted ancestral state in the brown algae is dioicy, transitions to monoicy have occurred frequently and independently in the different clades ³²,³³. The independent emergence of monoicous lineages from dioicous ancestors makes this group particularly interesting to examine the genomic consequences and mechanisms underlying the breakdown of dioicy.

Here, we explore multiple, repeated events of loss of dioicy (Figure 1) to investigate the molecular basis and level of convergence of the shifts to co-sexuality. We test the hypothesis that sexually dimorphic algae might be expected to have more sex-biased genes and, because dioicy is ancestral, we predicted that similar gene sets would be sex-biased across all the dioicous species. Contrary to our prediction, we demonstrate a lack of correlation between phenotypic sexual dimorphism and gene expression levels among dioicous brown algae. Ancestral state reconstruction indicated high turnover rates of sex-biased genes, yet independently recruited sex-biased genes shared similar functions across the species. To characterise the molecular changes associated with the evolution of monoicy, we then focused on modifications in gene expression patterns of orthologous genes that are specifically or preferentially expressed in haploid males and females of a dioicous species, when they function in a monoicous context. Male-biased genes were particularly concerned by both adaptive expression shifts and faster evolutionary rates associated with the transition to monoicy. Monoicous species displayed expression profiles that were more similar to those of the female of the closely related dioicous species than to the male. Finally, we identified a pronounced level of convergent gene expression changes associated with the emergence of co-sexuality, which were likely driven by selection.

## RESULTS

The present study examines sex-biased gene expression in dioicous brown algae and the gene expression changes associated with the transition from dioicy to monoicy. We based our analysis on transcriptomes sequenced from pairs of dioicous-monoicous species in four major clades of brown algae spanning app. 200 million years of evolution ³⁴. The transitions are predicted to have occurred at different times in the past (between 20 and 88 MY; Figure 1). Each pair represents an independent transition from dioicy to monoicy. We chose dioicous species with different levels of gamete dimorphism, reflecting the diverse levels of gamete dimorphism occurring across brown algae.
108 SEX-BIASED GENE EXPRESSION IN DIOICOUS BROWN ALGAE

109 Gene expression patterns in gametophytes of the eight brown algal species were measured by deep
110 sequencing (RNA-seq) of cDNA from male, female and co-sexual gametophytes. Transcript abundance
111 (measured as transcripts per million, TPM) was strongly correlated between biological replicates with r²
112 ranging from 0.89 to 0.99 (Supplementary Table 1). Counts of expressed genes (TPM>5th percentile
113 counts across all genes in at least one sample) identified a number of expressed genes that ranged from
114 13,180 to 27,391 (Supplementary Table 1).

115 Deseq2 was used to identify genes that were differentially expressed in each of the sexes of the dioicous
116 species ³⁵. The analysis retained only genes that displayed at least a 2-fold change expression level
117 between sexes (FC>2, p<0.05). Note that sex-linked genes (genes located in the sex-specific regions on
118 the V (male) and U (female) sex chromosomes; see methods), were removed from the set of sex-biased
119 genes and thus excluded from further analysis.

120 All four dioicous brown algae displayed substantial sex-biased gene expression, at least compared with
121 plants and other brown alga ¹³,¹⁵,³⁶ ranging from 12.71 % of the expressed genes in S. rigidula to 33.17% in
122 S. firma (Figure 2A-2B, Supplementary Tables 2-3). We found similar proportions of male-biased
123 compared with female-biased genes for the majority of the studied species (Figure 2A-2B) with the
124 exception of S. polyschides, where male-biased genes were more abundant than female-biased genes
125 (16.51% male-biased genes versus 9.39% female-biased genes; Chi²-test P < 2.2 x10⁻¹⁶).

126 SEX-BIASED GENE EXPRESSION AND PHENOTYPIC SEXUAL DIMORPHISM

127 To investigate the link between sex-biased gene expression and the level of sexual dimorphism, we
128 carried out morphometric measurements of male and female gametophytes complemented with
129 literature searches. These measurements allowed us to quantify the amount of phenotypic dimorphism
130 present in each of the four dioicous species (Supplementary Table 3, Figure 2C). In all dioicous species,
131 gamete size dimorphism was coherent with sexual differences in terms of gametophyte cell size
132 (Supplementary Table 3). For example, *D. herbacea* gametophytes presented marked sexual dimorphism
133 both at level of gamete size and gametophyte cell length, whereas *S. firma* was the species with least
134 sexual difference both in terms of gametophyte morphology and gamete size (Supplementary Table 4,
135 Figure 2C-D).

136 In animals, sexual differences at the phenotypic level are correlated with levels of sex-biased gene
137 expression ¹⁴,³⁷, but this correlation has not been found in plants ³⁶. We compared the differences in
138 gametophyte cell size between males and females with the proportion of sex-biased genes in each of the
139 four dioicous brown algal species. We detected no correlation between phenotypic sexual dimorphism
140 (gametophyte cell size) and the number of sex-biased genes (Figure 2E). For instance, *S. firma* was the
141 species that exhibited the highest level of sex-biased gene expression and nonetheless presented the
142 lowest level of phenotypic sexual dimorphism. Taken together, our observations indicate a considerable
143 level of sex-biased gene expression in the four dioicous species studied here, but the level of sex-biased
144 gene expression did not reflect the level of morphological dimorphism between males and females.
145 EVOLUTION OF SEX-BIASED GENE EXPRESSION IN DIOICOUS SPECIES

146 We next investigated how sex-biased gene expression has evolved by comparing the four dioicous brown
147 algal species. Orthofinder identified a total of 14,017 orthogroups (OGs) across the dioicous species, of
148 which 2,098 contained only one gene per species and therefore represented the set of 1:1:1:1 OGs. An
149 additional 2,778 OGs had a single member in each of three of the studied species (i.e., the gene was
150 missing in the fourth species). We considered that these 1:1:1:0 OGs, which likely represent single copy
151 ancestral genes that were lost in one of the species, also provide useful information about conservation
152 of sex-biased gene expression. Note that the 1:1:1:0 OGs could also represent OGs where one of the
153 genes is missing from one of the genome assemblies, particularly the draft genome assembly for S.
154 *rigidula*. Furthermore, we also included 1,085 orthogroups with a duplicated gene in a single species
155 (1:1:1:2 OGs) that aligned along more than 60% of their length, resulting in 5,961 'dioicous single-copy
156 orthologs' (DSOs; Supplementary Table 5, Extended Data Figure 1).

157 We then used maximum likelihood approaches to infer the ancestral states of sex-biased gene expression
158 across these dioicous species (Figure 3B). Our analysis identified very few genes that were predicted to be
159 ancestrally sex-biased, with the vast majority having evolved sex-bias at some point along the branches.
160 Among the 2,116 sex-biased DSOs in at least one species, only 43 (2.03%) were inferred to be sex-biased
161 in the last common ancestor of the four brown algal species (Figure 3). Accordingly, no DSOs were
162 consistently sex-biased across the four species (not different from what is expected by chance, exact test
163 multi-set intersection P = 0.506). A total of 139 OGs exhibited a bias in one species that was inconsistent
164 with the direction of bias observed in at least one other species (Supplementary Table 5).

165 Although the above analysis showed that sex-biased genes were not conserved among the four species, we
166 examined if sex-biased genes in different species were involved in similar functions, by comparing gene
167 ontology (GO) terms of sex-biased genes across species using Blast2Go³⁸. We detected significant
168 enrichment of GO terms for biological processes related to 'ion- and transmembrane transport' and
169 'cilium' often associated to male-biased genes across dioicous species. Conversely, the sets of female-
170 biased genes of all the studied species were enriched for GO terms related to oxidation/reduction
171 (Extended Data Figure 2, Supplementary Table 6). Taken together, our results indicate that whilst sex-
172 biased genes exhibited a high turnover rate during brown alga diversification, some of their predicted
173 functions were conserved across dioicous species.

174 We also asked whether sex-biased gene expression emerged in dioicous species as a result of random
175 expression evolution under low selective pressure for non-pleiotropic genes³⁶ or rather as a consequence
176 of sexual selection. To distinguish between these two possibilities, we computed phylogenetic
177 independent contrasts (PICs) of sex-biased genes across species in which those genes are not sex-biased
178 versus unbiased genes (Extended Data Figure 3). We found that PICs differed slightly between unbiased
179 genes and orthologs of sex-biased genes in species in which those genes are not sex-biased (Mann-
180 Whitney Rank test *P* = 0.0495). This result indicates that genes that evolved sex-bias may have done so
181 because they already experienced low constrains on their expression levels, possibly due to lower
182 pleiotropic expression patterns³⁶,³⁹, although we cannot exclude that sexual selection may be also
183 involved in the emergence of sex-biased gene expression in brown algae.
185 To study changes in sex-biased gene expression that accompany the transition from dioicy to monoicy,
186 we first identified single-copy orthologous genes for each of the four dioicous-monoicous sister species
187 pairs (pairwise single-copy OGs; PSOs, Figure 4A). We were able to infer between 6,109 and 11,953 PSOs
188 for each of the four pairs of species (Figure 4A; Supplementary Tables 7-11). PSOs were classified as being
189 sex-biased or unbiased by comparing male and female expression in each dioicous species (FDR<0.05,
190 FC>2). We then examined the patterns of expression of male-, female- and unbiased PSOs in dioicous
191 males and females and in the corresponding monoicous species.

192 In three out of the four species pairs, the levels of expression of sex-biased genes in the monoicous
193 species were similar to the values measured for orthologs in females of the corresponding dioicous
194 species (Figure 4A, Extended Data Figure 4). In these three pairs, male-biased genes were downregulated
195 in the monoicous species compared with males, and they displayed similar expression levels to male-
196 biased genes in females of the dioicous species, suggesting that de-masculinisation of gene expression of
197 the monoicous species counterpart occurred frequently. Female-biased genes were expressed at similar
198 mean levels in S. firma females compared with the corresponding monoicous species C. linearis. In the S.
199 polyschides-S. dermatodea pair of species, female-biased genes had a similar pattern in males and
200 monoicous. Both female- and male-biased genes in D. herbacea showed significantly different mean
201 expression levels compared with D. dudresnyi. In the S. rigidula/H. paniculata species pair, no significant
202 difference was detected between the expression of sex- biased and unbiased genes between the two
203 species. Note, however, that results for S. rigidula/H. paniculata were more difficult to interpret, as the
204 low number of sex-biased among the PSOs precluded robust statistical analysis.

205 We next investigated the gene expression profiles of monoicous species in order to test whether their
206 transcriptional patterns resemble those of their male or female dioicous counterparts. We computed the
207 Pearson product-moment coefficient of regressions of gene expression profiles (in log2(TPM+1)) of males
208 or females compared with that of the monoicous species within each species pair. We compared Pearson
209 correlation coefficients for both sex-biased genes and unbiased genes in both males or females,
210 considering sex-biased genes in males and females as independent groups. We also compared the
211 correlations of expression profiles with the orthologs of sex-biased and unbiased genes in the monoicous
212 species, separately for males and females. These groups of sex-biased versus unbiased genes being
213 expressed within the same individuals, we considered them as dependant groups in the *cocor* package⁴⁰.
214 Altogether, these analyses indicated that, with the exception of the *S. rigidula-H. paniculata* species pair,
215 the gene expression profiles of the monoicous species were significantly more similar to those of females
216 than they were to male profiles (Figure 4B; Extended Data Figure 5). Moreover, the close association
217 between female and monoicous expression profiles was observed for both sex-biased and unbiased
218 genes specifically in *Sacchoriza* and *Desmarestia* species pairs (Figure 4B, black asterisks at the top;
219 Extended Data Figure 6).

220 Interestingly, with the exception of the Ectocarpales species pair (*S. firma* - *C. linearis*), sex-biased gene
221 expression profiles diverged significantly less from the monoicous species than did that of the unbiased
222 genes (Figure 4B, S5). Overall, the expression profile similarity observed between females and monoicous
individuals were mainly driven by expression patterns of male-biased genes, with the exception of the *Desmarestia* species pair (Figure 4A, S5). We also noted that the highest similarity indexes for within species pairs were found for the species with the lowest level of sex-biased gene expression (*S. rigidula*), and the lowest similarity was observed for *S. firma*, the species with the highest level of sex-biased gene expression (Figure 4B). The relatively high proportion of sex-specific genes present in *S. firma* (Figure 2A) is unlikely to be the cause of the observed low similarity, because only 29 (0.49%) of the sex-specific genes are among the PSOs used for the Pearson similarity analysis.

Taken together, the above observations suggest that gene expression profiles of monoicous species tend to be more closely related to the females of the related dioicous species, and this similarity appears to be driven by sex-biased genes, in particular male-biased genes. The tendency to reproduce the female transcriptome in the monoicous species was repeatable in independent transitions to co-sexuality.

### IS SELECTION INVOLVED IN EXPRESSION CHANGES DURING TRANSITION TO MONOICY?

To examine whether changes in gene expression during transition to co-sexuality were the result of selective or neutral forces, we computed the degree of directional selection using $\Delta x$. This parameter evaluates the divergence in expression level in relation to the variation in expression level seen across replicates $^{11,36,41}$. We computed $\Delta x$ of the PSO sets, separately for each pair of species and reported the proportions of orthologs with an absolute $\Delta x > 1$, i.e., orthologs whose expression shift is attributable to directional selection (Supplementary Table 11, Figure 4C). Depending on the species pair, between 10.8% and 40.1% of unbiased genes exhibited expression shifts attributable to selection ($|\Delta x|>1$) (Supplementary Table 12). We then asked whether male- and female-biased genes were preferentially concerned by adaptive expression shifts during transitions to monoicy compared with unbiased genes. Figure 4C illustrates the proportion of orthologs with sex-bias displaying $|\Delta x| > 1$ (in other words, under putative directional selection) and how sex-biased genes are more likely to display $|\Delta x| > 1$ in comparison to the unbiased orthologs. Fisher's exact tests (asterisks in Figure 4C) showed that for three out of four species pairs, male-biased genes were indeed more likely to display $|\Delta x|>1$ than unbiased genes (Figure 4C). This was also the case for female-biased genes in *S. polyschides-S. dermatodea* pair of species (Fisher exact tests, $P<2.2 \times 10^{-16}$ in both sexes) and *D. herbacea-D. dudresnayi* pair (Fisher exact tests, $P = 3.9 \times 10^{-3}$ and $P = 1.8 \times 10^{-5}$ in females and males, respectively). In *S. polyschides-S. dermatodea* and *S. rigidula-H. paniculata*, female-biased genes showed lower levels of adaptive evolution of expression compared with unbiased genes (Supplementary Table 11, Figure 4C). Taken together, our observations indicate that male-biased genes preferentially exhibit a shift in expression during the transition to monoicy that may be explained by directional selection.

We also assessed if evolution of gene expression during the transition to monoicy has been driven by DNA sequence evolution, by using measures of sequence divergence ($d_N/d_S$). We computed $d_N/d_S$ for male-biased, female-biased and unbiased genes for each of the dioicy-monoicy species pairs. For all four pairs, male-biased genes consistently exhibited higher evolutionary rates than female-biased and unbiased genes, although this difference was significant only for the *S. polyschides-S. dermatodea* species pair (Figure 4D, Supplementary Table 13). As this is the 'youngest' species pair (Figure 1), it appears that the level of sequence divergence during transition to monoicy is not associated with the age of transition.
262

263

264

Taken together, our observations indicate that shifts from dioicy to monoicy involved modifications to transcriptional patterns (expression divergence) mostly at male-biased genes that were likely driven by selection but also coding sequence evolution.

265

CONVERGENT GENE EXPRESSION CHANGES DURING TRANSITION TO MONOICY

266

267

268

269

In order to assess the extent to which gene expression changes occurring during the transition to monoicy were shared across the four species pairs, we focused on the single-copy orthologs across the eight species, herein termed ‘All Single-copy Orthologs’ (ASO). We found a total of 1,708 ASO (following the same approach as for DSO, see methods).

270

271

272

273

274

275

276

277

278

279

280

Among the 1,708 ASOs, 718 were sex-biased in at least one dioicous species (Supplementary Tables 14, 15). Sex-biased genes were not over-represented among ASOs (Fisher exact test, *P* = 0.097). Sixty one percent of the ASOs (1,043 out of 1,708) exhibited a conserved pattern of expression across all monoicous species compared to the dioicous species. This proportion was significantly different from what was expected by chance (permutation tests, *P* = 0.0255, 10,000 permutations) suggesting convergent gene expression changes during transition to monoicy across all studied pairs of species. Decomposition of variance components for the 1,708 ASOs detected a clear pattern of grouping of monoicous species, further illustrating the extensive convergence of gene expression during the transition from dioicy towards monoicy (Figure 4E). Functional analysis of genes that are convergently expressed during the transition to monoicy highlighted terms such as nucleic acid metabolic processes and transmembrane transport (Extended Data Figure 7).

281

282

283

284

285

About half (527) of the 1,043 genes that were consistently differentially expressed in monoicous versus dioicous species had a |Δₓ| >1, which is significantly more in proportion than among the rest of the ASO (290 genes with |Δₓ| >1 among 665 ASO, Fisher exact test *P* = 0.00543). This observation indicates that convergent gene expression changes may be associated with directional selection during the switch to monoicy.

286

287

288

289

290

291

292

We next tested whether sexual selection potentially occurring in males and females of dioicous species would be relaxed in monoicous individuals. This would translate by a reduction of purifying selection resulting in increased sequence divergence (increased *d*<sub>N</sub>*/d*<sub>S</sub>). Convergent genes (i.e., genes exhibiting a convergent pattern of gene expression in monoicous species) tended to exhibit faster divergence rates compared with non-convergent genes although the difference was not significant (permutation *t*-test, *P* = 0.0566; Extended Data Figure 8). Noteworthy, male-biased (but not female-biased) genes showed significantly higher *d*<sub>N</sub>*/*d*<sub>S</sub> than unbiased genes (Supplementary Table 16).

293

294

295

296

297

298

A likelihood ratio test of branch models (after Benjamin-Hochberg correction for multiple testing), identified 689 orthologs under positive selection on monoicous branches, 404 of which exhibited convergent gene expression changes. Orthologs under positive selection were over-represented among genes with convergent gene expression (Fisher exact test, *P* = 0.025). Taken together, these observations suggest that directional selection plays a role in driving changes in expression patterns during transitions to co-sexuality.
299

DISCUSSION

300

SEXUAL DIMORPHISM AND SEX-BIASED GENE EXPRESSION ARE UNCOUPLED

301

302

303

304

305

306

307

308

309

310

311

312

313

314

315

316

317

318

319

320

321

Morphological and physiological differences between males and females are ultimately due to divergences between sex chromosomes in species with genetic sex determination<sup>27</sup>, but the majority of morphological sexual dimorphism is thought to be associated with autosomal sex-biased gene expression<sup>9-11</sup>. Thus, it would be expected that species showing more prominent differences in morphology between male and female would also be characterised by high levels of sex-biased gene expression, as has been shown to be the case in birds<sup>37</sup>. Our study, in contrast, revealed no correlation between the level of sex-biased gene expression and the degree of phenotypic sexual dimorphism in the brown algae studied here. Therefore, the link between gene expression evolution and sexual selection is uncertain for these organisms, and sexual selection is likely not to be the main driver of the sex-biased gene expression evolution. This observation may reflect a lower degree of sexual selection in the brown algae compared with animals. Brown algae have relatively low levels of sexual dimorphism<sup>15,29</sup> and are broadcast spawners so the opportunities for mate choice and/or mating competition are mainly constrained to interactions involving male and female gametes<sup>42</sup>. Consistent with the idea that gamete sexual selection may occur, it has been shown recently that in the absence of males, female gametes of brown alga populations lose their sexual morphological characteristics, e.g. female gametes produce lower levels of pheromone and engage in parthenogenesis more rapidly<sup>43</sup>. Noteworthy, sex-biased genes found in male and female gametophytes of the model brown alga *Ectocarpus* show more rapid rates of divergence across species (measured as dN/dS) compared to unbiased genes, and their accelerated evolution has been, at least partly, attributed to positive selection<sup>15</sup>. These observations suggest that sexual selection plays a role in the evolution of sex-biased genes in brown algae, but may not be the only driver of sex-biased gene expression in this group of organisms.

322

SEX-BIASED GENES EXHIBIT FUNCTIONAL CONVERGENCE

323

324

325

326

327

Although dioicy is predicted to be the ancestral sexual system in brown algae<sup>32</sup> our results clearly indicate that sex-bias in the expression of individual genes is neither ancestral nor convergent. We found a very limited level of shared (ancestral) sex-biased gene expression across the studied brown algal species, and instead our data is consistent with lineage-specific recruitment of sex-biased genes. Our observations emphasize therefore a substantial turnover of sex-biased expression among brown algal genes.

328

329

330

331

Interestingly, our study suggests that sex-biased expression may have emerged on genes that were experiencing lower selective constraints on their expression level, possibly due to lower pleiotropy, in addition to the potential effect of sex-specific selection occurring after the evolution of separate sexes. A similar situation has been described recently in plants<sup>36</sup> and animals<sup>39</sup>.

332

333

334

335

336

Although the dioicous brown algal species studied here shared very few sex-biased genes, we found some level of convergence in terms of sex-biased gene function, at least for a subset of the sex-biased genes. These include biological functions that were previously found to be enriched in *Ectocarpus* gametophytes<sup>15,44</sup>, further underscoring the conservation of sex-biased gene function and supporting primary sexual dimorphic roles. These functions may be associated with sex-specific biological processes. For example,
enrichment in oxidation-reduction functions may relate to the more conspicuous growth of female gametophytes producing larger gametes which secrete a sperm-attracting pheromone, whereas cilia and ion transport functions are likely associated to the production of fast-swimming, bi-flagellated sperm by male gametophytes. Considering that brown algae share an ancestral sex chromosome, and that genes within the non-recombining sex determining region play a role in sex<sup>45</sup>, one possibility is that sexual characteristics in these UV systems mainly involve genes within the SDR together with a relatively limited number of autosomal genes involved in primary sexual dimorphisms. In other words, differences between sexes arise mainly from the different physiological processes directly linked to the production of male or female gametes rather than extensive sexual selection, sexual specialization and/or sexual antagonism (i. e., secondary sexual dimorphism)<sup>12</sup>.

347 FATE OF SEX-BIASED GENE EXPRESSION DURING TRANSITION TO MONOICY

348 Our sampling of species distributed across the brown algae phylogeny, associating pairs of related
349 dioicous and monoicous species, allowed us to trace the fate of sex-biased gene expression during
350 independent events of transition from dioicy to monoicy. With the exception of one species pair, sex-
351 biased genes exhibited adaptive expression shifts during the transition to monoicy. Male-biased genes,
352 specifically, were the main drivers of gene expression changes during the transition to monoicy, while
353 unbiased genes exhibited limited changes in pattern of expression with the switch in sexual system. In the
354 model brown alga *Ectocarpus*, RNA-seq analysis of multiple tissues and life cycle stages indicated that sex-
355 biased genes have restricted patterns of expression, which is a proxy for limited pleiotropy<sup>15</sup>. Pleiotropy
356 is known to restrict gene evolution, imposing stricter functional constraints on pleiotropic genes<sup>39,46</sup>. The
357 reduced pleiotropy of sex-biased compared with unbiased genes may increase their propension to
358 adaptively shift towards their evolved optimal expression profile during evolutionary transitions, in this
359 case the transition to monoicy<sup>10,39,47</sup>.

360 Sex-biased genes in dioicous brown algae such as *Ectocarpus* sp. typically display high evolutionary rates 361 compared to unbiased genes due either to directional selection or relaxed purifying selection <sup>15</sup>. With the 362 transition to monoicy, increased relaxation of sex-specific purifying selection acting on sex-biased genes 363 may be expected, leading to increased rates of sequence evolution. Accordingly, male-biased genes for all 364 species pairs presented faster evolutionary rates (although not significant for all species) during the 365 switch to monoicy, compared with female-biased or unbiased genes. This observation points to a shared 366 process of sexually antagonistic selection within dioicous species, especially in males, that allowed for 367 faster evolutionary rates of male-biased genes when relaxed during the transition from dioicy to monoicy.

368 *CONVERGENT CHANGES DURING BREAKDOWN OF DIOCY*

369 Convergent evolution, where a similar trait evolves in different lineages, provides an opportunity to study 370 the repeatability of evolution. In the brown algae, co-sexuality has repeatedly emerged from uni-sexual 371 ancestors <sup>32</sup>. We found that more than half (61%) of the orthologs across the four pairs of species 372 displayed similar expression shifts concomitant with the transition to monoicy, indicating that common, 373 independently acquired mechanisms are associated with co-sexuality. Remarkably, a substantial number 374 of these convergent genes (38.7%) were under positive selection, underlying the idea that convergent 375 changes associated with the shift of sexual system may be driven by comparable evolutionary pressures
376

across these distant species. Monoicous gametophytes were more closely related to females of the
corresponding dioicous species counterpart, suggesting, as proposed in volvocine algae ⁴⁸, that monoicy
may have arisen from ancestral females.

377

378

379

In our study, the expression profiles of gametophytes of all four monoicous species resembled those of the female gametophytes of their dioicous counterpart. Moreover, sex-biased genes tended to maintain the level of expression they had in dioicous species, suggesting that they retained their ancestral function in the context of derived monoicy. When their expression shifted, sex-biased genes, and especially male-biased genes showed signs of selection acting on their expression level to a greater extent than it acted on unbiased genes. Together, our results demonstrate that common mechanisms underlie the transition to monoicy across distant brown algal lineages and suggest that independent events of loss of dioicy may have involved acquisition of genes related to male development by a female gametophyte. The work presented here establishes therefore a framework for understanding at the genomic level how co-sexual systems arise from ancestral haploid UV sexual systems in the brown algae.

380

381

382

383

384

385

386

387

388

389

# MATERIALS AND METHODS

390

## SAMPLE PREPARATION, RNA EXTRACTIONS AND SEQUENCING

391

The algal strains used and sequencing statistics and BioProject accession number are listed in Supplementary Table 1. Gametophytes of all eight species were cultured at 13 °C in autoclaved natural sea water (NSW) supplemented with half-strength Provasoli solution (PES; ⁴⁹) with a light:dark cycle of 12:12 h (20 µmol photons m⁻² s⁻¹) using daylight-type fluorescent tubes ⁵⁰. All manipulations were performed under a laminar flow hood in sterile conditions. Immature gametophytes (i.e., absence of sex-specific reproductive structures, oogonia or antheridia) of each strain were frozen in liquid nitrogen and kept at -80C until RNA extraction.

392

393

394

395

396

397

398

RNA from male and female pools was extracted from triplicate samples (each containing at least 800 individual gametophytes; except for *S. polyschides* and *S. dermatodea* where two replicates were used) using Qiagen RNA extraction Plant Mini kit. RNA quality and quantity were assessed using an Agilent 2100 bioanalyzer, associated with an RNA 6000 Nano kit. For each replicate, the RNA was quantified and cDNA was synthesized using an oligo-dT primer. The cDNA was fragmented, cloned, and sequenced by Fasteris (CH-1228 Plan-les-Ouates, Switzerland), using Illumina Hi-seq2000 for Saccorhiza and Desmarestia species; by Genome Quebec using an or Nextgen6000 for Halopteris and Chordaria species; and by Genoscope using Illumina Hi-seq 4000 for Sphacelaria and Sphaerotrichia species (see Supplementary Table 1 for details).

399

400

401

402

403

404

405

406

407

## TRANSCRIPTOME ASSEMBLIES AND GENE SET PREDICTIONS

408

Predicted gene sets were constructed for each species base on genome and transcriptome assemblies. In order to filter out potential contamination, first round assembled contigs were blasted against the NCBI non-redundant (nr) protein database using diamond v 0.9.21 ⁵¹ and reads that mapped on contigs with non-eukaryotic taxons were removed using blobtools v 1.0.1 ⁵². De novo transcriptomes were assembled using Trinity (*Saccorhiza polyschides*, *Saccorhiza dermatodea*, *Desmarestia dudresnayi*, *Desmarestia*

409

410

411

412
413 *herbacea* female, *Halopteris paniculata*, *Sphacelaria rigidula*) or rnaSPADES v 3.12.0 (*Chordaria linearis*,
414 *Sphaerotrichia firma*) with kmer size of 55.

415 All genomes were soft-masked using Repeatmasker v 4.0.9 after building a de novo transposable
416 elements and repeats database with RepeatModeler v 1.0.8 53. BRAKER2 54 and PASA (for *Desmarestia*
417 *herbacea* 55, using input predicted protein from the reference species *Ectocarpus* sp.
418 (EctsiV2_prot_LATEST.tfa 56, were used to predict gene sets used for all downstream analyses.

419 The final assemblies are available in NCBI (BioProject accession number PRJNA733856). Transcriptome
420 completeness was assessed using BUSCO v3 eukaryote gene set as reference (Odb9). Transcripts that had
421 DNA data support for only one sex (potentially sex-linked) were tested with PCR using at least 4 males and
422 4 females per species and removed from the sex-biased gene analysis. PCR primers are detailed in
423 Supplementary Table 17.

424 EXPRESSION QUANTIFICATION AND INFERENCE OF SEX-BIASED GENES

425 RNAseq reads adaptors were trimmed using trimmomatic v0.38 57 which was also used for read-quality
426 filtering: reads were removed if the leading or trailing base had a Phred score <3, or if the sliding window
427 Phred score, averaged over four bases, was <15. Reads shorter than 36 bases were discarded (as well as
428 pair of reads, if one of the pair was <36 bases long). Trimmomatic-processed RNAseq reads from each
429 library were used to quantify gene expression with kallisto v 0.44.0 58 using 31 bp-long kmers and
430 predicted transcript of each species. RNAseq libraries were composed of stranded (--fr-stranded or --rf-
431 stranded option) single-end reads (--single option) or paired-end reads (Supplementary Table 1). A gene
432 was considered expressed in a given species and/or a given sex when at least one library displayed an
433 expression level (in TPM) above the 5th percentile of TPM distribution across all genes and libraries within
434 a species and sex. Following 59, transcript abundances were then summed up within genes and multiplied
435 by the total library size, using the tximport package 35 to obtain the expression level for each gene in
436 transcripts per million reads (TPM).

437 Estimates of sex-biased gene expression in dioicous species were obtained using read count matrices as
438 input for the DESeq2 package (Love et al., 2014) in R 3.6.3. *P*-values were corrected for multiple testing
439 using Benjamini and Hochberg's algorithm in DESeq2, applying an adjusted *P*-value cut-off of 0.05 for
440 differential expression analysis. In addition, only genes with a minimum of 2-fold change expression level
441 between sexes were retained as sex-biased.

442 QUANTIFICATION OF PHENOTYPIC SEXUAL DIMORPHISM

443 Individual gametophytes from each of the strains were isolated in sea water and observed using an
444 inverted transmitted light microscope DMi8 (Leica) with the LAS X software. Between 269 and 556 cells
445 (348 cells on average per sex and per species) across five different gametophytes per species were
446 individually measured using Fidji 60. We used *t*-tests to compare cell length between groups. The
447 difference in mean cell length between sexes of dioicous species was computed and used as a proxy for
448 phenotypic sexual dimorphism. To investigate the relationship between phenotypic sexual dimorphism
449 and extent in sex-biased expression, phenotypic dimorphism was regressed against the fraction of sex-
450 biased genes in R, corrected for phylogeny using pylogenetic generalized least square method as
451 implemented in the nlme R package⁶¹.

452 ORTHOLOGY AND EVOLUTIONARY RATES WITHIN SPECIES PAIRS

453 We inferred pairwise single-copy orthologs (PSO) within the four species pair using Orthofinder with
454 default parameters⁶². We used kallisto v 0.44.0 to quantify expression levels for PSO within species pairs.

455 In order to infer the potential role of selection in expression changes between dioicous and monoicous
456 species we computed Δₓ. To summarize, we calculated Δₓ = *d* / *r* with *d* and *r* respectively given by:

457 $d=\text{Mean } X_{\text{dioicous}}-\text{Mean } X_{\text{monoicous}}/\text{Mean } X_{\text{dioicous}}$

458 and

459 $r=[X_{\text{dioicous}}]^{\text{high}}-[X_{\text{dioicous}}]^{\text{low}}/\text{Mean } X_{\text{dioicous}}$

460 where X is the expression level measured in TPM, 'High' and 'Low' represent the maximum and minimum
461 values. Δₓ was computed separately for females and males of the dioicous species, and for male-biased
462 genes, female-biased genes and unbiased genes also separately. Orthologs with |Δx| > 1 and a minimum
463 expression fold-change between sister species of 1.5 were considered to have had a significant
464 evolutionary expression shift. Fisher exact tests were computed to detect whether female-biased genes
465 (FBG) and male-biased genes (MBG) were more likely to show an absolute value of Δₓ > 1 compared to
466 unbiased genes.
467

467 Orthologous proteins between species pairs were aligned with MAFFT v7.453 63 the alignments were
468 curated with Gblocks v0.91b 64 and back-translated to nucleotides using translatorX 65. We used these
469 nucleotide alignments as input for phylogenetic analysis by maximum likelihood (PAML4, CodeML, 66) to
470 infer pairwise dN/dS (ω) with F3x4 model of codon frequencies. We retained orthologs with 0 < dS < 2 as
471 valid for further analysis. We compared species and sexes evolutionary rates separately for female-
472 biased, male-biased and unbiased genes, using permutation t-tests in R with 100,000 permutations.
473

473 EVOLUTION OF SEX-BIASED GENE EXPRESSION

474 We inferred a single orthologous gene set for the four dioicous species (DSO) using Orthofinder with
475 default parameters. Following the methods used in⁶⁷ we included in the DSOs the orthogroups genes that
476 were 1:1:1:0, likely due to situations in which a single-copy ancestral gene was lost in a single species. To
477 further account for gene prediction errors, we also included orthogroups with a single species presenting
478 two-genes that aligned on more than 60% of their length as duplicate genes. In the latter case, the
479 longest duplicated sequence was retained for further analysis.
480

480 A well resolved phylogeny of the Pheaophyceae was used as reference gene tree ³⁴ to infer where sex-
481 biased gene expression evolved along the phylogenetic tree. We coded DSO as either male-biased,
482 female-biased or unbiased for each species and used the ape package 68 in R to reconstruct the discrete
483 ancestral state *via* maximum likelihood. Proportions of ancestral genes in each category were plotted as
484 pie-charts on tree nodes and gain/loss of bias were reported on each branch. We further tested the
485 significance of overlap between sex-biased genes identified within dioicous species with exact multi-set
486 intersection test implemented in the SuperExactTest package v 1.0.7 in R⁶⁹.

487 We computed absolute standardized phylogenetic independent contrasts (PICs) among dioicous species,
488 using the *ape* package in R. Mean PICs were compared using Mann-Whitney Rank tests, between
489 unbiased genes and sex-biased genes with their expression measured in species in which they were not
490 sex-biased.

491 We inferred expression profile similarity index between monoicous species and males and females of
492 dioicous species within pairs as the Pearson correlation coefficient of PSO expression levels in log2(TPM +
493 1). This analysis was performed for all expressed genes, and separately for MB, FB and unbiased genes.
494 We compared Pearson coefficients of regression within each species pair, using the cocor package ⁴⁰,
495 considering gene expression profiles of males and females as independent gene sets. We also compared
496 SBG with unbiased genes within sexes, considering these gene sets as dependent. We report the *P*-value
497 based on Fisher's *z* or, when possible, Silver, Hittner and May's modification of Dunn and Clark's *z*.
498 Pearson's coefficients were plotted for each species pair.

499 CONVERGENT EXPRESSION CHANGES

500 Convergent changes associated with transitions to monoicy were investigated on single-copy orthologs
501 inferred across the eight studied species (termed 'All Single-copy Orthologs', ASO) following the same
502 methods as those used for the DSOs. Using this data set, we quantified gene expression with kallisto as
503 described above, and DESeq2 was used to infer orthologs significantly affected by sexual system but not
504 species pair (IfcShrink with "ashr" method, sexual system contrast ⁷⁰). Significance of the number of
505 convergent expression changes was tested with permutation tests (100,000 permutations). We used the
506 ComplexHeatmap package in R to visualize gene expression for each replicate. Orthogroups with
507 inconsistent sex-bias across different species (n=139) were removed from the dN/dS analysis of
508 convergent gene evolution.

509 Intersects between genes across PSO, DSO and ASO were represented using the UpSetR package v1.4.0
510
511

511 ASO EVOLUTIONARY RATES

512 Following the same process described for pairwise orthologs, we aligned and studied molecular
513 sequence divergence for all species orthologs (ASO) using CodeML. We used a 'two-ratio' branch
514 model (model = 2, Nssites = 0) to specifically study divergence on monoicous branches (foreground
515 branches). We compared ω values separately between sex-biased (male-biased and female-biased
516 genes) and unbiased genes with permutation t-tests (10,000 permutations). We also ran two
517 branch-site models in PAML to detect positive selection in foreground branches (model=2, Nsite=2,
518 ω=1 fixed (Null) or allowed to vary). Likelihood-ratio tests were used to compare the model of
519 selection with the null model in order to detect orthologs with sites under positive selection in the
520 monoicous branches. LRT P-values were corrected for multiple testing using Benjamini and
521 Hochberg's algorithm ⁷².
522 FUNCTIONAL ANNOTATION ANALYSIS

523 Predicted genes and orthogroups were blasted against the NCBI non-redundant (nr) protein database
524 with blast (v2.9.0). Functional annotation was performed using BLAST2GO 38, as well as the InterProScan
525 prediction of putative conserved protein domains 73. Gene set enrichment analysis was carried out
526 separately for each gene set using Fisher's exact test implemented in the TopGO package, with the
527 weight01 algorithm 74. Values were corrected for multiple testing using the Benjamini-Hochberg method
528 in order to control the false discovery rate. We investigated enrichment in terms of biological process
529 ontology and reported significant GO-terms with P-value < 0.01. All statistical analyses were performed in
530 R 3.6.3, plots were produced with ggplot2 in R (https://ggplot2.tidyverse.org/).

531 ACKNOWLEDGMENTS

532 We would like to thank J.R. Pannell for valuable comments on the present study, I. Theodorou for help
533 with microscopy, A. F. Peters for help with algal cultures and T. Broquet and H. Drost for comments and
534 advice on statistical methods. We are grateful to the Institut Français de Bioinformatique (ANR-11-INBS-
535 0013), BioGenouest and the Roscoff Bioinformatics platform ABiMS (http://abims.sb-roscoff.fr) for
536 providing computing and storage resources. This work was supported by an ERC grant to S.M.C (grant
537 agreement 864038), the France Génomique National Infrastructure project Phaeoexplorer (ANR-10-INBS-
538 09), the Centre National de la Recherche Scientifique and Sorbonne University.

539 DATA

540 Raw reads have been deposited in the SRA. BioProject accession number is PRJNA733856. Accession
541 codes are given in Supplementary Table 18.

542 AUTHORS CONTRIBUTIONS

543 S.M.C designed the study. Z.C conducted DNA and RNA extractions. C.C. prepared sequencing libraries.
544 O.G assembled transcriptomic gene sets of *S. rigidula*, *S. firma*, *C. linearis* and *H. paniculata*. A.L
545 assembled transcriptomic gene sets and identified V-linked genes of *D. herbacea*. G.G.C. assembled
546 transcriptomic gene sets of *D. dudresnayi*, *S. polyschides* and *S. dermatodea*. G.G.C conducted
547 bioinformatic and statistical analyses. S.M.C. and J.M.C contributed with resources. G.C. and S.M.C wrote
548 the manuscript. All authors approved the final version of the manuscript.

549 COMPETING INTERESTS

550 The authors declare no potential competing interests.

551 FIGURE LEGENDS

552 Figure 1. Diagram of the phylogeny of the eight species of brown algae investigated. Approximate
553 estimated age of nodes is based on³⁴; O. de Clerck pers. communication). A schematic view of typical
554 gamete size differences (female in red, male in blue) per species pair is presented. Dioicous species (D)
555 are marked in brown and monoicous species (M) in black.
Figure 2. Patterns of sexual dimorphism in dioicous brown algae. A) Pie charts representing the fractions of sex-biased genes among expressed genes (female-bias in red, male-bias in blue) in the four dioicous species. Gradients of colors represent the intensity of expression fold-change (FC), from 2FC difference to more than 15FC. The percentages are calculated based on the total number of expressed genes averaged across sexes. B) Comparison of gene expression levels, in log2(TPM+1), between males and females within dioicous species. Colour patterns follow the ones used in panel A, except for grey points which represent unbiased genes that presented a FC > 2. C) Scatterplots of the lengths of cells of immature gametophytes of dioicous species. The mean (point) and standard deviations (whiskers) are plotted per sex per species. Stars indicate significant difference between mean cell length, tested with two-sided t-tests. *0.01 < P < 0.05; **0.001 < P < 0.01; ***P < 0.001. D) Representative micrographs of male and female immature gametophytes viewed under an inverted light microscope for each dioicous species investigated. Micrographs show individual algae, representative of between 100-200 individuals grown in petri dishes. E) Linear regressions of the fraction of female- and male-biased genes (in red and blue, respectively) among the mean number of expressed genes across both sexes, against the mean difference in cell length recorded between the sexes (in µm), in the four dioicous species investigated. Linear regressions were fitted through phylogenetic generalized least square method, implemented in the R package "nlme". We report values of adjusted r² calculated with ANOVA.

Figure 3. Reconstruction of ancestral sex-biased gene sets across the four dioicous species. The number of inferred sex-biased genes (female-bias in red, male-bias in blue) at ancestral nodes as well as the inferred gain and loss of sex-biased genes along branches are displayed.

Figure 4. Evolution of sex-biased genes during transitions to monoicy. A) Comparison of gene expression levels within species pairs, in log2(TPM+1), using PSO gene sets. F: females. M: males. Mo: monoicous. The number of female-biased (FBG) and male-biased genes (MBG) among PSO are displayed. Note that only the sex-biased genes with a single-copy ortholog in the corresponding monoicous specie are displayed in the plots (in other words, the SBG represented in the plots are a subset of the SBG identified within each dioicous species). Boxes represent the interquartile range (25th and 75th percentiles) of the data, the line inside the box represents the median, whiskers represent the largest/smallest value within 1.5 times interquartile range above and below the 75th and 25th percentile, respectively. Statistical tests are permutation t-tests using 100,000 permutations. Paired two-sided t-test were used for comparisons between sexes of the same species (dioicous species). B) Comparisons of similarity index values (Pearson coefficients) between expression profiles (in log2(TPM+1)) of pairwise single copy orthologs (PSOs) between monoicous and dioicous species pairs. The figure represents male versus female similarity indexes in relation to the monoicous expression profiles. Note that similarity index are represented separately for sex-biased genes in females (red) and in males (blue), as well as for unbiased genes averaged across sexes in the dioicous species (black). Pearson coefficients were plotted for each species pair in increasing order of the proportion of SBG among expressed genes of dioicous species (x axis). Stars in the top panel represent significant differences between Pearson coefficients, taking into account the correlations between compared gene sets, using the cocor package in R. Red and blue stars indicate a significant difference between female (red stars) or male sex-biased genes (blue stars) Pearson coefficients with unbiased genes coefficients. Top panel black stars indicate a significant difference of
Pearson coefficients of unbiased genes between males and females. Significant differences of coefficients between sex-biased genes in females and males are indicated directly on the plot. * : 0.01 < P < 0.05; ** : 0.001 < P < 0.01; *** : P < 0.001. See also Figure S4 and S5. C) Fraction of female-biased genes (red), male-biased genes (blue) and non-biased genes (grey) with an absolute value of ΔX >1 and a fold-change > 1.5, calculated within species pairs (on PSO). The percentages are calculated on the total number of orthologs in each category. Down-regulated genes in the monoicous species are represented below the y=0 line (ΔX < -1), upregulated genes in the monoicous species are represented above the y=0 line (ΔX >1). Stars indicate a significant over-representation of female-biased or male-biased genes with an absolute ΔX >1 compared with the proportion of unbiased genes with ΔX >1, tested using Fisher exact tests. D) Sequence divergence, measured as dN/dS (ω), between dioicous and monoicous species calculated within species pairs (PSO). Boxes represent the interquartile range (25th and 75th percentiles) of the data, the line inside the box represents the median, whiskers represent the largest/smallest value within 1.5 times interquartile range above and below the 75th and 25th percentiles, respectively. Statistical tests are permutation two-sided t-tests using 100,000 permutations, p-values are displayed in parentheses. * : 0.01 < P < 0.05; ** : 0.001 < P < 0.01; *** : P < 0.001. E) Principal component analysis (PCA) plot of all the RNA-seq samples, using ASOs. Monoicous species are plotted in orange, female samples in red and male samples in blue.

613
614 SUPPLEMENTARY TABLE LEGENDS

615 Supplementary Table 1. Brown algae species used in the study and summary of gene expression data sets.

616 Supplementary Table 2. Summary statistics of sex-biased gene expression. Fractions of sex-biased genes
617 are calculated over the number of expressed genes in the sex of bias.

618 Supplementary Table 3. Gene expression values in transcripts per million (TPM; from kallisto) and
619 expression bias (from Deseq2) for each of the genes across the four dioicous species.

620 Supplementary Table 4. Summary of sexually dimorphic traits in the four dioicous brown algae species
621 investigated in the present study.

622 Supplementary Table 5. Orthogroups belonging to the dioicous single-copy orthologs gene set (DSO) and
623 their corresponding gene and bias per species.

624 Supplementary Table 6. Description of GO-term enrichment of sex-biased genes identified in each
625 dioicous species. *P*-values correspond to Fisher exact tests (FDR < 0.05). Terms found in two or more
626 species are highlighted in bold.

627 Supplementary Table 7. Summary description of all pairwise single-copy orthologs (PSO).

628 Supplementary Table 8. Orthogroups belonging to the pairwise single-copy orthologs genes et (PSO)
629 within the Tilopteridales species pair. Mean expression level in TPM across replicates as well as the bias
630 status in the dioicous species are reported.

631 Supplementary Table 9. Orthogroups belonging to the pairwise single-copy orthologs gene set (PSO)
632 within the Desmarestiales species pair. Mean expression level in TPM across replicates as well as the bias
633 status in the dioicous species are reported.

634 Supplementary Table 10. Orthogroups belonging to the pairwise single-copy orthologs gene set (PSO)
635 within the Sphacelariales species pair. Mean expression level in TPM across replicates as well as the bias
636 status in the dioicous species are reported.

637 Supplementary Table 11. Orthogroups belonging to the pairwise single-copy orthologs gene set (PSO)
638 within the Ectocarpales species pair. Mean expression level in TPM across replicates as well as the bias
639 status in the dioicous species are reported.

640 Supplementary Table 12. Summary statistics of Δₓ within the four species pairs. Male-biased genes (MBG)
641 and female-biased genes (FBG) that were significantly more likely or less likely to present |Δₓ| > 1 were
642 highlighted in green and purple, respectively (Fisher's exact tests).

643 Supplementary Table 13. *P*-values of permutation *t*-tests (100,000 permutations) of sequence divergence
644 data (dN/dS) calculated within species pair, between female-, male-biased and unbiased genes.

645 Supplementary Table 14. Summary description of all single-copy orthologs (ASO)
646 Supplementary Table 15. Orthogroups belonging to the all single-copy orthologs gene set (ASO). Mean
647 expression level in TPM across replicates as well as the bias status in the dioicous species are reported.

648 Supplementary Table 16. *P*-values of permutation *t*-tests (10,000 permutations) of sequence divergence
649 data (dN/dS), calculated specifically for monoicous branches (branch model) across ASOs, between
650 female-, male-biased and unbiased genes. Significant difference of divergence with unbiased genes are
651 put in bold.

652 Supplementary Table 17. Primers used to test candidate sex-linked contigs in the different brown algal
653 species.

654 Supplementary Table 18. Accession references.

655

656 REFERENCES

657 1. Charlesworth, D. Theories of the evolution of dioecy. in *Gender and Sexual Dimorphism in Flowering
658 Plants* (eds. Geber, M. A., Dawson, T. E. & Delph, L. F.) 33–60 (Springer Berlin Heidelberg, 1999).
659 doi:10.1007/978-3-662-03908-3_2.

660 2. Charlesworth, B. & Charlesworth, D. A Model for the evolution of dioecy and gynodioecy. *Am. Nat.*
661 **112**, 975–997 (1978).

662 3. Charlesworth, D. Evolution of plant breeding systems. *Curr. Biol. CB 16*, R726-735 (2006).

663 4. Barrett, S. C. H. Understanding plant reproductive diversity. *Philos. Trans. R. Soc. Lond. B. Biol. Sci.*
664 **365**, 99–109 (2010).

665 5. Charnov, E., Bull, J. J. & Maynard Smith, J. Why be an hermaphrodite? *Nature* **263**, 125–126 (1976).

666 6. Meagher, T. R. Linking the evolution of gender variation to floral development. *Ann. Bot.* **100**, 165–
667 176 (2007).

668 7. Ghiselin, M. T. The evolution of hermaphroditism among animals. *Q. Rev. Biol.* **44**, 189–208 (1969).

669 8. Bachtrog, D. *et al*. Are all sex chromosomes created equal? *Trends Genet.* **TIG 27**, 350–357 (2011).

670 9. Mank, J. E. & Ellegren, H. Are sex-biased genes more dispensable? *Biol. Lett.* **5**, 409–412 (2009).

671 10. Parsch, J. & Ellegren, H. The evolutionary causes and consequences of sex-biased gene expression.
672 *Nat. Rev. Genet.* **14**, 83–87 (2013).
673 11. Zemp, N. et al. Evolution of sex-biased gene expression in a dioecious plant. *Nat. Plants* 2, 16168 (2016).

674

675 12. Charlesworth, D. Does sexual dimorphism in plants promote sex chromosome evolution? *Environ. Exp. Bot.* **146**, (2017).

676

677 13. Cossard, G. G., Toups, M. A. & Pannell, J. R. Sexual dimorphism and rapid turnover in gene expression in pre-reproductive seedlings of a dioecious herb. *Ann. Bot.* mcy183–mcy183 (2018) doi:10.1093/aob/mcy183.

678

679

680 14. Harrison, P. W. et al. Sexual selection drives evolution and rapid turnover of male gene expression. *Proc. Natl. Acad. Sci. U. S. A.* **112**, 4393–4398 (2015).

681

682 15. Lipinska et al. Sexual dimorphism and the evolution of sex-biased gene expression in the brown alga *Ectocarpus*. *Mol. Biol. Evol.* **32**, 1581–1597 (2015).

683

684 16. Avise, J. C. & Mank, J. E. Evolutionary perspectives on hermaphroditism in fishes. *Sex. Dev. Genet.* *Mol. Biol. Evol. Endocrinol. Embryol. Pathol. Sex Determ. Differ.* **3**, 152–163 (2009).

685

686 17. Denver, D. R., Clark, K. A. & Raboin, M. J. Reproductive mode evolution in nematodes: insights from molecular phylogenies and recently discovered species. *Mol. Phylogenet. Evol.* **61**, 584–592 (2011).

687

688 18. Lloyd, D. G. The maintenance of gynodioecy and androdioecy in angiosperms. *Genetica* **45**, 325–339 (1975).

689

690 19. Pannell, J. Mixed genetic and environmental sex determination in an androdioecious population of *Mercurialis annua*. *Heredity* **78**, 50–56 (1997).

691

692 20. Schaefer, H. & Renner, S. S. A three-genome phylogeny of Momordica (Cucurbitaceae) suggests seven returns from dioecy to monoecy and recent long-distance dispersal to Asia. *Mol. Phylogenet. Evol.* **54**, 553–560 (2010).

693

694

695 21. Cossard, G. G., Gerchen, J. F., Li, X., Cuenot, Y. & Pannell, J. R. The rapid dissolution of dioecy by experimental evolution. *Curr. Biol. CB* (2021) doi:10.1016/j.cub.2020.12.028.

696
697 22. Kafer, J., Marais, G. A. B. & Pannell, J. R. On the rarity of dioecy in flowering plants. *Mol. Ecol.* **26**, 1225–1241 (2017).

699 23. Crossman, A. & Charlesworth, D. Breakdown of dioecy: models where males acquire cosexual functions. *Evol. Int. J. Org. Evol.* **68**, 426–440 (2014).

701 24. Ehlers, B. K. & Bataillon, T. ‘Inconstant males’ and the maintenance of labile sex expression in subdioecious plants. *New Phytol.* **174**, 194–211 (2007).

703 25. Umen, J. & Coelho, S. Algal sex determination and the evolution of anisogamy. *Annu. Rev. Microbiol.* **73**, 267–291 (2019).

705 26. Coelho, S. M., Gueno, J., Lipinska, A. P., Cock, J. M. & Umen, J. G. UV chromosomes and haploid sexual systems. *Trends Plant Sci.* **23**, 794–807 (2018).

707 27. Bachtrog, D. et al. Sex determination: why so many ways of doing it? *PLoS Biol.* **12**, e1001899 (2014).

708 28. Mignerot, L. & Coelho, S. M. The origin and evolution of the sexes: Novel insights from a distant eukaryotic linage. *C. R. Biol.* **339**, 252–257 (2016).

710 29. Luthringer, R., Cormier, A., Peters, A. F., Cock, J. M. & Coelho, S. M. Sexual dimorphism in the brown algae. *Perspectives in Phycology* **1**, 11–25 (2015).

712 30. Villarreal, J. C. & Renner, S. S. Correlates of monoicy and dioicy in hornworts, the apparent sister group to vascular plants. *BMC Evol. Biol.* **13**, 239 (2013).

714 31. Bringloe, T. T. et al. Phylogeny and evolution of the brown algae. *Crit. Rev. Plant Sci.* **39**, 281–321 (2020).

716 32. Heesch, S. et al. Evolution of life cycles and reproductive traits: insights from the brown algae. *bioRxiv* (2019) doi:10.1101/530477.

718 33. Cánovas, F. G., Mota, C. F., Serrão, E. A. & Pearson, G. A. Driving south: a multi-gene phylogeny of the brown algal family Fucaceae reveals relationships and recent drivers of a marine radiation. *BMC Evol. Biol.* **11**, 371 (2011).
34. Kawai, H., Hanyuda, T., Draisma, S. G. A., Wilce, R. T. & Andersen, R. A. Molecular phylogeny of two unusual brown algae, *Phaeostrophion irregulare* and *Platysiphon glacialis*, proposal of the Stschapoviales ord. nov. and Platysiphonaceae fam. nov., and a re-examination of divergence times for brown algal orders. *J. Phycol.* **51**, 918–928 (2015).

35. Love, M. I., Huber, W. & Anders, S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. *Genome Biol.* **15**, 550 (2014).

36. Scharmann, M., Rebelo, A. G. & Pannell, J. R. High rates of evolution preceded shifts to sex-biased gene expression in *Leucadendron*, the most sexually dimorphic angiosperms. *bioRxiv* 2021.01.12.426328 (2021) doi:10.1101/2021.01.12.426328.

37. Pointer, M. A., Harrison, P. W., Wright, A. E. & Mank, J. E. Masculinization of gene expression Is associated with exaggeration of male sexual dimorphism. *PLOS Genet.* **9**, e1003697 (2013).

38. Conesa, A. & Götz, S. Blast2GO: A comprehensive suite for functional analysis in plant genomics. *Int. J. Plant Genomics* **2008**, 619832 (2008).

39. Papakostas, S. et al. Gene pleiotropy constrains gene expression changes in fish adapted to different thermal conditions. *Nat. Commun.* **5**, 4071 (2014).

40. Diedenhofen, B. & Musch, J. cocor: a comprehensive solution for the statistical comparison of correlations. *PloS One* **10**, e0121945 (2015).

41. Ometto, L., Shoemaker, D., Ross, K. G. & Keller, L. Evolution of Gene Expression in Fire Ants: The Effects of Developmental Stage, Caste, and Species. *Mol. Biol. Evol.* **28**, 1381–1392 (2010).

42. Evans, J. P. & Lymbery, R. A. Sexual selection after gamete release in broadcast spawning invertebrates. *Philos. Trans. R. Soc. B Biol. Sci.* **375**, 20200069 (2020).

43. Hoshino, M., Okino, T. & Kogame, K. Parthenogenetic female populations in the brown alga *Scytosiphon lomentaria* (Scytosiphonaceae, Ectocarpales): decay of a sexual trait and acquisition of asexual traits. *J. Phycol.* (2018) doi:10.1111/jpy.12812.
745 44. Lipinska *et al*. Development of PCR-based markers to determine the sex of kelps. *PLoS ONE* **10**, e0140535 (2015).

747 45. Lipinska, A. P. *et al*. Multiple gene movements into and out of haploid sex chromosomes. *Genome Biol.* **18**, 104 (2017).

749 46. Duret, L. & Mouchiroud, D. Determinants of substitution rates in mammalian genes: expression pattern affects selection intensity but not mutation rate. *Mol. Biol. Evol.* **17**, 68–74 (2000).

751 47. Kimura, M. & Ohta, T. On some principles governing molecular evolution. *Proc. Natl. Acad. Sci. U. S. A.* **71**, 2848–2852 (1974).

753 48. Yamamoto, K. *et al*. Three genomes in the algal genus Volvox reveal the fate of a haploid sex-determining region after a transition to homothallism. *Proc. Natl. Acad. Sci. U. S. A.* **118**, (2021).

755 49. Starr, R. & Zeikus, J. UTEX—The Culture Collection of Algae at the University of Texas at Austin. *J. Phycol.* **29**, 1–106 (2004).

757 50. Coelho, S. M. *et al*. How to cultivate *Ectocarpus*. *Cold Spring Harb. Protoc.* **2012**, 258–261 (2012).

758 51. Buchfink, B., Xie, C. & Huson, D. H. Fast and sensitive protein alignment using DIAMOND. *Nat. Methods* **12**, 59–60 (2015).

760 52. Laetsch, D. & Blaxter, M. BlobTools: Interrogation of genome assemblies [version 1; referees: 2 approved with reservations]. *F1000Research* **2017**, 61287 6:1287 (2017).

762 53. Flynn, J. M. *et al*. RepeatModeler2 for automated genomic discovery of transposable element families. *Proc. Natl. Acad. Sci. U. S. A.* **117**, 9451–9457 (2020).

764 54. Brůna, T., Hoff, K. J., Lomsadze, A., Stanke, M. & Borodovsky, M. BRAKER2: automatic eukaryotic genome annotation with GeneMark-EP+ and AUGUSTUS supported by a protein database. *NAR Genomics Bioinforma.* **3**, lqaa108 (2021).

767 55. Haas, B. J. *et al*. Improving the Arabidopsis genome annotation using maximal transcript alignment assemblies. *Nucleic Acids Res.* **31**, 5654–5666 (2003).
769 56. Cormier, A. et al. Re-annotation, improved large-scale assembly and establishment of a catalogue of
770 noncoding loci for the genome of the model brown alga *Ectocarpus*. *New Phytol.* **214**, 219–232
771 (2017).

772 57. Bolger, A. M., Lohse, M. & Usadel, B. Trimmomatic: a flexible trimmer for Illumina sequence data.
773 *Bioinforma. Oxf. Engl.* **30**, 2114–2120 (2014).

774 58. Bray, N. L., Pimentel, H., Melsted, P. & Pachter, L. Near-optimal probabilistic RNA-seq quantification.
775 *Nat. Biotechnol.* **34**, 525–527 (2016).

776 59. Soneson, C. & Robinson, M. D. Bias, robustness and scalability in single-cell differential expression
777 analysis. *Nat. Methods* **15**, 255–261 (2018).

778 60. Schindelin, J. et al. Fiji: an open-source platform for biological-image analysis. *Nat. Methods* **9**, 676–
779 682 (2012).

780 61. Pinheiro, J., Bates, D., DebRoy, S. S. & Sarkar, D. NLme: linear and nonlinear mixed effects models. *R
781 Package Version* **31-110 3**, 1–113 (2013).

782 62. Emms, D. M. & Kelly, S. OrthoFinder: solving fundamental biases in whole genome comparisons
783 dramatically improves orthogroup inference accuracy. *Genome Biol.* **16**, 157 (2015).

784 63. Katoh, K. & Standley, D. M. MAFFT Multiple Sequence Alignment Software Version 7: Improvements
785 in Performance and Usability. *Mol. Biol. Evol.* **30**, 772–780 (2013).

786 64. Talavera, G. & Castresana, J. Improvement of phylogenies after removing divergent and ambiguously
787 aligned blocks from protein sequence alignments. *Syst. Biol.* **56**, 564–577 (2007).

788 65. Abascal, F., Zardoya, R. & Telford, M. J. TranslatorX: multiple alignment of nucleotide sequences
789 guided by amino acid translations. *Nucleic Acids Res.* **38**, W7-13 (2010).

790 66. Yang, Z. PAML 4: phylogenetic analysis by maximum likelihood. *Mol. Biol. Evol.* **24**, 1586–1591 (2007).

791 67. Lipinska, A. P. et al. Rapid turnover of life-cycle-related genes in the brown algae. *Genome Biol.* **20**, 35
792 (2019).
793 68. Paradis, E. & Schliep, K. ape 5.0: an environment for modern phylogenetics and evolutionary analyses in R. *Bioinformatics* **35**, 526–528 (2019).

795 69. Wang, M., Zhao, Y. & Zhang, B. Efficient Test and Visualization of Multi-Set Intersections. *Sci. Rep.* **5**, 16923 (2015).

797 70. Stephens, M. False discovery rates: a new deal. *Biostatistics* **18**, 275–294 (2017).

798 71. Conway, J. R., Lex, A. & Gehlenborg, N. UpSetR: an R package for the visualization of intersecting sets and their properties. *Bioinforma. Oxf. Engl.* **33**, 2938–2940 (2017).

800 72. Benjamini, Y. & Hochberg, Y. On the Adaptive Control of the False Discovery Rate in Multiple Testing with Independent Statistics. *J. Educ. Behav. Stat.* **25**, 60–83 (2000).

802 73. Quevillon, E. et al. InterProScan: protein domains identifier. *Nucleic Acids Res.* **33**, W116-120 (2005).

803 74. Alexa, A. & Rahnenfuhrer. topGO: Enrichment Analysis for Gene Ontology. (2019).

804