# Sex-limited experimental evolution drives transcriptomic divergence in a hermaphrodite

Aivars Cīrulis ID <fnref n="1" />,<fnref n="2" />,<fnref n="3" />, Anna K. Nordén<fnref n="1" />, Allison M. Churcher<fnref n="4" />, Steven A. Ramm<fnref n="5" />,<fnref n="6" />, Kira S. Zadesenets<fnref n="7" />, and Jessica K. Abbott ID <fnref n="1" />,*

## Abstract

The evolution of gonochorism from hermaphroditism is linked with the formation of sex chromosomes, as well as the evolution of sex-biased and sex-specific gene expression to allow both sexes to reach their fitness optimum. There is evidence that sexual selection drives the evolution of male-biased gene expression in particular. However, previous research in this area in animals comes from either theoretical models or comparative studies of already old sex chromosomes. We therefore investigated changes in gene expression under 3 different selection regimes for the simultaneous hermaphrodite *Macrostomum lignano* subjected to sex-limited experimental evolution (i.e. selection for fitness via eggs, sperm, or a control regime allowing both). After 21 and 22 generations of selection for male-specific or female-specific fitness, we characterized changes in whole-organism gene expression. We found that female-selected lines had changed the most in their gene expression. Although annotation for this species is limited, gene ontology term and Kyoto Encyclopedia of Genes and Genomes pathway analyses suggest that metabolic changes (e.g. biosynthesis of amino acids and carbon metabolism) are an important adaptive component. As predicted, we found that the expression of genes previously identified as testis-biased candidates tended to be downregulated in the female-selected lines. We did not find any significant expression differences for previously identified candidates of other sex-specific organs, but this may simply reflect that few transcripts have been characterized in this way. In conclusion, our experiment suggests that changes in testis-biased gene expression are important in the early evolution of sex chromosomes and gonochorism.

**Key words:** *Macrostomum lignano*, experimental evolution, sexual selection, sex-biased gene expression, evolution of go-
nochorism, hermaphrodite, sex chromosome evolution.

<footnote n="1">¹Department of Biology, Lund University, 223 62 Lund, Sweden</footnote>

<footnote n="2">²Laboratory of Microbiology and Pathology, Institute of Food Safety, Animal Health and Environment "BIOR," Riga LV-1076, Latvia</footnote>

<footnote n="3">³Faculty of Biology, University of Latvia, Riga LV-1004, Latvia</footnote>

<footnote n="4">⁴Department of Molecular Biology, National Bioinformatics Infrastructure Sweden, Umeå University, 901 87 Umeå, Sweden</footnote>

<footnote n="5">⁵Department of Evolutionary Biology, Bielefeld University, 33615 Bielefeld, Germany</footnote>

<footnote n="6">⁶UMR 6553 ECOBIO, Université de Rennes, 35042 Rennes, France</footnote>

<footnote n="7">⁷Department of Molecular Genetics, Cell Biology and Bionformatics, The Federal Research Center Institute of Cytology and Genetics SB RAS, 630090 Novosibirsk, Russian Federation</footnote>

<footnote>*Corresponding author: E-mail: jessica.abbott@biol.lu.se.</footnote>

<footnote>Accepted: December 23, 2023</footnote>
## Significance

The evolution of separate sexes from a hermaphroditic ancestor can evolve either gradually, via small changes in sex-specific gene expression, or rapidly, via sterility mutations in either sex role, leading to the evolution of sex chromosomes. However, this has never been observed in real time. We therefore carried out sex-limited experimental evolution in a hermaphroditic worm to study these early stages in the evolution of separate sexes. After 20 generations of evolution toward females or males, we found that worms from the female-limited selection lines have diverged the most in gene expression. The response was mainly by downregulating genes associated with male function, suggesting that sex-specific gene expression can evolve rapidly in transitions from hermaphroditism to separate sexes.

## Introduction

The evolution of separate sexes from a hermaphrodite ancestor either results gradually by decreasing investment in one sex function or rapidly by the evolution of sex chromosomes via sterility mutations (Bachtrog et al. 2014). Sex chromosomes generally arise from a pair of autosomes, which acquire a major sex-determining gene (or genes) by mutation (most commonly the duplication of a downstream gonadal factor undergoing transcriptional rewiring [Rafati et al. 2020]). During the transition from hermaphrodite to separate-sexed populations, at least 2 mutations are believed to be needed; one that causes sterility in the female function (resulting in the production of males) and one that causes sterility in the male function (resulting in the production of females). These mutations may occur in the same or separate genes. It is unlikely that both mutations arise simultaneously, so populations transitioning to separate sexes are expected to go through an intermediate step of androdioecy (males and hermaphrodites) or gynodioecy (females and hermaphrodites), of which gynodioecy is more common (Beukeboom and Perrin 2014). A frequent type of sex-determining mutation is recessive male sterility, which drives the evolution of an XY system if followed by fixation of a dominant female sterility mutation. Conversely, if the initial male sterility mutation is dominant, then a ZW system can evolve, by the acquisition of a recessive female sterility mutation. Once this has happened, beneficial sex-specific genes become linked to the sex-determining region on the Y or W. Recombination is then suppressed to avoid breaking up favorable gene combinations, and heteromorphic sex chromosomes eventually evolve (Beukeboom and Perrin 2014).

Sex chromosomes play an important role in the evolution of sexual dimorphism and facilitate sex-limited or sex-biased gene expression (Grath and Parsch 2016). The degree of sex bias in expression depends on the developmental stage and condition of the organism and often varies across tissues (Grath and Parsch 2016). Unsurprisingly, sex bias in expression is usually strongest in the sexually dimorphic gonads of sexually mature and well-conditioned male and female adults (Grath and Parsch 2016). It is well established that sexual selection shapes the rapid evolution of sex-biased genes and

male-biased genes in particular (Parsch and Ellegren 2013; Veltsos et al. 2017), although genetic drift can also play a role, especially in small populations (Prentout et al. 2023). However, much of what we know of early sex chromosome evolution and the subsequent evolution of sex-biased gene expression in animals is through mathematical models or comparative studies of sex chromosome systems that can be hundreds of thousands to millions of years old (Charlesworth 2019).

The evolution of sexual dimorphism can be investigated in a manipulative way using experimental evolution to observe changes in real time. For example, a number of previous experimental evolution studies have aimed to modify the strength of sexual selection. These studies have been performed in both gonochorists and hermaphrodites, in order to study the genomic and phenotypic response to the altered selection pressure (Hollis et al. 2014; Innocenti et al. 2014; Janicke et al. 2016; Veltsos et al. 2017; Bonel et al. 2018; Wiberg et al. 2021). In one such study in *Drosophila melanogaster*, sexual selection was relaxed in lines with monogamous mating, and gene expression was then compared with control polygamous lines (Hollis et al. 2014). The monogamous lines had the increased expression of female-biased genes and decreased expression of male-biased genes compared to control lines in both sexes, suggesting that sexual selection shapes the expression of sex-biased genes. Another similar study examined the expression profile of female flies only, finding similar changes in gene expression in monogamous females compared with the control lines, which occurred mostly in the ovary (Innocenti et al. 2014). In contrast, data from *D. pseudoobscura* showed masculinization of the transcriptome under monogamous conditions, although the effect was dependent on tissue, sex, and courtship status (Veltsos et al. 2017). Using a different experimental system, Joag et al. (2016) selected males of the bulb mite (*Rhizoglyphus robini*) for 1 of the 2 alternative male morphs in this species (armored fighter males and feminized scrambler males) per selection regime found expression differences in females. These studies demonstrate that altered sex-specific selection pressures have the potential to shape sex differences in expression.

Experimental evolution through sex-specific selection has also been performed in hermaphrodites (Janicke et al. 2016; Bonel et al. 2018). Female-specific selection (via
enforced monogamy) in the hermaphroditic snail *Physa acuta* seems to have resulted in the accumulation of deleterious mutations affecting fitness in the male role as well as juvenile fitness (Bonel et al. 2018). Conversely, male-specific selection (via increased mating group size while equalizing female fitness) did not result in differences compared with the control group (Bonel et al. 2018). These results suggest that sexual selection is important in eliminating deleterious alleles through selection on male fitness, even in hermaphroditic species. While enforced monogamy in *Macrostomum lignano* did not lead to heritable changes in sex allocation measured as testes area against the total area of both gonads, changes in the male copulatory organ (stylet) and sperm morphology were observed (Janicke et al. 2016). These studies show that experimentally manipulating the strength of sexual-specific and sex-specific selection (in particular, relaxing selection on male fitness) can shed light on the evolutionary dynamics of sex-biased genes (or sex function-biased genes in hermaphrodites).

Another approach to understanding the evolution of sex bias in gene expression due to the evolution of sex chromosomes is to carry out sex-limited evolution specifically of sex chromosomes. An experimental evolution study by Abbott et al. (2020) showed that if the X chromosome is inherited only through males in *D. melanogaster*, then the upregulation of male-benefit and the downregulation of female-benefit sexually antagonistic genes evolve. The response to selection involved genes on all chromosomes, many of which were linked to metabolism. However, relatively little is known about how genes acquire sex-specific functions on nascent sex chromosomes since this process is difficult to study in real time. Even young sex chromosome systems are often several million years old (Charlesworth 2019; Rafati et al. 2020; Edvardsen et al. 2022). We therefore aimed to gain a deeper insight into the very earliest stages in the evolution of sex-specific gene expression, using a sex-limited experimental evolution system in the simultaneously hermaphroditic flatworm *M. lignano* (Nordén et al. 2023). More specifically, we have previously subjected 4 replicate populations to male-limited selection, 4 to female-limited selection, and 4 to the control regime in which fitness continued to be gained through both the male and female functions. The selection regimes were achieved using a green fluorescent protein (GFP)-transformed line (Wudarski et al. 2017) in which this genetic marker was used as a proxy for a sex-determining gene to impose sex function-limited evolution. This essentially means that the GFP gene functions as a dominant sterility mutation in one sex function (either a female sterility mutation in the male-limited selection lines or a male sterility mutation in the female-limited selection lines, but without actually causing sterility), mimicking the evolution of a sex-determining locus on an autosome. The male-limited selection lines thus serve as a model for the evolution of a male-dominant sex chromosome system (XY), and the

female-limited selection lines serve as a model for the evolution of a female-dominant sex chromosome system (ZW), meaning that in each generation the GFP-expressing worms from those replicates could gain fitness only through their male or their female sex function, but not both.

A previous study of these lines found evidence of an evolutionary response on the phenotypic level after 14 generations in the female-selected lines, with increased fitness in the female sex function and decreased fitness in the male sex function compared with both the control and male-limited selection lines (Nordén et al. 2023). There were no differences in gonad size between selection regimes, so this suggests that the fitness differences that were observed were not related purely to differences in morphology. Here, we build on these results by investigating the evolutionary response to the selection regime on the genomic level, specifically changes in gene expression. We characterized the whole-organism expression profile after 21 and 22 generations of selection and compared the observed changes between the selection regimes, especially in relation to previously identified sex organ-biased genes (Arbore et al. 2015; Lengerer et al. 2018; Ramm et al. 2019). We predicted that the expression of putative female-specific genes (especially ovary) would increase and expression of putative male-specific genes (especially testes) would decrease in the female-selected lines and vice versa in the male-selected lines. Our results partially support this scenario, and since many of the genes identified from our study had not previously been characterized as having sex function-specific expression or function, we may have uncovered several novel candidates involved in sex function specialization.

## Results

### Differentially Expressed Transcripts and Their Gene Pathways

We were mainly interested in divergence between the female-selected and male-selected lines, and indeed 395 transcripts were differentially expressed between these lines (Fig. 1). However, the largest number of differentially expressed transcripts (1175) was between the female-selected lines and controls (Fig. 1). This was somewhat surprising, since we expected greatest divergence between the male-selected and female-selected lines, but the modest number of differentially expressed transcripts between the male-selected lines and control lines (318) may suggest that these 2 selection regimes are not well differentiated, consistent with previous phenotypic results (Nordén et al. 2023). Finally, we found that changes in 234 transcripts were specific to the female-selected lines and 84 were specific to the male-selected lines as they differ from each other as well as the controls (Fig. 1). There was no evidence that changes in gene expression were observed more often than expected in the GFP scaffold,
<figure><img src="image_1.png" /><figcaption>**Fig. 1.**—Venn diagram for significantly differentially expressed transcripts. The largest number of differentially expressed transcripts was between the female-selected lines and controls (1,175 in total), where transcripts were mainly downregulated in the female-selected lines relative to the control lines. The difference between the male-selected and female-selected lines was smaller (395), and transcripts were mainly upregulated in the male-selected lines relative to the control lines. Finally, the male-selected lines differed from the controls the least (318), and again in this comparison, transcripts are mainly down-regulated in the male-selected lines relative to the control lines. The barplots are a schematic illustration of the most common direction of change between the selection regimes based on the pattern of significant differences and do not reflect actual values.</figcaption></figure>

either when comparing among lines (supplementary fig. S6, Supplementary Material online; paired Wilcoxon signed ranks test MvsF comparison: V = 14,811, P-value = 0.9916; FvsC comparison: V = 16,228, P-value = 0.1616; MvsC comparison: V = 13,299, P-value = 0.237) or in comparison with the rest of the genome (unpaired Wilcoxon rank sum test FvsC: W = 86,084,530, P-value = 0.539; MvsC: W = 84,163,464, P-value = 0.9751), nor was there any evidence that the GFP scaffold was an outlier in terms of its average change in expression; instead, changes have occurred genome wide. Almost all transcripts have changed more than 2-fold in expression, suggesting that these changes are not only biologically but also statistically significant (Fig. 2). We found several classes of transcripts that are particularly interesting in the context of the evolution of sexual dimorphism: (i) 3 transcripts (TRINITY_DN112097_c0_g1_i11, TRINITY_DN114758_c0_g1_i2, and TRINITY_DN102787_c0_g1_i5), which were significantly upregulated in the female-selected lines and significantly downregulated in the male-selected lines, with controls being intermediate and significantly different from both of them, (ii) 234 transcripts, which were significantly differentially expressed in the female-selected lines but do not differ between the control and male-selected lines, and (iii) 84 transcripts, which were significantly differentially expressed in the

male-selected lines but do not differ between the control and female-selected lines (Fig. 1).

To better understand the biological relevance of these differentially expressed transcripts, we performed gene ontology (GO) term and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses. Although annotation for this species is limited, GO term analysis (supplementary table S1, Supplementary Material online) suggests that metabolic changes (e.g. biosynthesis of amino acids and carbon metabolism) may potentially have been an important adaptive component, and this result seems to be confirmed by the KEGG pathway analysis (supplementary table S5, Supplementary Material online).

## Sex Function-Specific Organ Expression

We found changes in some of the 3,360 transcripts previously identified (Ramm et al. 2019) as putatively testis-biased, and these were upregulated in the male-selected lines and downregulated in the female-selected lines more often than expected by chance, although the difference was only marginally significant (Table 1). There was however a robust pattern of downregulation of testis-biased transcripts in the female-selected lines compared with the control lines (Table 1). Additionally, we checked if there was a
<figure><img src="image_2.png" /><figcaption>**Fig. 2.**—Volcano plots of statistically significantly differentially expressed transcripts between FvsC a), MvsC b), and MvsF c). Vertical lines indicate whether the transcripts are also likely to be biologically significant (i.e. have at least a 2-fold change in expression). The very few transcripts, which are not passing this threshold, are colored yellow.</figcaption></figure>

correlation between fold changes in this study and those in Ramm et al. (2019) to see if the changes in our dataset are consistent with a "plasticity-first" scenario, where genes that are plastic in expression in the ancestral population would respond more to the selection regimes. However, the correlation was not significant (supplementary fig. S3, Supplementary Material online), not supporting this hypothesis. There was no difference in expression between the selection regimes in the female-specific cement gland and antrum genes, or male-specific stylet and prostate genes identified by Lengerer et al. (2018) and Weber et al. (2018) (all *P* > 0.05).

A gene known to be sex function-biased that shows a pattern that is consistent with predictions based on function, e.g. a testis-biased transcript that is upregulated in MvsF and MvsC, but downregulated in FvsC, is a perfect candidate for further exploration. We found one such promising candidate gene for female fitness, which corresponds to hypothetical protein BOX15_Mlig029807g1 in the genome assembly published by Wudarski et al. (2017) and annotated as transmembrane protein 72 (human: TMEM72, e = 5e-06) and ovary-biased under transcript ID RNA815_9198 in the transcriptome published by Ramm et al. (2019). In

humans, TMEM72 is located on autosome 10 and is kidney specific, with some expression in the colon and testes as well (https://www.ncbi.nlm.nih.gov/gene/643236). Since it is ovary biased (Ramm et al. 2019) in *M. lignano* and upregulated in the female-selected worms and downregulated in the male-selected worms according to our analysis, this indicates a functional divergence of this protein during animal evolution.

By performing a Principal Component Analysis (PCA) on testes-biased and ovary-biased transcripts, we identified the main transcripts driving differences between the selection regimes (supplementary fig. S2, Supplementary Material online). We observed that Principal Component (PC) 1 for ovaries and PC2 for testes differentiate to some extent the selection lines as expected, where the male-selected and female-selected lines tend to be at the extremes, with control in the middle. The main transcripts driving the difference for testes are TRINITY_DN111916_c0_g1_i10.1, TRINITY_DN92387_c0_g1_i3, and TRINITY_DN87487_c0_g2_i4 and for ovaries areTRINITY_DN110581_c0_g1_i4.1, TRINITY_DN118329_c1_g1_i1, and TRINITY_DN97896_c1_g1_i5 (supplementary fig. S2, Supplementary Material online). However, none of the PCs was found to significantly differ
<caption>**Table 1**</caption>

Sex-specific organ expression

<table><thead><tr><th>Organ</th><th>Upregulated<sup>a</sup></th><th>Downregulated<sup>a</sup></th><th>P-value<sup>b</sup></th></tr></thead><tbody><tr><td colspan="4">Male-selected versus female-selected lines</td></tr><tr><td>Testis (3360)</td><td>13</td><td>4</td><td>0.049</td></tr><tr><td>Ovary (323)</td><td>1</td><td>2</td><td>N/A</td></tr><tr><td>Tail (366)</td><td>0</td><td>0</td><td>N/A</td></tr><tr><td>Gonad (likely ovary-biased) (127)</td><td>2</td><td>1</td><td>N/A</td></tr><tr><td colspan="4">Female-selected versus control lines</td></tr><tr><td>Testis (3360)</td><td>9</td><td>52</td><td>1.803e-08</td></tr><tr><td>Ovary (323)</td><td>2</td><td>1</td><td>N/A</td></tr><tr><td>Tail (366)</td><td>6</td><td>1</td><td>0.125</td></tr><tr><td>Gonad (likely ovary-biased) (127)</td><td>0</td><td>2</td><td>N/A</td></tr><tr><td colspan="4">Male-selected versus control lines</td></tr><tr><td>Testis (3360)</td><td>4</td><td>16</td><td>0.012</td></tr><tr><td>Ovary (323)</td><td>1</td><td>2</td><td>N/A</td></tr><tr><td>Tail (366)</td><td>1</td><td>0</td><td>N/A</td></tr><tr><td>Gonad (likely ovary-biased) (127)</td><td>1</td><td>0</td><td>N/A</td></tr></tbody></table>

between the selection lines (supplementary tables S3 and S4, Supplementary Material online).

## Discussion

Here, we examined the RNA expression profile of a hermaphrodite subjected to sex-limited selection for over 20 generations and found >1,000 differentially expressed transcripts between the selection regimes. Our results show that the male-selected and female-selected lines have both diverged compared with the controls. However, the female-selected lines seem to have responded the most to selection (in terms of the number of transcripts with significantly different expression compared with the controls; Fig. 1), which is in line with results from a previous fitness assay (Nordén et al. 2023). Many of the transcripts were downregulated in the MvsC and FvsC comparisons (i.e. downregulated in M and F relative to C; Fig. 1), and the ones that could be annotated were often linked to metabolic pathways (supplementary tables S5 and table S1, Supplementary Material online). These changes are unlikely to be driven by allometric effects of difference in organ size, since there is no evidence that the lines differ consistently in testes or ovary sizes (Nordén et al. 2023), and results were qualitatively similar when mean organ sizes per line were included in the analysis (data not shown). There is no evidence for

gonochorism of *M. lignano* in natural populations, and the entire genus seems to be exclusively hermaphroditic, so it is unlikely that there are any preexisting loci in this species, which have the potential to act as sex-determining loci. Instead, since it is known that *M. lignano* is highly plastic in terms of morphological response to mating opportunities (Schärer and Ladurner 2003), it is likely that the response that we observed is a result of selection on standing genetic variation in sex allocation and/or plasticity. Unfortunately, the annotation information is limited since *M. lignano* is only distantly related to many standard model organisms. This means that only 40% of transcripts have a putative function (Wasik et al. 2015) and makes comparison with other species challenging. Nevertheless, when we looked at the expression of genes associated with sex-specific organs, then testes-specific transcripts showed more divergence than the ovary-specific transcripts (Table 1), which we discuss in more detail below. Besides selection, of course genetic drift seems to have played a role, since differences in gene expression overall (data not shown) or specifically of organ-biased transcripts (supplementary fig. S2, Supplementary Material online) are very different among the different populations within the selection regime, the result very similar to observations in the plant genus *Silene* (Prentout et al. 2023).

## Differences in Expression Associated With Sex-Specific Organs

Sex-specific gene expression has been shown to vary widely between tissue types and is expected to be most different in the gonads (Grath and Parsch 2016). We (Nordén et al. 2023) and others (Vellnow et al. 2018) have also previously shown that larger testes result in higher male fitness returns. Because of this, we expected to see changes in the expression of putative testes-biased, ovary-biased, gonad-biased, and tail-biased transcripts reported in Ramm et al. (2019). Our results show that several of these transcripts do indeed show differential expression in our study, particularly testes-specific ones. The fact that testes-biased genes were overrepresented and show differences between the selection regimes compared with ovaries (Table 1), even when correcting for the overrepresentation of testes-biased genes compared with the other organs, is consistent with the fact that sexual selection usually has a larger impact on male than female reproductive function (Beukeboom and Perrin 2014). In 2 of the 3 possible pairwise comparisons, testes-associated genes were preferentially upregulated or downregulated in the expected direction (i.e. upregulated in the male-selected lines compared with the female-selected lines and downregulated in the female-selected lines compared with the control lines [Table 1]). Moreover, these results align with preliminary results for sex allocation in worms starved for 2 wk. There, the female-selected lines did not increase testes size when

<footnote><sup>a</sup>Numbers in brackets indicate the total number of transcripts. Transcripts are upregulated or downregulated in the first selection line compared with the second (e.g. 13 testes-biased transcripts are upregulated and 4 downregulated in the male-selected worms compared with the female-selected worms).</footnote>

<footnote><sup>b</sup>P-values were obtained using binom.test function in R (R Core Team 2020), where the minimal sample size is 6; thus, comparisons below 6 transcripts are not applicable for statistical testing and are marked as N/A.</footnote>
provided with a mating partner, which the control and male-
selected lines did (Li et al., unpublished data). In *M. lignano*,
testes also mature prior to ovaries, which is consistent with
high fitness returns via the male sex function in small
individuals (Vizoso and Schärer 2007). It would therefore be
interesting to study changes in gene expression during devel-
opment and in response to altered resource availability and/or
social environment in our experimental evolution lines in f/u-
ture. We also investigated changes in prostate-specific,
antrum-specific, stylet-specific, and cement gland-specific
transcripts identified by Lengerer et al. (2018) and prostate-
specific transcripts identified by Weber et al. (2018), but
found no evidence of any differences. This may simply be
due to a lack of power since there were many more testes-
biased transcripts compared with the other organs (3,360
vs. <400) or because we used a very conservative approach
for comparisons between transcriptome datasets. However,
there are also plausible biological explanations for these
null results. The lack of an effect in prostate-associated and
stylet-associated transcripts (which are related to male func-
tion) might be explained by relatively weak selection in the
male-selected lines (see next section) or that the stylet is a
fixed morphological structure, and thus, differential gene ex-
pression is expected mainly during the development of the
stylet. It is less clear why there was no evidence of a change
in antrum-associated and cement gland-associated tran-
scripts (which are expected to be related to female function),
since the female-selected lines showed the most evidence of
fitness changes (Nordén et al. 2023) and also had the highest
number of differentially expressed transcripts in the pairwise
comparisons (Fig. 1). However, perhaps the antrum is not an
important organ in which to look for changes in gene expres-
sion in our setup, because its main function is to receive sperm
during mating, and so this organ might rather be important to
the outcome of sperm competition (Vizoso et al. 2010). In our
experimental setup, opportunities for sperm competition
and/or cryptic female choice are likely to have been rather
low since individuals were limited to 2 partners (see below).

Potential Confounding factors

There are 3 main potential confounding factors in this experi-
ment, 2 of which are related to the selection protocol. First,
the experimental crossing protocol has likely resulted in lower
levels of sperm competition compared with the ancestral
population. It has been shown that *M. lignano* can mate
with up to 10 different individuals within a day when kept
in large groups (Janicke and Schärer 2009). However, in
our experimental protocol, selected individuals are only pro-
vided with 2 mating partners. This setup was chosen to min-
imize handling time per individual well and allows larger total
population sizes of selected worms, while still maintaining
minimal levels of sperm competition (Nordén et al. 2023).
Sex allocation is therefore expected to be more female-biased

in our setup compared with the ancestral population, since
when sperm competition is low, it is sufficient to produce
relatively few sperm to still be successful (Brauer et al.
2007). This may explain why the male-selected lines seem
to show a weak response to selection compared with the
female-selected lines (the MvsF comparison had the lowest
amount of uniquely differentially expressed transcripts as
well as the lower expression of testes-biased transcripts com-
pared with the controls; see Fig. 1 and Table 1).

Second, use of the GFP marker as a proxy for a sex-
determining locus may inadvertently select for higher GFP ex-
pression, potentially affecting the expression of loci linked to
the GFP locus that do not have any effect on sex function–
specific fitness. We therefore compared GFP expression be-
tween our selection lines and the ancestral BAS1 stock and
found no significant differences (supplementary fig. S4,
Supplementary Material online). Variation in GFP expression
may also come from variation in copy number, since it is
known that inbred laboratory populations can have a high
frequency of abnormal karyotypes and that the most com-
mon variant is a duplication of the large chromosome, where
the GFP locus is believed to be located (Marie-Orleach et al.
2014). We therefore also investigated variation in karyotype
in the selection lines in order to look for any evidence of a
possible bias between the selection regimes. However, we
did not find statistically significant differences in the occur-
rence of abnormal karyotypes between the selection lines
(supplementary table S2, Supplementary Material online). It
is therefore unlikely that the expression differences we have
detected are a spurious consequence of selection on the
GFP expression.

Finally, plasticity in sex allocation is likely to be an important
factor when working with *M. lignano*. Optimal conditions in-
crease sex bias in gonochorists (Wyman et al. 2010), but the
opposite is expected to be true in hermaphrodites, where
low resource availability can reveal tradeoffs between the
sex functions (Schärer et al. 2005). In our experimental proto-
col, individuals were fed ad libitum and kept in a constant
group size as adults, potentially minimizing the effect of any
genes associated with plasticity in sex allocation, in contrast
to a previous transcriptome study, which attempted to identify
sex function-biased genes based on their response to mating
group size (Ramm et al. 2019). However, it seems that plasticity
in expression has not facilitated the adaptive response ob-
served here, because our transcript fold changes did not correl-
ate with the fold changes of corresponding transcripts from
Ramm et al. (2019; supplementary fig. S3, Supplementary
Material online).

Comparisons With Similar Studies

Other experimental evolution studies have manipulated the
intensity of sexual selection in the laboratory and examined
the response at the genomic level, particularly in *Drosophila*.
In general, these experiments support the idea that monogamy (i.e. reduced intensity of sexual selection) results in the evolution of feminized gene expression (Innocenti et al. 2014). Moreover, it has been shown in several species that monogamy decreases the expression of seminal fluid genes involved in postcopulatory sperm competition, which is very important in *M. lignano* (Ramm 2020). This is consistent with our speculation that the experimental protocol used here seems to be more powerful to elicit a response on the female side. To our knowledge, no similar intraspecific genomic studies have been carried out in hermaphrodites, although there are a few that have manipulated the intensity of sexual selection and examined the phenotypic response. In a study of sex-specific sexual selection in the hermaphroditic snail *P. acuta*, there was evidence of an overall decrease in juvenile survival in the monogamous female-selected lines but not in the polygamous male-selected lines, suggesting that strong sexual selection in this species is beneficial through purging of deleterious mutations in the male sex function (Bonel et al. 2018). In another study in *M. lignano*, it was demonstrated that enforced monogamy does not lead to heritable changes in sex allocation (Janicke et al. 2016). One possible explanation for this is that sex allocation is highly plastic in response to social environment in *M. lignano*, which suggests that the heritability of sex allocation may be low. Our results provide an intriguing contrast to both these previous studies. Unlike the *P. acuta* results (Bonel et al. 2018), we did not find any evidence of reduced offspring survival in the female-selected lines (Nordén et al. 2023), while still finding evidence of a putatively adaptive response to selection via downregulation of testes-biased genes (albeit marginally when comparing F and M lines, but more robustly when comparing F and C lines) more often than expected by chance in the female-selected lines (Table 1). Moreover, a recent comparative study in the plant genus *Silene*, which has undergone transition to separate sexes, showed that indeed positive selection strongly drives the downregulation of male-biased genes in females (Prentout et al. 2023). This suggests that the response to selection in our experiment is not simply a result of nonadaptive accumulation of deleterious mutations.

## Implications for the Evolution of Separate Sexes

Despite some similarities to studies that have manipulated the intensity of sexual selection, our main aim with this experiment was to try to gain insights into the early evolution of sex chromosomes and separate sexes. In nature, the most common pathway toward the evolution of separate sexes is the initial accumulation of male deleterious mutations leading to gynodioecy. Indeed, there are a few studies of changes in gene expression in recent transitions to separate sexes in plants (Zemp et al. 2018; Massonnet et al. 2020; Prentout et al. 2023). Zemp et al. (2018) showed that most changes occurred in female *Silene latifolia*, with the

largest changes in male-biased genes, which indicates loss of male function in evolving females. Some male-biased genes also decreased in expression in males compared with sister species (*S. vulgaris*) hermaphrodites, results very similar to ours. Interestingly, although our selection protocol should be most efficient on loci that are physically linked to the GFP locus, we found evidence of a genome-wide response to selection. This is somewhat in contrast to classical predictions about the evolution of sex chromosomes, where loci related to sexual specialization are expected to build up around the sex-determining locus (Massonnet et al. 2020). However, it is consistent with other experiments that have selected on a limited part of the genome (Abbott et al. 2020; Manat 2021) and found evidence of genome-wide changes in expression. This may reflect differences in time scale, where the response to experimental evolution is likely to be determined by the available standing genetic variation (Kawecki et al. 2012), while longer term sexual specialization during sex chromosome evolution may be driven by de novo mutations to a greater extent (Massonnet et al. 2020).

Our results are also comparable with a study of male-limited X chromosome evolution in *Drosophila* (Abbott et al. 2020), which found that male-limited selection resulted in genome-wide masculinization of expression (i.e. upregulation of male-biased genes and downregulation of female-biased genes), as well as a clear signature of changes in genes linked to metabolism. These results are qualitatively similar to ours, where several lines of evidence suggest overall feminization of expression in the female-selected lines (i.e. preferential downregulation of testes-associated genes [Table 1] and involvement of metabolism-related genes in the response to selection [supplementary table S5, Supplementary Material online]). Involvement of metabolic genes is perhaps not unexpected, as sex differences in metabolism are often substantial and can be linked to sex chromosomes (Link and Reue 2017; Ellison and Bachtrog 2019; Abbott et al. 2020; Arnqvist et al. 2022). However, an intriguing aspect here is the fact that metabolism may be linked to sex-specific fitness in different ways in hermaphrodites than in gonochorists. First, while metabolic differences in gonochorists are probably mainly related to body size (Buchholz et al. 2001) and/or various aspects of sexual selection, in hermaphrodites, sex function-specific selection on metabolism may act more through the production of gametes (Schärer et al. 2015). Second, it is known that sex allocation is strongly influenced by resource availability in this and other hermaphroditic species (Vizoso and Schärer 2007), as well as social environment (Marie-Orleach et al. 2017; Ramm et al. 2019). Resource availability and mating group size may therefore have important interactive effects on metabolic activity. Specifically, previous results in this species suggest increased investment into the female function (egg production) and decreased investment into the male function and somatic maintenance when the pool of potential mating partners is
small (Janicke et al. 2013), all of which is likely to impact selection on metabolic genes. Our results therefore add to the growing pool of evidence that metabolism is a key trait in the evolution of sex differences (Arnqvist et al. 2022), including during transitions from hermaphroditism to separate sexes.

## Conclusion

In summary, we found evidence that sex-limited selection in a hermaphrodite results in changes in gene expression that seems likely to be adaptive. Testes-biased genes generally changed in the expected direction (i.e. upregulated in the male-selected lines and downregulated in the female-selected lines), and annotations suggest that metabolic changes may have been an important adaptive component. These results are all consistent with previous data from other species. There was no significant change in expression in transcripts associated with the cement gland, antrum, stylet, or prostate genes, which may partly be due to changes in the intensity of sexual selection because of the experimental protocol. There was no evidence that using the GFP marker as a proxy for a sex-determining locus resulted in any off-target effects and no clear evidence that plasticity has played an important role in the response to selection, suggesting that the main limitation of the experimental protocol was reduced opportunities for sperm competition. In addition, we have identified hundreds of candidate sex function-specific transcripts that would be valuable to study in more detail in the future.

Our results may also have broader implications for transitions from hermaphrodites to separate sexes in nature. Comparative studies have shown that the most common pathway toward the evolution of separate sexes is the initial accumulation of male-deleterious mutations leading to gynodioecy (females and hermaphrodites; Weeks et al. 2006; Beukeboom and Perrin 2014; Zemp et al. 2018; Massonnet et al. 2020; Prentout et al. 2023). Here, we show that sexual specialization in expression is likely to be an important first step after establishment of a sex-determining locus, consistent with theory and previous work in plants (Massonnet et al. 2020).

## Materials and Methods

### Study Species and Experimental Evolution Protocol

*M. lignano* (Macrostomorpha, Rhabditophora, Platyhelminthes, and Lophotrochozoa) is a small (ca. 1.5 mm in length and 0.3 mm in width, containing approximately 25,000 cells) transparent, free-living, regenerating, obligate outcrossing, and simultaneously hermaphroditic flatworm occurring in the intertidal fine sand of the Mediterranean Sea (Ladurner et al. 2005, 2000; Zadesenets et al. 2016). See the **Supplementary Material** online for more details about the structure of the sex organs in this species.

The diploid genome of *M. lignano* consists of 8 chromosomes, 1 pair of which is noticeably larger than the other 3 (Egger and Ishida 2005). However, different karyotypes (e.g. 2n = 9 or 10) have also been found in inbred laboratory populations (Zadesenets et al. 2016, 2017, 2020). The large chromosome itself is a fused duplicate of the 3 small chromosomes, meaning that the regular karyotype is a cryptic tetraploid and additional copies of this chromosome represent pentaploid (2n = 9) and hexaploid (2n = 10) individuals (Zadesenets et al. 2017). These worms do not show noticeable abnormalities in the phenotype and are able to sire viable offspring (Zadesenets et al. 2017). It has a low level of DNA methylation (Wasik et al. 2015), but many transcripts are trans-spliced (Wudarski et al. 2017), increasing transcriptome complexity. The predicted number of genes is 20,000 (Wasik et al. 2015).

In the laboratory, worms are typically kept in populations of 100 individuals in glass Petri dishes filled with 32‰ f/2 medium (Andersen et al. 2005) and fed ad libitum with the diatom alga *Nitzschia curvilineata* in a constant environment of 20 °C and 60% humidity on a 14:10 h light:dark cycle (Ladurner et al. 2005). The generation time in the laboratory takes several weeks. These worms copulate frequently—several times in an hour by performing both sex functions at the same time and lay 1 to 2 eggs per day (Schärer et al. 2004).

Selection lines were started in May 2014, and results in this study are from material collected in generations 21 and 22. Details on the setup and maintenance of the selection lines can be found elsewhere (Nordén et al. 2023). Briefly, all lines originate from the same genetic background (crossing of the transgenic GFP line BAS1 to the wild-type line LS2, to produce worms heterozygous for the GFP gene, which were then backcrossed again to the LS2, creating again heterozygous GFP-positive worms used for the selection and GFP-negative worms used as a "source" population, which provided mating partners for GFP-positive individuals; **supplementary fig. S7, Supplementary Material online**). Each of the selection regimes (male-limited, female-limited, and control) was replicated 4 times, resulting in 12 populations in total (M1-4, F1-4, and C1-4, respectively). Replicate numbers indicate day and time of culturing, not degree of genetic relatedness, which is similar across all populations.

Each selected generation comes from 48 GFP (+) worms crossed to twice as many GFP (−) mates from the corresponding source population. One generation of the experimental evolution protocol takes 5 wk. In the first week, each GFP (+) worm is crossed with 2 GFP (−) worms from their respective "source" population in individual well plates (24-well tissue culture plates; for more detail, see Nordén et al. [2023]). In the second week, worms are moved to new plates for egg laying. Here, the focal worms are allowed to lay eggs for the female-selected lines (selecting for fitness of the GFP [+] worms through the female function, i.e. eggs), while the GFP (−) mating partners are allowed to lay eggs for the male-
selected lines (selecting for fitness of the GFP [+] worms through the male function, i.e. sperm). In the control regime, half of the worms are treated in the same way as the female-selected lines and half in the same way as the male-selected lines. During the third week, parental worms are removed, and offspring left to hatch and grow. In the fourth week, juveniles (≤100) are moved to feeding dishes (glass Petri dishes with algae and f/2 medium) until it is time to start the next generation. This setup means that selection should be most efficient on loci that are located close to the GFP locus.

## Karyotyping

To see if our experimental populations harbor abnormal karyotypes and if whether they were evenly distributed among the selection regimes, we determined the karyotypes of 8 worms per population (32 per selection regime). The single-worm karyotyping was carried out by an air-drying technique as previously described (Zadesenets et al. 2016). For each worm specimen, we checked at least 10 complete metaphases to detect its chromosome number. For all specimens, we found no cases of mosaicism, i.e. all analyzed metaphases in one worm consisted of the same chromosome number. Microimages of chromosome spreads were captured and analyzed using a CCD camera installed on an Axioplan 2 imaging microscope equipped with filtercube #49 (ZEISS, Germany) and using the ISIS4 software package (MetaSystems GmbH, Germany) at the Inter-institutional Shared Center for Microscopic Analysis of Biological Objects (Institute of Cytology and Genetics SB RAS, Novosibirsk, Russian Federation). We tested for a nonrandom distribution of abnormal karyotypes using a binomial generalized linear model, with counts of normal and abnormal karyotypes as the dependent variable, and replicate population and selection regime as the independent variables. Note that since we only have one set of counts per combination of population and selection regime, it is not possible to test for an interaction between them. This analysis was included to check that our selection procedure has not inadvertently resulted in high frequencies of abnormal karyotypes in any selection regime or replicate population. Results are presented in the Supplementary Material online (Karyotypes; supplementary fig. S3, Supplementary Material online).

## RNA Extraction and Sequencing

In both generations 21 and 22, we collected and pooled 40 individuals per replicate from among the worms that had not been randomly selected for producing the next generation (except for replicate F3 in generation 21, when only 35 individuals were available), resulting in 24 samples in total. Additionally, we collected 3 samples, each of 40 individuals, from the BAS1 population to serve as a reference transcriptome for the "ancestral" gene expression and thus enable us to test if there was any change in GFP expression in the

female-selected, male-selected, and control lines compared with the ancestor. All samples were homogenized with a pipette in buffer RLT, before being stored at -80 °C for approximately a month. RNA was extracted from each sample using the RNA purification protocol for animal tissues with RNeasy Mini Kit spin column technology (Qiagen, Hilden, Germany) and DNase I treatment prior to elution in nuclease-free water. RNA quality and quantity was assessed with an Agilent 2100 BioAnalyzer (Agilent Technologies, Santa Clara, CA, USA). Interestingly, *M. lignano* 28S rRNA seems to contain an endogenous "hidden break," as observed in most insects (Winnebeck et al. 2010), thus upon standard denaturation, the hydrogen bonds are disrupted. This results in the release of 2 similar-sized fragments that both migrate together with 18S rRNA, creating one rRNA peak, which makes it impossible to acquire a RNA integrity number value. The resulting rRNA profile should therefore not be interpreted as degraded. As an alternative, the denaturation step can be excluded, resulting in dual peaks (supplementary fig. S1, Supplementary Material online). At this stage, 3 samples were excluded due to low RNA concentration (F1 in generation 21, and C1 and M3 in generation 22). After quality control, the remaining samples were sent for cluster generation and 125 cycles paired-end sequencing with the HiSeq2500 system and v4 sequencing chemistry (Illumina Inc., San Diego, CA, USA) at Science for Life Laboratory (SciLifeLab, SNP&SEQ Technology Platform, Uppsala, Sweden).

## Bioinformatics and Statistics

The quality and properties of raw sequencing data were examined using the nf-core/rnaseq pipeline (rnaseq-1.3, Nextflow version 19.07.0; Ewels et al. 2019). The genome features table from the Wudarski et al. (2017) version of the whole genome assembly Mlig_3_7 was converted to GTF format using Cufflinks version 2.2.1 (Trapnell et al. 2010) and used to run the pipeline. This initial analysis revealed that, though the data quality was in general good, a large proportion of the sequence reads mapped to multiple locations (only about 40% of the reads mapped uniquely to one locus). This result is most likely explained by redundancy within the genome assembly due to a recent whole genome duplication event in this species (forming the large chromosome), as well as differences in the chromosome numbers between lines (the DV1 line used for the genome assembly Mlig_3_7 has an extra copy of the large chromosome). To overcome this limitation and increase the number of reads available for downstream analysis, a transcriptome-based approach was used to generate a line-specific reference. First, reads were trimmed, and remaining adaptors were removed using Trimmomatic (version 0.36; Bolger et al. 2014). A de novo transcriptome assembly using the 3 ancestral GFP line samples was then constructed using Trinity version 2.4.0 (Grabherr et al. 2011). This resulted in a 531-Mb
transcriptome assembly containing 889,873 transcripts, which is considerably higher than the number of predicted genes for this species. Isoforms, noncoding RNA, lowly expressed genes, heterozygosity within the sequenced population, and sequencing errors likely account for the large number of transcripts in the transcriptome assembly (Babarinde and Hutchins 2022). For the purpose of this study however, the number of transcripts should not play a major role since lowly expressed transcripts and transcripts with few/no reads that map will not be carried forward in the differential expression analysis. Bowtie2 (version 2.3.4.3) was used to map the raw reads back to the assembly, yielding a mapping rate of ~83%. A BUSCO (version 2.0.1; Seppey et al. 2019) search of the metazoan_odb9 dataset returned 85.1% complete genes, 16.5% single genes, and 68.8% duplicated genes. Mapping the data to the *M. lignano* genome assembly rather than the de novo transcriptome resulted in similar numbers of DE transcripts (data not shown), which suggests that the choice of reference does not have a large impact on the downstream results.

To annotate the transcriptome, we performed a BLASTX (version 2.6.0+; Altschul et al. 1990) search against Swiss-Prot (accessed on May 27, 2019) using an e-value cutoff of 1e-10. Only hits that were at least 25% identical over a minimum alignment length of 100 were retained. InterProScan (version 5.30-69.0; Jones et al. 2014) was used to classify transcripts into protein families. This search included the Pfam, ProDom, and SuperFamily member databases. GO terms and InterPro entry information were included in these results. To obtain Entrez gene ids for downstream analyses in EdgeR, a BLASTX search against the mouse Ensembl proteins (release 38) was done using an e-value cutoff of 1e-10. The corresponding Entrez gene ids were retrieved using the mygene python package (version 3.1.0; https://pypi.org/project/mygene/).

To quantify transcript expression, Salmon (version 0.9.1) was used in quasi-mapping mode with 100 bootstrap replicates (Patro et al. 2017). A kmer length of 31 was used for indexing. Sample C3 in generation 22 was excluded from the differential expression analysis because of quality issues. Differentially expressed transcripts between selection regimes (male-limited, female-limited, and control) were identified with EdgeR (version 3.24.3; Robinson et al. 2010). Transcripts that had fewer than 0.5 counts-per-million in 3 or more samples were excluded from the analysis. Transcripts with an absolute log2-transformed fold change >1 and adjusted P-value < 0.05 (after Benjamini-Hochberg false discovery rate [FDR] correction) were considered significantly differentially expressed between the selection regimes. We also analyzed expression data using relative allocation to testes and ovaries as a covariate to control for any selection regime differences in total amount of testis/ ovary tissue, but this produced qualitatively similar results in the downstream analyses, so these results are not

reported. This result is not unexpected, since the lines do not differ significantly in mean testis and ovary size measured under standardized conditions (Nordén et al. 2023; Li et al., unpublished data). Volcano plots using EdgeR output tables were also created with ggplot2 (Wickham 2016) and library (ggrepel). Due to its role as an analog of a sex-determining locus, one of the predictions that can be made from models of sex chromosome evolution is that differentiation among selection lines should be found preferentially in loci located near the GFP locus. We therefore tested for differentiation among selection lines for all transcripts located on the GFP scaffold using paired Wilcoxon signed-rank tests. We also tested whether the mean change in expression in M and F lines compared with C lines was higher on the GFP scaffold than for the rest of the genome using unpaired rank sum Wilcoxon tests (including all transcripts and not only the significant ones). Finally, we checked whether the GFP scaffold was an outlier compared with the rest of the genome by examining its place in the distribution of average change in expression per scaffold (i.e. whether it fell outside the 2.5 and 97.5 percentiles of the distribution). Note that the GFP scaffold comprises only a small part of the large chromosome, but since we do not know where other scaffolds are located, we cannot check for linkage effects on the scaffolds next to the GFP scaffold.

GO and KEGG pathway analyses were conducted using the goana function in the EdgeR package. The results of the GO term enrichment and KEGG pathway analyses were adjusted for multiple testing using FDR correction, and terms with an adjusted P-value < 0.05 were considered statistically significant. Additionally, we looked at the standardized fold change in previously identified sex organ-biased transcripts (cement glands, antrum, prostate, and stylet [Lengerer et al. 2018], as well as for prostate in Weber et al. [2018] and ovary, testes, gonad [most likely ovary biased], and tail [Ramm et al. 2019]) in each selection regime as blastn queries (version 2.6.0+) to search in our transcriptome assembly. Note that all of the transcripts in Ramm et al. (2019) were originally identified as ovary, testes, and tail specific in Arbore et al. (2015). We elected to compare our data primarily with results from Ramm et al. (2019) since this study examined how expression changed according to social group size, and we felt that such socially sensitive genes could be considered more robust organ-specific candidates and likely to be implicated in sex-specific functions. Sequences that were at least 95% identical over a minimum of 500 bases were considered matches. Transcripts with organ biased expression profiles linked to the ovaries, cement glands, and antrum were considered female-function biased, and transcripts linked to the testes, prostate, and stylet were considered male-function biased. We also analyzed the transcripts expressed in both gonads (however they are more likely ovary biased) and tail gene groups, since many genes in these regions are known to have sex-related functions (Arbore et al. 2015; Ramm et al. 2019). However, as we used stringent
criteria and did not include transcripts that were shorter than 500 bases, for example, and the mapping was not always one to one, that complicates comparisons between the different studies. Exact binomial tests were performed to test for over-representation of specific categories of transcripts (e.g. testis biased, ovary biased) among the transcripts that had significantly different expression between the selection regimes or our and Ramm et al. (2019) studies. We checked whether there was in situ hybridization data available for any of these transcripts but only found one gonad-specific transcript TRINITY_DN115795_c0_g1_i10 (RNA815_1618.1), which is downregulated in FvsC and upregulated in MvsF and has in situ hybridization data from Arbore et al. (2015), showing expression in ovaries. To test the correlation of fold changes between our dataset and Ramm et al. (2019) transcripts from Table 1, we used Pearson correlation (library(stringi)). Moreover, PCA analyses were carried out on all sex-biased transcripts identified by Ramm et al. (2019) using the prcomp function and then visualized with ggplot2 (Wickham 2016) in R to test for an axis of overall feminization/masculinization, which could separate the selection regimes to some extent, and see which transcripts might be driving that separation. PC statistical significance was tested by creating a linear mixed model (library(lmerTest)) dependent on selection regime and controlled for population with significance testing using ana-va() and Tukey test (library(multcomp)). Graphs and statistics were obtained using R version 3.6.3 (R Core Team 2020).

## Supplementary Material

Supplementary material is available at *Genome Biology and Evolution* online.

## Acknowledgments

We thank Lukas Schärer for contributing with important feedback when planning the study, Tomas Johansson for the RNA sample assessment on the Agilent 2100 BioAnalyzer, and Vignesh Ramnath for optimizing specific transcript TPM script in python. We would also like to thank National Genomics Infrastructure (NGI) at SciLifeLab for providing sequencing service. Bioinformatics support was received through National Bioinformatics Infrastructure Sweden (NBIS) long-term support, which provided support for A.C. through the Knut and Alice Wallenberg Foundation. A.C. was additionally supported by Jörgen Lindströms Stipendiefond. The computations were enabled by resources provided by the National Academic Infrastructure for Supercomputing in Sweden (NAISS) and the Swedish National Infrastructure for Computing (SNIC) at UPPMAX partially funded by the Swedish Research Council through grant agreements no. 2022-06725 and no. 2018-05973. K.S.Z. was supported by project FWNR-2022-0015. This study was funded by grants to J.K.A from Maja och Erik Lindqvists forskningsstiftelse,

Nilsson-Ehle-donationerna, the Crafoord Foundation (20120628 and 20140644), VR (2011-05679 and 2015-04680), and ERC-StG-2015-678148.

## Data Availability

The raw sequence reads have been deposited in NCBI's Gene Expression Omnibus (GEO) and are accessible through GEO Series accession number GSE242254 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE242254).

## Literature Cited

* Abbott JK, Chippindale AK, Morrow EH. The microevolutionary response to male-limited X-chromosome evolution in *Drosophila melanogaster* reflects macroevolutionary patterns. *J Evol Biol.* 2020;33(6):738–750. https://doi.org/10.1111/jeb.13618.

* Altschul SF, Gish W, Miller W, Myers EW, Lipman DJ. Basic local alignment search tool. *J Mol Biol.* 1990;215(3):403–410. https://doi.org/10.1016/S0022-2836(05)80360-2.

* Andersen R, Berges J, Harrison P, Watanabe MM. Recipes for freshwater and seawater media. In: Andersen RA, editor. Algal culturing techniques. Amsterdam: Elsevier; 2005. p. 429–538.

* Arbore R, Sekii K, Beisel C, Ladurner P, Berezikov E, Schärer L. Positional RNA-Seq identifies candidate genes for phenotypic engineering of sexual traits. *Front Zool.* 2015;12(1):14. https://doi.org/10.1186/s12983-015-0106-0.

* Arnvqvist G, Ronn J, Watson C, Goenaga J, Immonen E. Concerted evolution of metabolic rate, economics of mating, ecology, and pace of life across seed beetles. *Proc Natl Acad Sci U S A.* 2022;119(33): e2205564119. https://doi.org/10.1073/pnas.2205564119.

* Babarinde IA, Hutchins AP. The effects of sequencing depth on the assembly of coding and noncoding transcripts in the human genome. *BMC Genomics.* 2022;23(1):487. https://doi.org/10.1186/s12864-022-08717-z.

* Bachtrog D, Mank JE, Peichel CL, Kirkpatrick M, Otto SP, Ashman TL, Hahn MW, Kitano J, Mayrose I, Ming R, et al. Sex determination: why so many ways of doing it? *PLoS Biol.* 2014;12(7):e1001899. https://doi.org/10.1371/journal.pbio.1001899.

* Beukeboom LW, Perrin N. 2014. The evolution of sex determination. Oxford: Oxford University Press.

* Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data. *Bioinformatics.* 2014;30(15):2114–2120. https://doi.org/10.1093/bioinformatics/btu170.

* Bonel N, Noël E, Janicke T, Sartori K, Chapuis E, Ségard A, Meconcelli S, Pélissié B, Sarda V, David P. Asymmetric evolutionary responses to sex-specific selection in a hermaphrodite. *Evolution.* 2018;72(10): 2181–2201. https://doi.org/10.1111/evo.13565.

* Brauer VS, Schärer L, Michiels NK. Phenotypically flexible sex allocation in a simultaneous hermaphrodite. *Evolution.* 2007;61:216–222. https://doi.org/10.1111/j.1558-5646.2007.00018.x.

* Buchholz AC, Rafii M, Pencharz PB. Is resting metabolic rate different between men and women? *Br J Nutr.* 2001;86(6):641–646. https://doi.org/10.1079/bjn2001471.

* Charlesworth D. Young sex chromosomes in plants and animals. *New Phytol.* 2019;224(3):1095–1107. https://doi.org/10.1111/nph.16002.

* Edvardsen RB, Wallerman O, Furmanek T, Kleppe L, Jern P, Wallberg A, Kjærner-Semb E, Mæhle S, Olausson SK, Sundström E, et al. Heterochiasmy and the establishment of gsdf as a novel sex determining gene in Atlantic halibut. *PLoS Genet.* 2022;18(2): e1010011. https://doi.org/10.1371/journal.pgen.1010011.
Egger B, Ishida S. Chromosome fission or duplication in *Macrostomum lignano* (Macrostomorpha, Plathelminthes)—remarks on chromosome numbers in ‘archoophoran turbellarians’. J Zool Syst Evol Res. 2005;43(2):127–132. https://doi.org/10.1111/j.1439-0469.2005.00300.x

Ellison C, Bachtrog D. Recurrent gene co-amplification on Drosophila X and Y chromosomes. PLoS Genet. 2019;15(7):e1008251. https://doi.org/10.1371/journal.pgen.1008251.

Ewels PA, Peltzer A, Fillinger S, Alneberg J, Patel H, Wilm A, Garcia MU, Di Tommaso P, Nahnsen S. nf-core: Community curated bioinformatics pipelines. bioRxiv 610741. https://doi.org/10.1101/610741, 16 April 2019, preprint: not peer reviewed.

Grabherr MG, Haas BJ, Yassour M, Levin JZ, Thompson DA, Amit I, Adiconis X, Fan L, Raychowdhury R, Zeng Q, et al. Full-length transcriptome assembly from RNA-Seq data without a reference genome. Nat Biotechnol. 2011;29(7):644–652. https://doi.org/10.1038/nbt.1883.

Grath S, Parsch J. Sex-biased gene expression. Annu Rev Genet. 2016;50(1):29–44. https://doi.org/10.1146/annurev-genet-120215-035429.

Hollis B, Houle D, Yan Z, Kawecki TJ, Keller L. Evolution under monogamy feminizes gene expression in *Drosophila melanogaster*. Nat Commun. 2014;5(1):3482. https://doi.org/10.1038/ncomms4482.

Innocenti P, Flis I, Morrow EH. Female responses to experimental removal of sexual selection components in *Drosophila melanogaster*. BMC Evol Biol. 2014;14(1):239. https://doi.org/10.1186/s12862-014-0239-3.

Janicke T, Marie-Orleach L, De Mulder K, Berezikov E, Ladurner P, Vizoso DB, Schärer L. Sex allocation adjustment to mating group size in a simultaneous hermaphrodite. Evolution. 2013;67(11): 3233–3242. https://doi.org/10.1111/evo.12189.

Janicke T, Sandner P, Ramm SA, Vizoso DB, Schärer L. Experimentally evolved and phenotypically plastic responses to enforced monogamy in a hermaphroditic flatworm. J Evol Biol. 2016;29(9): 1713–1727. https://doi.org/10.1111/jeb.12910.

Janicke T, Schärer L. Determinants of mating and sperm-transfer success in a simultaneous hermaphrodite. J Evol Biol. 2009;22(2): 405–415. https://doi.org/10.1111/j.1420-9101.2008.01660.x.

Joag R, Stuglik M, Konczal M, Plesnar-Bielak A, Skrzynecka A, Babik W, Radwan J. Transcriptomics of intralocus sexual conflict: gene expression patterns in females change in response to selection on a male secondary sexual trait in the bulb mite. Genome Biol Evol. 2016;8(8):2351–2357. https://doi.org/10.1093/gbe/evw169.

Jones P, Binns D, Chang HY, Fraser M, Li W, McAnulla C, McWilliam H, Maslen J, Mitchell A, Nuka G, et al. InterProScan 5: genome-scale protein function classification. Bioinformatics. 2014;30(9): 1236–1240. https://doi.org/10.1093/bioinformatics/btu031.

Kawecki TJ, Lenski RE, Ebert D, Hollis B, Olivieri I, Whitlock MC. Experimental evolution. Trends Ecol Evol. 2012;27(10):547–560. https://doi.org/10.1016/j.tree.2012.06.001.

Ladurner P, Rieger R, Baguñà J. Spatial distribution and differentiation potential of stem cells in hatchlings and adults in the marine platyhelminth *Macrostomum* sp.: a bromodeoxyuridine analysis. Dev Biol. 2000;226(2):231–241. https://doi.org/10.1006/dbio.2000.9867.

Ladurner P, Schärer L, Salvenmoser W, Rieger RM. A new model organism among the lower bilateria and the use of digital microscopy in taxonomy of meiobenthic Platyhelminthes: *Macrostomum lignano*, n. sp. (Rhabditophora, Macrostomorpha). J Zool Syst Evol Res. 2005;43(2): 114–126. https://doi.org/10.1111/j.1439-0469.2005.00299.114.

Lengerer B, Wunderer J, Pjeta R, Carta G, Kao D, Aboobaker A, Beisel C, Berezikov E, Salvenmoser W, Ladurner P. Organ specific gene expression in the regenerating tail of *Macrostomum lignano*. Dev Biol. 2018;433(2):448–460. https://doi.org/10.1016/j.ydbio.2017.07.021.

Link JC, Reue K. Genetic basis for sex differences in obesity and lipid metabolism. Annu Rev Nutr. 2017;37(1):225–245. https://doi.org/10.1146/annurev-nutr-071816-064827.

Manat Y. The genomic basis of the response to female-limited X-chromosome evolution. Lund: Lund University; 2021.

Marie-Orleach L, Janicke T, Vizoso DB, Eichmann M, Schärer L. Fluorescent sperm in a transparent worm: validation of a GFP marker to study sexual selection. BMC Evol Biol. 2014;14(1):148. https://doi.org/10.1186/1471-2148-14-148.

Marie-Orleach L, Vogt-Burri N, Mouginot P, Schlatter A, Vizoso DB, Bailey NW, Schärer L. Indirect genetic effects and sexual conflicts: partner genotype influences multiple morphological and behavioral reproductive traits in a flatworm. Evolution. 2017;71(5): 1232–1245. https://doi.org/10.1111/evo.13218.

Massonnet M, Cochetel N, Minio A, Vondras AM, Lin J, Muyle A, Garcia JF, Zhou Y, Delledonne M, Riaz S, et al. The genetic basis of sex determination in grapes. Nat Commun. 2020;11(1):2902. https://doi.org/10.1038/s41467-020-16700-z.

Nordén AK, Ramm SA, Abbott JK. Rapid evolution of sex role specialization in a hermaphrodite under sex-limited selection. Evolution. 2023;77(4): 1066–1076. https://doi.org/10.1101/2022.04.21.489077.

Parsch J, Ellegren H. The evolutionary causes and consequences of sex-biased gene expression. Nat Rev Genet. 2013;14(2):83–87. https://doi.org/10.1038/nrg3376.

Patro R, Duggal G, Love MI, Irizarry RA, Kingsford C. Salmon provides fast and bias-aware quantification of transcript expression. Nat Methods. 2017;14(4):417–419. https://doi.org/10.1038/nmeth.4197.

Prentout D, Muyle A, Zemp N, El Filali A, Boussau B, Touzet P, Widmer A, Kafer J, Marais GA. Evolution of sex-biased gene expression during transitions to separate sexes in the Silene genus. bioRxiv 560480. https://doi.org/10.1101/2023.10.02.560480, 03 October 2023, preprint: not peer reviewed.

Rafati N, Chen J, Herpin A, Pettersson ME, Han F, Feng C, Wallerman O, Rubin CJ, Péron S, Cocco A, et al. Reconstruction of the birth of a male sex chromosome present in Atlantic herring. Proc Natl Acad Sci U S A. 2020;117(39):24359–24368. https://doi.org/10.1073/pnas.2009925117.

Ramm SA. Seminal fluid and accessory male investment in sperm competition: seminal fluid and sperm competition. Phil Trans R Soc B: Biol Sci. 2020;375(1813):20200068. https://doi.org/10.1098/rstb.2020.0068rstb20200068.

Ramm SA, Lengerer B, Arbore R, Pjeta R, Wunderer J, Giannakara A, Berezikov E, Ladurner P, Schärer L. Sex allocation plasticity on a transcriptome scale: socially sensitive gene expression in a simultaneous hermaphrodite. Mol Ecol. 2019;28(9):2321–2341. https://doi.org/10.1111/mec.15077.

R Core Team. 2020. R: A Language and Environment for Statistical Computing. [accessed 2023 Dec 14]. https://www.r-project.org/.

Robinson MD, McCarthy DJ, Smyth GK. Edger: a bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. 2010;26(1):139–140. https://doi.org/10.1093/bioinformatics/btp616.

Schärer L, Janicke T, Ramm SA. Sexual conflict in hermaphrodites. Cold Spring Harb Perspect Biol. 2015;7(1):a017673. https://doi.org/10.1101/cshperspect.a017673.

Schärer L, Joss G, Sandner P. Mating behaviour of the marine turbellarian *Macrostomum* sp.: these worms suck. Mar Biol. 2004;145(2): 373–380. https://doi.org/10.1007/s00227-004-1314-x.

Schärer L, Ladurner P. 2003. Phenotypically plastic adjustment of sex allocation in a simultaneous hermaphrodite. Proc R Soc Lond B. 270(1518):935–941. https://doi.org/10.1098/rspb.2002.2323

Schärer L, Sandner P, Michiels NK. Trade-off between male and female allocation in the simultaneously hermaphroditic flatworm
Macrostomum sp. J Evol Biol. 2005:18(2):396-404. https://doi.org/10.1111/j.1420-9101.2004.00827.x.

Seppey M, Manni M, Zdobnov EM. BUSCO: assessing genome assembly and annotation completeness. Methods Mol Biol. 2019:1962: 227-245. https://doi.org/10.1007/978-1-4939-9173-0_14.

Trapnell C, Williams BA, Pertea G, Mortazavi A, Kwan G, van Baren MJ, Salzberg SL, Wold BJ, Pachter L. Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation. Nat Biotechnol. 2010:28(5): 511-515. https://doi.org/10.1038/nbt.1621.

Vellnow N, Marie-Orleach L, Zadesenets KS, Schärer L. Bigger testes increase paternity in a simultaneous hermaphrodite, independently of the sperm competition level. J Evol Biol. 2018:31(2):180-196. https://doi.org/10.1111/jeb.13212.

Veltsos P, Fang Y, Cossins AR, Snook RR, Ritchie MG. Mating system manipulation and the evolution of sex-biased gene expression in Drosophila. Nat Commun. 2017:8(1):2072. https://doi.org/10.1038/s41467-017-02232-6.

Vizoso DB, Rieger G, Schärer L. Goings-on inside a worm: functional hypotheses derived from sexual conflict thinking. Biol J Linnean Soc. 2010:99(2):370-383. https://doi.org/10.1111/j.1095-8312.2009.01363.x.

Vizoso DB, Schärer L. Resource-dependent sex-allocation in a simultaneous hermaphrodite. J Evol Biol. 2007:20(3):1046-1055. https://doi.org/10.1111/j.1420-9101.2007.01294.x.

Wasik K, Gurtowski J, Zhou X, Ramos OM, Delás MJ, Battistoni G, El Demerdash O, Falcatori I, Vizoso DB, Smith AD, et al. Genome and transcriptome of the regeneration- competent flatworm, *Macrostomum lignano*. Proc Natl Acad Sci U S A. 2015:112(40): 12462-12467. https://doi.org/10.1073/pnas.1516718112.

Weber M, Wunderer J, Lengerer B, Pjeta R, Rodrigues M, Schärer L, Ladurner P, Ramm SA. A targeted in situ hybridization screen identifies putative seminal fluid proteins in a simultaneously hermaphroditic flatworm. BMC Evol Biol. 2018:18(1):81. https://doi.org/10.1186/s12862-018-1187-0.

Weeks SC, Benvenuto C, Reed SK. 2006. When males and hermaphrodites coexist: a review of androdioecy in animals. Integr Comp Biol. 46(4):449-513. https://doi.org/10.1093/icb/icj048.

Wiberg RAW, Veltsos P, Snook RR, Ritchie MG. Experimental evolution supports signatures of sexual selection in genomic divergence. Evol Lett. 2021:5(3):214-229. https://doi.org/10.1002/evl3.220.

Wickham H. Ggplot2: elegant graphics for data analysis. New York (NY): Springer-Verlag; 2016.

Winnebeck EC, Millar CD, Warman GR. Why does insect RNA look degraded? J Insect Sci. 2010:10:159. https://doi.org/10.1673/031.010.14119.

Wudarski J, Simanov D, Ustyantsev K, de Mulder K, Grelling M, Grudniewska M, Beltman F, Glazenburg L, Demircan T, Wunderer J, et al. Efficient transgenesis and annotated genome sequence of the regenerative flatworm model *Macrostomum lignano*. Nat Commun. 2017:8(1):2120. https://doi.org/10.1038/s41467-017-02214-8.

Wyman MJ, Agrawal AF, Rowe L. Condition-dependence of the sexually dimorphic transcriptome in *Drosophila melanogaster*. Evolution. 2010:64(6):1836-1848. https://doi.org/10.1111/j.1558-5646.2009.00938.x.

Zadesenets KS, Jetybayev IY, Schärer L, Rubtsov NB. Genome and karyotype reorganization after whole genome duplication in free-living flatworms of the genus Macrostomum. Int J Mol Sci. 2020:21(2):680. https://doi.org/10.3390/ijms21020680.

Zadesenets KS, Schärer L, Rubtsov NB. New insights into the karyotype evolution of the free-living flatworm *Macrostomum lignano* (Platyhelminthes, Turbellaria). Sci Rep. 2017:7(1):6066. https://doi.org/10.1038/s41598-017-06498-0.

Zadesenets KS, Vizoso DB, Schlatter A, Konopatskaia ID, Berezikov E, Schärer L, Rubtsov NB. Evidence for karyotype polymorphism in the free-living flatworm, *Macrostomum lignano*, a model organism for evolutionary and developmental biology. PLoS One. 2016:11(10):e0164915. https://doi.org/10.1371/journal.pone.0164915.

Zemp N, Widmer A, Charlesworth D. Has adaptation occurred in males and females since separate sexes evolved in the plant *Silene latifolia*? Proc R Soc B: Biol Sci. 2018:285(1883):20172824. https://doi.org/10.1098/rspb.2017.2824.

Associate editor: Andrea Betancourt