# Identification of General Patterns of Sex-Biased Expression in Daphnia, a Genus with Environmental Sex Determination

Cécile Molinier, Céline Reisser, Peter Fields, Adeline Ségard, Yan Galimov,
Christoph R. Haag

► To cite this version:

Cécile Molinier, Céline Reisser, Peter Fields, Adeline Ségard, Yan Galimov, et al.. Identification of General Patterns of Sex-Biased Expression in Daphnia, a Genus with Environmental Sex Determination. G3, 2018, 8 (5), pp.1523-1533. 10.1534/g3.118.200174. hal-04202122

HAL Id: hal-04202122

https://hal.science/hal-04202122v1

Submitted on 12 Dec 2024

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
# Identification of General Patterns of Sex-Biased Expression in *Daphnia*, a Genus with Environmental Sex Determination

Cécile Molinier,*,¹,² Céline M.O. Reisser,*,†,*,¹ Peter Fields,§ Adeline Ségard,* Yan Galimov,** and Christoph R. Haag*,†

*Centre d'Ecologie Fonctionnelle et Evolutive (CEFE)-Unité Mixte de Recherche 5175, Centre National de la Recherche Scientifique (CNRS)-Université de Montpellier-Université Paul-Valéry Montpellier-Ecole Pratique des Hautes Etudes (EPHE), 1919 Route de Mende, 34293 Montpellier Cedex 5, France, †Université de Fribourg, Ecology and Evolution, Ch. du Musée 10, 1700 Fribourg, Switzerland, ‡IFREMER Centre du Pacifique, UMR 241 EIO, Labex CORAIL, BP 49, 98719 Taravao, Tahiti, Polynésie Française, §Universität Basel, Zoology Institute, Evolutionary Biology, Vesalgasse 1, 4051 Basel, Switzerland, and **Koltsov Institute of Developmental Biology RAS ul. Vavilova 26, 119334 Moscow, Russia

ORCID IDs: 0000-0003-2959-2524 (P.F.); 0000-0002-8817-1431 (C.R.H.)

**ABSTRACT** *Daphnia* reproduce by cyclic-parthenogenesis, where phases of asexual reproduction are intermit-
ted by sexual production of diapause stages. This life cycle, together with environmental sex determination,
allow the comparison of gene expression between genetically identical males and females. We investigated
gene expression differences between males and females in four genotypes of *Daphnia magna* and compared
the results with published data on sex-biased gene expression in two other *Daphnia* species, each representing
one of the major phylogenetic clades within the genus. We found that 42% of all annotated genes showed sex-
biased expression in *D. magna*. This proportion is similar both to estimates from other *Daphnia* species as well as
from species with genetic sex determination, suggesting that sex-biased expression is not reduced under
environmental sex determination. Among 7453 single copy, one-to-one orthologs in the three *Daphnia*
species, 707 consistently showed sex-biased expression and 675 were biased in the same direction in all
three species. Hence these genes represent a core-set of genes with consistent sex-differential expres-
sion in the genus. A functional analysis identified that several of them are involved in known sex de-
termination pathways. Moreover, 75% were overexpressed in females rather than males, a pattern that
appears to be a general feature of sex-biased gene expression in *Daphnia*.

KEYWORDS

Differential gene
expression
males
females
*Daphnia*
ESD
Genetics of Sex

Patterns of gene expression often differ strongly between male and
female individuals of the same species (Ellegren & Parsch 2007). In part,
this difference is driven by genes on sex chromosomes, which show a

particularly strong tendency for sex-biased (or sex-differential) expres-
sion (Ellegren & Parsch 2007; Bergero & Charlesworth 2009; Grath &
Parsch 2016). However, sex-biased expression also occurs in many
autosomal genes, and their products also fundamentally contribute to
differences between male and female phenotypes (Ellegren & Parsch
2007; Grath & Parsch 2016; Wright *et al.* 2017). A particularly inter-
esting case is that of species with environmental sex determination
(ESD), where the same genotype may develop into a male or female,
depending on environmental cues. Pure ESD species do not have sex
chromosomes, and sex differentiation entirely relies on autosomal
genes (Bull 1985). In GSD species, on the other hand, sex chromosomes
contain a particularly high number of sex-biased genes (Mank 2009;
Grath & Parsch 2016), and some of the genetic differences between
sexes (e.g., sex-specific genomic regions or allelic differences) cause
further, downstream expression differences, including for autosomal
genes (Yang *et al.* 2006; Wijchers & Festenstein 2011). Both these

<footnote>Copyright © 2018 Molinier et al.
doi: https://doi.org/10.1534/g3.118.200174
Manuscript received November 7, 2017; accepted for publication February 25,
2018; published Early Online March 13, 2018.
This is an open-access article distributed under the terms of the Creative Commons
Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/),
which permits unrestricted use, distribution, and reproduction in any medium,
provided the original work is properly cited.
Supplemental material available at Figshare: https://doi.org/10.25387/
g3.5965285
¹Equal contribution
²Corresponding author: Cécile Molinier, Centre d'Ecologie Fonctionnelle et
Evolutive CEFE UMR 5175, 1919, route de Mende, 34293 Montpellier Cedex 5,
France, e-mail: cecile.molinier@cefe.cnrs.fr</footnote>
observations suggest that species with ESD may have a lower number of genes with sex-biased expression compared to species with GSD. Alternatively, however, species with ESD may show a higher number of genes with sex-specific expression than species with GSD, because no genetic differences exist between sexes, and their entire sex-specific phenotypes are by differential gene expression (Grath & Parsch 2016, Mank 2017).

Among the species with ESD that have so far been investigated for sex-biased expression (Torres Maldonado et al. 2002; Shoemaker et al. 2007; Yatsu et al. 2016; Radhakrishnan et al. 2017), we find two species of the genus *Daphnia* (Colbourne et al. 2011; Huylmans et al. 2016). *Daphnia* reproduce by cyclical parthenogenesis: during the asexual phase of the life cycle, mothers clonally produce sons or daughters, but this asexual phase is intermitted by sexual reproduction, leading to the production of diapause eggs, which give rise to female hatchlings. Males and females are morphologically distinct (Scourfield & Harding 1966), and the sex of the clonally produced offspring is determined by environmental factors such as shortened photoperiod and/or increased population density (Roulin et al. 2013; Korpelainen 1990; Hobaek & Larsson 1990). Specifically, the production of males is induced by a hormone emitted by the mother in response to environmental conditions (Olmstead & Leblanc 2002). Moreover, male production can also be induced experimentally by adding the hormone analog methyl farnesoate (MF) to the culture medium at a precise moment of the ovarian cycle corresponding to 48 to 72h after moulting (Olmstead & Leblanc 2002).

Since the publication of the *D. pulex* genome (Colbourne et al. 2011), the amount of genomic and transcriptomic resources for the genus has markedly increased (Routtu et al. 2014; Xu et al. 2015; Dukić et al. 2016; Orsini et al. 2016; Huylmans et al. 2016; Giraudo et al. 2017; Lynch et al. 2017; Spanier et al. 2017; Toyota et al. 2017; Ye et al. 2017; Herrmann et al. 2018). Previous studies have investigated sex-biased gene expression in two *Daphnia* species (Colbourne et al. 2011; Eads et al. 2007; Huylmans et al. 2016). The two species, *D. pulex* and *D. galeata* each belong to one of the two major phylogenetic groups within the subgenus *Daphnia* (Colbourne & Hebert 1997; Ishida et al. 2006; Adamowicz, Petrusek & Colbourne 2009). These studies reported a high number of genes with sex-biased expression in both species and a preponderance of female-biased genes (i.e., genes overexpressed in females) compared to male-biased genes.

The genus *Daphnia* contains a second subgenus, *Ctenodaphnia*, which notably contains the species *D. magna*. This species is not only one of the major genomic model organisms of the genus (Miner et al. 2012, GenBank accession number: LRGB00000000.1), but also for studies on sex differentiation under ESD and other sex-related traits, including local adaptation in male production, evolution of genetic sex determination (which occurs in some genotypes of *D. magna*, not investigated here), and uniparental reproduction (Kato et al. 2011; Galimov, Walser & Haag 2011; Svendsen et al. 2015; Reisser et al. 2016; Roulin et al. 2013). Here we present an analysis of sex-biased gene expression in *D. magna*, based on RNA-sequencing of males and females of four different genotypes (males and females of the same genotypes are members of the same clone). The four genotypes were used as biological replicates, that is, the expression of individual genes was classified as “sex-biased” only if the bias was consistent among all genotypes tested (i.e., if the overall pairwise test with four replicates was significant). While the primary aim was to study sex-biased gene expression in this important model organism, we also compare our results to those from the previous studies on *D. pulex* and *D. galeata*. The aims of this comparison were to identify general patterns of sex-biased gene expression in the genus, and to identify a core-set of genes with consistent sex-biased expression within the genus.

# MATERIAL AND METHODS

## Study design and origin of clones

We carried out RNA sequencing on adult males and females of *D. magna*, reared under standard culturing conditions (see below). We used four different genotypes which originated from a single population in Moscow, Russia (N55°45'48.65", E37°34'54.00"), as biological replicates. One library was prepared per genotype and sex, resulting in a total of eight libraries. Furthermore, each library consisted of eight technical (experimental) replicates, that is, eight clonal individuals of the same genotype and sex, pooled together into the library. Hence, a total of 64 individuals were raised for the experiment. Technical replicates were used to reduce variation due to small differences in environmental conditions (light, temperature, food, etc.) on gene expression. Such small environmental differences may be caused, for instance, by different positions of individuals within the culture tubes in the culture chamber.

## Laboratory culturing

Gravid parthenogenetic females were transferred individually to standard culturing conditions: a single individual in a 50 mL Falcon tube containing 20 mL of artificial medium for *Daphnia* (Klüttgen et al. 1994), fed with 150 µL of algae solution (50 million of cells of *Scenedesmus* sp. per mL), and kept at 19° under a 16:8 hr light-dark photoperiod. Each technical replicate was reared under these standard conditions during two pre-experimental clonal generations to remove potential material effects (Gorbi et al. 2011). To that end, one randomly selected offspring of the second clutch was transferred to a new tube to start the next clonal generation. Third-generation offspring were used for RNA sequencing.

Third generation males were produced by placing second generation females in a medium containing 400 µM of methyl farnesoate (“MF”, Echelon Biosciences) just prior the production of their first clutch, as determined by well-swollen ovaries (Olmstead & Leblanc 2002). This ensured that the sex of their second clutch offspring was determined in the presence of MF. Otherwise, males were treated in the same way as described for the females. In particular, the newborn males were transferred back to standard medium, just as the third-generation females, and, throughout, culture media were exchanged daily for both males and females.

## Sampling

Just before the third-generation females released their first clutch, all individuals (males and females) were transferred individually to a 1.5 mL well on a culture plate, where they were kept for about three days before being sampled. Since RNA was extracted from whole individuals, no food was added during the last 12 h before sampling in order to minimize algal RNA contamination (most of which will be digested and hence degraded after 12 h). The period without food was kept relatively short to minimize starvation-dependent gene regulation.

To remove as much culture medium as possible, the individuals were blotted with absorbing paper (previously sterilized with UV radiation for 30 min), and then transferred to a 1.5 mL tube that was directly immersed in liquid nitrogen. Directly after the flash-freezing, three volumes of RNAlater ICE solution were added to preserve RNA, and samples were placed at -80°.

## RNA extraction, library preparation, and RNA-sequencing

The eight technical replicates (individuals of the same genotype and sex) were pooled, resulting in two samples (one male, one female) per

<figure><img src="image_2.png" /></figure>
biological replicate. Total RNA extraction and purification was carried out for each of the eight samples following the protocol of the *Daphnia* genomic consortium (DGC; DGC, Indiana University October 11, 2007), using Trizol Reagents and Qiagen RNEasy Mini Kit. The extracted and purified RNA samples were transferred to -80° and shipped on dry ice to the BSSE Genomic Sequencing Facilities, University of Basel, Switzerland.

Two lanes of cDNA corresponding to the eight biological replicates were constructed by the Department of *Biosystems* Science and Engineering (D-BSSE). The eight samples were individually labeled using TruSeq preparation kits. Each library (2 lanes) was sequenced on an Illumina NextSequation 500 sequencer with 76 cycles with the stranded paired-end protocol.

## Quality control and filtering

Read quality was assessed with FastQC v.0.10.1 (http://www.bioinformatics.babraham.ac.uk/projects/fastqc). Paired-end sequences were subjected to adapter trimming and quality filtering using trimmomatic v.0.36 (Bolger, Lohse & Usadel 2014). After trimming of adapter sequences, terminal bases with a quality score below three were removed from both ends of each read. Then, using the sliding window function and again moving in from both sides, further 4 bp-fragments were removed while their average quality score was below 15.

## Mapping and counting

Filtered reads were mapped to the publicly available *D. magna* genome assembly (NCBI database; Assembly name: daphmag2.4; GenBank assembly accession: GCA_001632505.1, Bioprojects accession: PRJNA298946), consisting of 28'801 scaffolds, 38'559 contigs and a total sequence length of 129'543'483 bp, as well as a genome annotation with 26'646 genes, using the RNA-Seq aligner STAR (Dobin et al. 2013) using default settings. The raw counts (number of mapped reads per transcript per sample) were obtained with the software featureCounts (Liao, Smyth & Shi 2014), a fast tool for counting mapped paired-end reads. Counts were summarized at the gene level using the annotation file.

## Differential gene expression

The analysis of differential gene expression was carried out with DESeq2 (version 1.10.1, Love, Huber & Anders 2014) implemented in R (R Core Team 2017). Raw read counts were used as input data, and the subsequent analyses used the normalizations of read counts as performed by DESeq2, which is currently considered best practice for the analysis of RNA-sequencing data (Conesa et al. 2016; Schurch et al. 2016). The males vs. females comparison was carried out with a two-factor design taking into account clone identity and sex. All *p*-values were adjusted for multiple testing with the Benjamini-Hochberg method, as implemented in DESeq2. Genes were considered differentially expressed (DE) if they had an adjusted *p*-value < 0.05 (False discovery rate, FDR = 5%). The degree of sex bias was determined by the fold-change (abbreviated FC) difference between the treatments. DE genes were classified into five groups: <twofold, >twofold, 2- to fivefold, 5- to 10-fold and >10-fold difference in expression (absolute changes rather than log-transformed changes). To obtain a broader overview of the expression profiles of the significantly DE genes, we performed a hierarchical clustering of the sex-DE genes, as implemented in DESeq2. The normalization in DESeq2 does not account for transcript length, hence it is possible that some differential exon usage (that could ultimately result in the existence of different isoforms) could be mistakenly interpreted as differential expression. However, because different transcripts

of most genes differ only weakly to moderately in size (Chern *et al.* 2006), normalization by transcript length (which has been criticized for other reasons, Dillies *et al.* 2013) would not strongly affect the inferred levels of fold-change in expression levels. Therefore, inferences of differential expression should typically be robust, at least in the class of genes with a greater than twofold change.

## Comparison of sex-biased gene expression With *D. galeata* and *D. pulex*

Protein sequences of *D. magna* (v2.4 GenBank: LRG00000000), *D. pulex* (version JGI060905: http://wfleabase.org/release1/current_release/fasta/) and *D. galeata* (http://server7.wfleabase.org/genome/Daphnia_galeata/) reference genomes were used as input in the software OrthoFinder (Emms & Kelly 2015), a fast method for inferring orthologous groups from protein sequences with enhanced accuracy. These correspond to 26'646, 30'940, and 33'555 protein sequences, respectively. We also used another software, OrthoMCL (Li, Stoeckert & Roos 2003) for inferring orthologs. Because the results were qualitatively and quantitatively similar, we present only the results of OrthoFinder here. For further analysis, we retained only those genes that were identified by OrthoFinder as single copy, one-to-one orthologs in all three species. We then compared this list with our list of sex-DE genes (adjusted *p* < 0.05), as well as the lists of sex-DE genes in *D. galeata* and *D. pulex* (Colbourne *et al.* 2011; Huylmans *et al.* 2016). The R package VennDiagram (Chen & Boutros 2011) was used for visualization of sex-DE genes in the three species.

We then focused on the core subset of 675 orthologous genes that were found to be consistently sex-DE in all three species and used BLAST2GO (version 4.1.9, (Conesa *et al.* 2016) to perform a functional annotation. The protein sequences of the *D. magna* genes in question were annotated using the NCBI nr database, allowing for 20 output alignments per query sequence with an *e*-value threshold of 0.001. The mapping and annotation steps implemented in BLAST2GO were run with default settings. Additionally, InterPro IDs from InterProScan were merged to the annotation to improve accuracy. Graphical representation of GO categories belonging to the 675 core-genes was obtained using the R package metacoder (Foster, Sharpton & Grünwald 2017). The BLAST2GO output files were searched for the terms "sex determination" (GO accession number 0007530), "sex differentiation" (0007548), "male sex differentiation" (0046661), "female gonad development" (0008585), "male gonad development" (0008584), "female sex determination" (0030237), "male sex determination" (0030238), and "female sex differentiation" (0046660). In addition, a universal list of genes involved in sex determination and sex differentiation pathways was established by searching the entire Genbank database for genes associated with the terms "sex determination/differentiation". After removal of redundancies, we obtained a list of 541 genes (hereafter referred to as the "NCBI list of genes"). We compared the annotations of our 675 core-genes (as determined by the BLAST2GO analysis) with this list to identify any shared genes. Finally, we performed an enrichment analysis of GO terms for the 675 core genes with consistent sex-biased expression in all three species, taking the 26'646 genes of *D. magna* as the reference GO composition. This was done using the GOatools Python script (https://github.com/tanghaibao/goatools), which performs Fisher's exact tests for differences in frequencies of GO-terms between the two lists (with Bonferroni correction). Enriched GO categories were summarized by a reduction of the complexity and level of GO terms (medium; allowed similarity = 0.7). UniprotKB was used to determine Gene Ontology Biological Process and Molecular Functions that were over-represented among genes DE between sexes and visualized with REVIGO (Supek *et al.* 2011).
## Data availability

Table S1 contains the initial numbers of raw reads and the percentages retained during the filtering steps. All genes found differentially expressed between sexes in *D. magna* are provided in Table S2. Table S3 contains data of sex-biased gene expression in the three *Daphnia* species investigated and their ortholog groups from OrthoFinder. Table S4 contains results of the annotated genes with sex-biased expression in all three species that were found in the NCBI list of genes involved in sex determination/differentiation. Raw Illumina sequence reads from this study have been deposited in the NCBI SRA database under the BioProject accession number PRJNA437676, and SRA accession number SRP134480. Supplemental material available at Figshare: https://doi.org/10.25387/g3.5965285.

## RESULTS

### Data quality

The RNA sequencing of the eight libraries resulted in a total of 1.59 billion raw reads. An average of 99.01% of raw reads passed the quality control. After end-trimming, an average of 93.03% aligned to the reference genome, resulting in an average of 81 million aligned reads per library, which constitutes a robust data basis for differential gene expression analyses. Table S1 shows the percentages of reads retained at each step in each of the samples.

### Sex-biased gene expression in *Daphnia magna*

We found a high number of genes that were DE between males and females with a total of 11'197 out of 26'646 genes being DE (adjusted *p* < 0.05), of which 8384 genes showed at least a twofold change (Table 1). The strong sexually dimorphic expression patterns can be visualized in the expression heatmap of the 8384 DE genes with more than two-fold expression difference between the sexes (Figure 1). Overall, a slight, but significant (*p* < 0.0001) majority of those genes were male-biased rather than female-biased (Table 1). This male-bias was found for all categories, except for the genes with a weak (< twofold) sex bias (Table 1). The list of all sex-biased genes can be found in supplemental data (Table S2).

### Comparisons of sex-biased gene expression among the three *Daphnia* species

The software OrthoFinder identified 7453 single copy, one-to-one orthologs present in all three species (Table 2). Among these, 5707 (76.5%) were sex-DE in at least one species, and 707 genes (9.5%) were sex-DE in all three species (Figure 2, Table 2). Only 32 of these 707 genes (4.5%) showed a different direction of bias in one of the species. The remaining 675 genes were biased in the same direction in all three species, and we therefore refer to these genes as the core-set of sex-DE genes in *Daphnia*. Among the genes of the core-set, 75% were female-biased (Figure 2), and genes with a strong expression-bias between sexes were more likely to be included in this core-set (Figure 3). Genes that showed significant sex-biased expression in only two out of the three species showed very similar patterns: a high proportion showed consistent bias (i.e., in the same direction in both species), and there was an excess of female-biased compared to male-biased

<figure><img src="image_3.png" /><figcaption>**Figure 1** Heatmap showing the normalized expression levels of the sex-DE genes (*p* < 0.05) with at least a two-fold expression difference between males and females. Each line represents a gene and each column specific sample (genotype and sex), with relative expression levels indicated by color (from highly overexpressed, red to highly underexpressed, blue, as indicated by the scale to the right). The dendrogram above the sample columns indicates clustering according to the Euclidean distance matrix implemented in the pheatmap R.</figcaption></figure>

genes (Figure 2). The excess of female-biased genes was even observed among the genes with sex-biased expression in only one species (Table S3). This was not only the case in *D. pulex* and *D. galeata*, for which an excess of female-biased genes had been reported earlier (Huylmans et al. 2016), but also in *D. magna*, where this result contrasts with the slight excess of male-biased genes found when all 26'646 genes were considered (as opposed to only the 7453 genes, for which single-copy, one-to-one orthologs could be identified in the other two species). The list of the 7453 orthologs, as well as the data on sex-biased expression for the three species is given in the supplementary data (Table S3).

### Functional analysis of the core-set of sex-DE genes

Among the core-set of 675 orthologous genes that were consistently sex-DE in all three species, 592 had an annotated function. The major GO categories of these genes are shown in Figure 4. This figure highlights that the largest fraction of genes belongs to the categories “cellular process”, “metabolic process”, “single organism process” and “biological regulation”. The results of the GO enrichment analysis are shown in Figure 5. Enriched terms are linked to “RNA binding” processes,

<table><caption>**Table 1** Number of significantly (adjusted *p* < 0.05) sex-DE genes in *Daphnia magna*, for different degrees of bias as well as percentage of genes with male-biased expression</caption><thead><tr><th></th><th>All</th><th>&lt; twofold</th><th>2- to fivefold</th><th>5- to 10-fold</th><th>&gt; 10-fold</th><th>&gt; twofold</th></tr></thead><tbody><tr><th>Number of sex-DE genes</th><td>11197</td><td>2813</td><td>6005</td><td>1409</td><td>970</td><td>8384</td></tr><tr><th>Percentage male-biased</th><td>53%</td><td>38%</td><td>57%</td><td>64%</td><td>53%</td><td>58%</td></tr></tbody></table>

<figure><img src="image_4.png" /></figure>
<table><caption>Table 2 Number of the one-to-one orthologous genes found by OrthoFinder in the three *Daphnia* species</caption><thead><tr><th>Species</th><th>Number of protein sequences in the genome</th><th>Number of Sex-biased genes<sup>a</sup></th><th>Proportion of sex-biased genes</th><th>Number of 1-to-1 orthologs</th><th>Number of sex-biased 1-to-1 orthologs</th><th>Proportion of sex-biased 1-to-1 orthologs</th></tr></thead><tbody><tr><td><i>D. magna</i></td><td>26646</td><td>11192</td><td>0.42</td><td>7453</td><td>4818</td><td>0.65</td></tr><tr><td><i>D. pulex</i></td><td>30940</td><td>6393</td><td>0.21</td><td>7453</td><td>2812</td><td>0.38</td></tr><tr><td><i>D. galeata</i></td><td>33555</td><td>5842</td><td>0.17</td><td>7453</td><td>1384</td><td>0.19</td></tr></tbody></table>

known to play a key role in post-transcriptional gene regulation (Glisovic et al. 2008; Cléry & Allain 2013) and, more generally, terms linked to “RNA”. Of the 592 genes only one gene has a GO term linked to sex determination or sex differentiation (which is neither more nor less than expected by chance): The gene “peptidyl-prolyl cis-trans isomerase FKBP4” has a female-biased expression in all three species and its GO term includes “male sex differentiation”. Moreover, among the same 592 core genes with a functional annotation, 14 were listed in the NCBI list of genes known to be involved in sex determination or differentiation pathways in other species (Table S4).

DISCUSSION

Sex-biased gene expression in *Daphnia magna*

We found a very high number of genes being DE between males and females in *D. magna*. This result is largely congruent with the previous studies on other *Daphnia* species (Colbourne et al. 2011; Huylmans et al. 2016), but the overall number and proportion of genes that show sex-biased expression is higher in *D. magna* than in the other two species. Indeed, the proportion of genes with sex-biased expression in *D. magna* is similar to that reported for species with genetic sex determination (GSD). For instance, two early, microarray-based studies on *Drosophila melanogaster* found twofold or greater expression differences between sexes in 30–40% of all genes (Parisi et al. 2004, Innocenti & Morrow 2010). A more recent study based on RNA-sequencing found that about two-third of genes showed sex-biased expression in multiple *Drosophila* species (genes with a less than twofold change included). In *D. magna*, the proportion among all genes is 42% (31% with a greater than twofold change). However, considering only the set of single-copy, orthologous genes (which likely contain a lower proportion of annotation errors, see below), the proportion is 65% (41% with a greater than twofold change). Our results therefore suggest that these two systems, which are comparable in terms of body size and due to the fact that whole, adult animals were sampled, have similar proportions of sex-biased genes. More studies on comparable ESD-GSD species pairs will have to be investigated to determine the generality of this conclusion.

A recent study on *D. pulex* found that the annotation of the genome used here (Colbourne et al. 2011) likely contained a non-negligible fraction of pseudogenes or other false positives (Ye et al. 2017), and that the number of genes in the genome may be closer to 20'000 than the initially estimated ~31'000. It is possible that the current estimates of the total number of genes in the *D. magna* (~27'000) and *D. galeata* (~34'000) genomes may also be overestimates. Using these genomes as reference in a differential gene expression analysis may have affected the estimate of the proportion of differentially expressed genes only if the proportion of falsely annotated genes is different among the DE genes than among the non-DE genes. It is unclear if any such bias exists. However, two independent lines of evidence suggest that, if anything, such a bias has led to an underestimation of the proportion of sex-DE

genes: First, the proportion of sex-DE genes was higher among the single-copy orthologs (which less likely contain annotation error) than among all genes in two out of three *Daphnia* species. Second, *D. pulex* and *D. galeata*, which both have higher estimated number of genes than *D. magna*, have a lower estimated proportion (~20%) of sex-DE genes. Yet, the differences in the proportions of sex-DE genes among the three *Daphnia* species may also be explained by differences in methodology and statistical power. The *D. pulex* data (Colbourne et al. 2011) are based on a microarray study, a methodology known to be less sensitive for lowly-expressed genes than RNA-sequencing (Harrison, Wright & Mank 2012). The *D. galeata* study (Huylmans et al. 2016) was based on RNA-sequencing but used only two clonal lines as biological replicates. As mentioned by Huylmans et al. (2016), this may have led to a rather low statistical power to detect sex-DE genes, especially among the considerable number of genes that showed expression differences between the two clones.

When comparing the proportion of genes with sex-biased expression with other studies, it is important to remember that we performed RNA-sequencing on whole animals and hence included all tissues present at that the time of sampling in adult males and females. Patterns of sex-specific gene expression are known to be tissue-specific in many cases (Ellegren & Parsch 2007; Toyota et al. 2017), with strongest differences being found in brains (at least in mammals) and, unsurprisingly, in gonad tissues (Mank 2009). Thus, it is difficult to compare our results with those in larger animals, where studies have mostly been carried out on specific tissues (e.g., 54.5% of genes were found to be sex-DE in *Mus musculus* liver, (Yang et al. 2006)). Moreover, genes that do not show sex-biased expression may still differ in their expression patterns among tissues (Yang et al. 2006). Hence, when sampling whole animals, some of these genes may be identified as sex-biased because different tissues may occur in different proportions in males vs. females, for instance due to anatomical differences between sexes. It is possible that a part of the genes that were found to be sex-biased in our study and in other studies based on whole animals (e.g., *Drosophila*) are explained by such effects (i.e., sex-biased expression of these genes may be a consequence rather than a cause of the phenotypic differences between sexes (Mank 2017)).

Another factor that potentially contributes to an overestimation of the number of sex-DE genes is the fact that males were produced by artificially treating their mothers with the juvenile hormone analog MF (Huylmans et al. 2016). However, in a separate study we found that MF exposure changes expression levels of a much lower number of genes (only a few 100s) than were DE between sexes (Molinier et al., in prep). Moreover, the males used in our experiment were exposed to MF only for three days when they were still oocytes inside the ovaries of their mothers up to one day after they were released from the brood pouch. It is thus highly unlikely that large proportions of the sex-DE genes are in fact explained by the effects of early MF exposure and would not have shown up in naturally produced males (which also involves exposure to a natural juvenile hormone produced by their mother).

<footnote><sup>a</sup>As identified in the single species studies (current study, Colbourne et al. 2011, Huylmans et al. 2016).</footnote>
<figure><img src="image_5.png" /><figcaption>B</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>C</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 2** A. Venn diagram showing the number of sex-DE genes among the 7453 one-to one orthologs in each of the three species of *Daphnia*. B. Number of sex-DE genes being biased in the same vs. opposite directions. Panels from left to right: genes being sex-DE all three species (707 genes), genes being sex-DE only in *D. magna* and *D. pulex* (1402 genes), *D. magna* and *D. galeata*, (349 genes), and *D. pulex* and *D. galeata* (142 genes). C. Number of female-biased and male-biased genes in each of the four categories depicted in Figure 2 B. (only genes being biased in the same direction).</figcaption></figure>

<figure><img src="image_8.png" /></figure>
<figure><img src="image_9.png" /><figcaption>**Figure 3** Proportion of genes with different degrees of sex-bias (the degree of sex-bias is summarized in four categories of fold change). Panels from left to right: genes being sex biased in all three species (707 genes), in two species (1751 genes), and only in *D. magna* (2360 genes). Only the 7453 genes, for which single-copy, one-to-one orthologs could be identified in all three species were considered for this analysis.</figcaption></figure>

Contrary to previous findings in *D. galeata* and *D. pulex* (which both show a clear excess of female-biased genes over male-biased genes), we observed a slight excess of male-biased genes in *D. magna*. It is currently difficult to say whether this difference between studies reflects a biological reality (i.e., difference between species within the genus) or whether it may be explained by some methodological differences between the studies. Interestingly, however, a strong excess of female-biased genes was recovered also in *D. magna* when investigating the subset of single-copy genes for which one-to-one orthologs could be identified in the other species. The excess of female-biased genes might be a general feature in the genus *Daphnia*, at least for single-copy genes that are sufficiently conserved for orthologs to be identified across the major sub-clades of the genus. On the other hand, the strong excess of female-biased single-copy one-to-one orthologs, suggests that the remaining genes show a substantial excess of male-biased expression, at least in *D. magna*. The remaining genes likely include many paralogs and other less conserved genes, for which the identification of orthologs is difficult. Hence, different evolutionary rates of genes with male-biased vs. female-biased expression could drive the observed patterns. Genes with male-biased expression evolve faster than female-biased genes in *D. pulex* and *D. galeata*, as well as in Drosophila, *Caenorhabditis*, and several mammal species (see reviews by Ellegren & Parsch 2007 and Parsch & Ellegren 2013). Faster evolution of genes with male-biased expression might be explained by positive selection being more common in these genes (Ellegren & Parsch 2007), and may lead to ortholog identification being more difficult in these genes, which may explain or at least contribute to the observed difference between the two sets of genes in *D. magna*.

## Comparisons of sex-biased gene expression among the three *Daphnia* species

The 7453 single-copy, one-to-one orthologs present in all three species represent less than 30% of all genes used in the *D. magna* analysis. However, as pointed out above, the number of genes predicted by the current genome annotation used in this study may in fact be a rather strong overestimation of the true number of genes present in the species. Second, the software identified a considerably higher number of orthogroups, which, however, also include paralogs. We decided to

restrict the analysis to single-copy, one-to-one orthologs because interpretation of expression patterns in paralogs is much less straightforward. For instance, if a sex-DE gene is single-copy in one species but has two paralogs in another species due to duplication after speciation, it is difficult to say which one of the two genes is more homologous in function, and hence should also be sex-DE if the gene belongs to the core-set of genes with sex-biased expression in the genus (Koonin 2005).

Only a low percentage of the genes found to be sex-biased in two and especially in three species were biased in different directions. Moreover, genes with a higher expression bias were more likely to be sex-DE in two or three species than just in one. While the latter observation may in part be explained by issues of statistical power (genes with a low degree of sex-bias have a lower probability to be detected), both observations nonetheless suggest that the 675 orthologous genes that were found to be consistently sex-DE in all three species indeed represent a robust core-set of sex-biased genes in the genus *Daphnia*. It is likely that some genes that were DE between a pair of species but not in all species should also have been included in this core-set, as differential expression may have been non-significant in one species just due to a lack of statistical power. Indeed, the three studies differ in methodology (microarray vs. RNA-Sequencing), number of biological replicates, aspects of data analysis, etc. (see above). In addition, the quality of the genome assemblies and annotations used to analyze these data may also differ between species. These differences may also explain some of the between-species differences in the number and proportion of sex-DE genes.

## Functional analysis of the core-set of sex-DE genes

Our study identified a core-set of genes for which sex-biased expression is probably conserved in the genus *Daphnia*. Hence these genes may play a fundamental role in determining and maintaining male vs. female phenotypes in this genus with environmental sex determination. The functional analysis of these genes identified the gene "peptidyl-prolyl cis-trans isomerase FKBP4" with GO term "male sex differentiation", an immunophilin protein with peptidylprolyl isomerase and co-chaperone activities. It is a component of steroid receptors heterocomplexes and may play a role in
<figure><img src="image_10.png" /><figcaption>**Figure 4** Composition and hierarchical organization of GO terms associated to the 675 genes with consistent sex-biased gene expression in all three species. nobs: number of genes in the given GO category.</figcaption></figure>

the intracellular trafficking of hetero-oligomeric forms of steroid hormone
receptors between cytoplasm and nuclear compartments. Steroid receptors
initiate the signal transduction for steroid hormones, including sexual
hormones such as estrogen and androgen (e.g., Voigt et al. 2009). Their
role in sex dimorphism, also in species with environmental sex determi-
nation, thus makes sense, though the role of this particular gene has not yet
been investigated in *Daphnia*. The 14 additional genes whose functional
annotation matched of the descriptors of the genes on the NCBI list of
genes involved in sex determination/differentiation pathways may repre-
sent further fundamental genes involved in sex determination or sex dif-
ferentiation in *Daphnia*. They contain functions known to be implicated in
sexual development in other species, such as the "Beta-catenin 1" which is
a key transcriptional regulator of the canonical Wnt-signaling pathway,
known to be implicated in female reproductive development in mammals
(Bernard & Harley 2007, Liu, Bingham & Parker 2008). We also found the
"fibroblast growth factor receptor", receptor of Fgfs (fibroblast growth
factors), whose function may be involved in sex determination and re-
productive system development in many species and appears to be highly
conserved (Colvin et al. 2001). Finally, the gene "ovarian tumor" is also
known to be involved and required in the determination of the sexual
identity of female germ cells (Pauli, Oliver & Mahowald 1993). Our results

suggest that these genes are involved in maintaining phenotypic differences
between sexes, also at the adult stage, at last in *Daphnia*. A large proportion
of the identified core genes have unknown function. Therefore, we need
more functional annotations, especially also on more closely related spe-
cies, as well as sex-specific expression data on earlier developmental stages
before we can obtain a clear mechanistic picture of sex determination and
sex differentiation in *Daphnia* (Kato *et al.* 2011).

In conclusion, our study provides data on sex-biased gene expression
for the model organism *D. magna* and for *Daphnia* in general, specif-
ically by identifying a core-set of sex-DE genes in all three major sub-
clades of the genus. More generally, our results suggest that the
proportion of genes with sex-biased expression in ESD species is not
lower than in species with GSD.

ACKNOWLEDGMENTS

We thank the Zoo of Moscow and N. I. Skuratov for sampling permits,
and David Frey for help with maintenance of indoors cultures. We
thank the Department of Biosystem Science and Engineering of the
ETH Zurich, in particular C. Beisel and I. Nissen for Illumina
sequencing, and we gratefully acknowledge support by M-P. Dubois,

<figure><img src="image_11.png" /></figure>
<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 5** Enrichment analysis of GO terms among the core-set of genes compared to the entire list of *D. magna* genes. Shown are over-represented (Bonferroni-corrected *p* < 0.05) terms and functional categories (the latter distinguished by color) for molecular functions (A) and biological processes (B). The size of each rectangle is proportional to the -log(p-value) for its category.</figcaption></figure>

the platform Service des Marqueurs Génétiques en Ecologie at CEFE, and the genotyping and sequencing facilities of the Institut des Sciences de l'Evolution-Montpellier and the Labex Centre Méditerranée Environnement Biodiversité (CeMEB). We thank Mathilde Cordellier and John Colbourne for sharing data and information and two anonymous reviewers for their comments. We thank the University of Fribourg, the Montpellier Bioinformatic Biodiversity plateform and the Labex CeMEB for access to high-performance computing clusters. The sequence data for the *D. magna* genome

project V2.4 were produced by The Center for Genomics and Bioinformatics at Indiana University and distributed via wFleaBase in collaboration with the *Daphnia* Genomics Consortium (project supported in part by NIH award 5R24GM078274-02 "*Daphnia* Functional Genomics Resources"). Our work was supported by the Swiss National Science Foundation (Grant no. 31003A_138203), the Russian Foundation of Basic Research, and the European Union (Marie Curie Career Integration Grant PCIG13-GA-2013-618961, DamanNMP). YG was supported by RFBR grant 16-04-01579.
LITERATURE CITED

Adamowicz, S., A. Petrusek, and J. Colbourne, 2009 The scale of divergence: a phylogenetic appraisal of intercontinental allopatric speciation in a passively dispersed freshwater zooplankton genus. Mol. Phylogenet. Evol. 50: 423-436. https://doi.org/10.1016/j.ympev.2008.11.026

Bergero, R., and D. Charlesworth, 2009 The evolution of restricted recombination in sex chromosomes. Trends Ecol. Evol. 24(2): 94-102. https://doi.org/10.1016/j.tree.2008.09.010

Bernard, P., and V. Harley, 2007 Wnt4 action in gonadal development and sex determination. The international journal of biochemistry & cell. https://doi.org/10.1016/j.biocel.2006.06.007

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

Bull, J. J., 1985 Sex determining mechanisms: An evolutionary perspective. Experientia 41(10): 1285-1296. https://doi.org/10.1007/BF01952071

Chen, H., and P. Boutros, 2011 VennDiagram: a package for the generation of highly-customizable Venn and Euler diagrams in R. BMC Bioinformatics 12: 35. https://doi.org/10.1186/1471-2105-12-35

Chern, T. M., E. Van Nimwegen, C. Kai, J. Kawai, P. Carninci et al., 2006 A simple physical model predicts small exon length variations. PLoS Genet. 2: 606-613.

Cléry, A., and F. Allain, 2013 From structure to function of rna binding domains.

Colbourne, J. & Hebert, P., 1997 Evolutionary origins of phenotypic diversity in Daphnia. Molecular evolution and.

Colbourne, J. K., M. E. Pfeffer, D. Gilbert, W. K. Thomas, A. Tucker et al., 2011 The ecorespective genome of Daphnia pulex. Science 331: 555-561. https://doi.org/10.1126/science.1197761

Colvin, J., R. Green, J. Schmahl, B. Capel, and D. Ornitz, 2001 Male-to-female sex reversal in mice lacking fibroblast growth factor 9. Cell.

Conesa, A., P. Madrigal, S. Tarazona, D. Gomez-Cabrero, A. Cervera et al., 2016 A survey of best practices for RNA-seq data analysis. Genome Biol. 17(1): 13 (erratum: Genome Biol. 17: 181). https://doi.org/10.1186/s13059-016-0881-8

Dillies, M. A., A. Rau, J. Aubert, C. Hennequet-Antier, M. Jeanmougin et al., 2013 A comprehensive evaluation of normalization methods for Illumina high-throughput RNA sequencing data analysis. Brief. Bioinform. 14(6): 671-683. https://doi.org/10.1093/bib/bbs046

Dobin, A., C. A. Davis, F. Schlesinger, J. Drenkow, C. Zaleski et al., 2013 STAR: ultrafast universal RNA-seq aligner. Bioinformatics 29(1): 15-21. https://doi.org/10.1093/bioinformatics/bts635

Dukić, M., D. Berner, M. Roesti, C. R. Haag, and D. Ebert, 2016 A high-density genetic map reveals variation in recombination rate across the genome of Daphnia magna. BMC Genet. 17(1): 137. https://doi.org/10.1186/s12863-016-0445-7

Eads, B. D., J. K. Colbourne, E. Bohuski, and J. Andrews, 2007 Profiling sex-biased gene expression during parthenogenetic reproduction in Daphnia pulex. BMC Genomics 8(1): 464. https://doi.org/10.1186/1471-2164-8-464

Ellegren, H., and J. Parsch, 2007 The evolution of sex-biased genes and sex-biased gene expression. Nat. Rev. Genet. 8(9): 689-698. https://doi.org/10.1038/nrg2167

Emms, D., and S. Kelly, 2015 OrthoFinder: solving fundamental biases in whole genome comparisons dramatically improves orthogroup inference accuracy. Genome Biol. 16(1): 157. https://doi.org/10.1186/s13059-015-0721-2

Foster, Z. S. L., T. J. Sharpton, and N. J. Grünwald, 2017 Metacoder: an R package for visualization and manipulation of community taxonomic diversity data. PLoS Comput Biol. 13: e1005404. https://doi.org/10.1371/journal.pcbi.1005404

Galimov, Y., B. Walser, and C. R. Haag, 2011 Frequency and inheritance of non-male producing clones in Daphnia magna: Evolution towards sex specialization in a cyclical parthenogen? J. Evol. Biol. 24(7): 1572-1583. https://doi.org/10.1111/j.1420-9101.2011.02288.x

Giraudo, M., M. Douville, G. Cottin, and M. Houde, 2017 Transcriptomic, cellular and life-history responses of Daphnia magna chronically exposed to benzotriazoles: Endocrine-disrupting potential and molting effects. PLoS One, 12: e0171763. https://doi.org/10.1371/journal.pone.0171763

Glisovic, T., J. Bachorik, J. Yong, and G. Dreyfuss, 2008 RNA-binding proteins and post-transcriptional gene regulation. FEBS Lett. 582(14): 1977-1986. https://doi.org/10.1016/j.feb.2008.03.004

Gorbi, G., F. Moroni, S. Sei, and V. Rossi, 2011 Anticipatory maternal effects in two different clones of Daphnia magna in response to food shortage. J. Limnol. 70(2): 222. https://doi.org/10.4081/jlimnol.2011.222

Grath, S., and J. Parsch, 2016 Sex-Biased Gene Expression. Annu. Rev. Genet. 50(1): 29-44. https://doi.org/10.1146/annurev-genet-120215-035429

Harrison, P. W., A. E. Wright, and J. E. Mank, 2012 The evolution of gene expression and the transcriptome-phenotype relationship. Semin. Cell Dev. Biol. 23(2): 222-229. https://doi.org/10.1016/j.semcdb.2011.12.004

Herrmann, M., S. P. Ravindran, K. Schwenk, and M. Cordellier, 2018 Population transcriptomics in Daphnia: the role of thermal selection. Mol. Ecol. https://doi.org/10.1111/mec.14450

Hobaek, A., and P. Larsson, 1990 Sex Determination in Daphnia Magna. Ecology 71(6): 2255-2268. https://doi.org/10.2307/1938637

Huymans, A. K., A. L. Ezquerra, J. Parsch, and M. Cordellier, 2016 De novo transcriptome assembly and sex-biased gene expression in the cicalay parthenogenetic daphnia galeata. Genome Biol. Evol. 8(10): 3120-3139. https://doi.org/10.1093/gbe/evw221

Innocenti, P., and E. H. Morrow, 2010 The sexually antagonistic genes of Drosophila melanogaster. PLoS Biol. 8: e1000335. https://doi.org/10.1371/journal.pbio.1000335

Ishida, S., A. A. Kotov, and D. J. Taylor, 2006 A new divergent lineage of Daphnia (Cladocera: Anomopoda) and its morphological and genetical differentiation from Daphnia curvirostris Eylmann, 1887. Zool. J. Linn. Soc. 146: 385-405.

Kato, Y., K. Kobayashi, H. Watanabe, and T. Iguchi, 2011 Environmental sex determination in the branchiopod crustacean Daphnia magna: deep conservation of a Doublesex gene in the sex-determining pathway. PLoS Genet. 7: e1001345. https://doi.org/10.1371/journal.pgen.1001345

Klütten, B., U. Dülmer, M. Engels, and H. Ratte, 1994 ADaM, an artificial freshwater for the culture of zooplankton. Water Res. 28(3): 743-746. https://doi.org/10.1016/0043-1354(94)90157-0

Koonin, E. V., 2005 Orthologs, Paralogs, and Evolutionary Genomics. Annu. Rev. Genet. 39(1): 309-338. https://doi.org/10.1146/annurev.genet.39.073003.114725

Korpelainen, H., 1990 Sex ratios and conditions required for environmental sex determination in animals. Biol. Rev. Camb. Philos. Soc. 65(2): 147-184. https://doi.org/10.1111/j.1469-185X.1990.tb01187.x

Li, L., C. J. Stoeckert, and D. S. Roos, 2003 OrthoMCL: identification of ortholog groups for eukaryotic genomes. Genome Res. 13(9): 2178-2189. https://doi.org/10.1101/gr.1224503

Liao, Y., G. K. Smyth, and W. Shi, 2014 FeatureCounts: an efficient general purpose program for assigning sequence reads to genomic features. Bioinformatics 30(7): 923-930. https://doi.org/10.1093/bioinformatics/btt656

Liu, C., Bingham, N. & Parker, K., 2008 Sex-specific roles of β-catenin in mouse gonadal development. Human molecular. https://doi.org/10.1093/hmg/ddn362

Love, M. I., W. Huber, and S. Anders, 2014 Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 15(12): 550. https://doi.org/10.1186/s13059-014-0550-8

Lynch, M., R. Gutenkunst, M. Ackerman, K. Spitze, Z. Ye et al., 2017 Population genomics of Daphnia pulex. Genetics 206: 315-332.

Mank, J. E., 2009 Sex Chromosomes and the Evolution of Sexual Dimorphism: Lessons from the Genome. Am. Nat. 173(2): 141-150. https://doi.org/10.1086/595754

Mank, J. E., 2017 The transcriptional architecture of phenotypic dimorphism. Nat. Ecol. Evol. 1(1): 6. https://doi.org/10.1038/s41559-016-0006

Miner, B. E., L. De Meester, M. E. Pfeffer, W. Lampert, and N. G. Hairston, 2012 Linking genes to communities and ecosystems: Daphnia as an ecogenomic model. Proc. Biol. Sci. 279(1735): 1873-1882. https://doi.org/10.1098/rspb.2011.2404

Olmstead, A. W., and G. A. Leblanc, 2002 Juvenile hormone methyl farnesoate is a sex determinant in the crustacean Daphnia magna. J. Exp. Zool. 293(7): 736-739. https://doi.org/10.1002/jez.10162

<figure><img src="image_14.png" /></figure>
Orsini, L., D. Gilbert, R. Podicheti, M. Jansen, J. B. Brown *et al.*, 2016 *Daphnia magna* transcriptome by RNA-Seq across 12 environmental stressors. *Sci. Data* 3: 160030. https://doi.org/10.1038/sdata.2016.30

Parisì, M., R. Nuttall, P. Edwards, J. Minor, D. Naiman *et al.*, 2004 *A survey of ovary-, testis-, and soma-biased gene expression in Drosophila melanogaster adults. Genome Biol.* 5(6): R40. https://doi.org/10.1186/gb-2004-5-6-r40

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

Pauli, D., B. Oliver, and A. Mahowald, 1993 *The role of the ovarian tumor locus in Drosophila melanogaster germ line sex determination. Development* 119: 123–134.

R Core Team, 2017 *R: A language and environment for statistical computing*. R Foundation for Statistical Computing, Vienna, Austria. URL https://www.R-project.org/.

Radhakrishnan, S., R. Listerman, J. Neuwald, A. Severin, and N. Valenzuela, 2017 *Transcriptomic responses to environmental temperature by tur-tles with temperature-dependent and genotypic sex determination assessed by RNAseq inform the genetic architecture of embryonic gonadal development*. *PLoS One*, 12: e0172044. https://doi.org/10.1371/journal.pone.0172044

Reisser, C. M. O., D. Fasel, E. Hurlimann, M. Dukic, C. Haag-Liautard *et al.*, 2016 *Transition from environmental to partial genetic sex determina-tion in Daphnia through the evolution of a female-determining incipient W-chromosome. bioRxiv*.

Roulin, A. C., J. Routtu, M. D. Hall, T. Janicke, I. Colson *et al.*, 2013 *Local adaptation of sex induction in a facultative sexual crustacean: insights from QTL mapping and natural populations of Daphnia magna*. Mol. Ecol.* 22(13): 3567–3579. https://doi.org/10.1111/mec.12308

Routtu, J., M. D. Hall, B. Albere, C. Beisel, R. Bergeron *et al.*, 2014 *An SNP-based second-generation genetic map of Daphnia magna and its appli-cation to QTL analysis of phenotypic traits*. BMC Genomics 15(1): 1033. https://doi.org/10.1186/1471-2164-15-1033

Scourfield, D. J., and J. P. Harding, 1966 *A Key to the British Freshwater Cladocera, with Notes on their Ecology*, pp. 55, Ed. 3. Science Publication No. 5. Freshwater biological Association.

Schurch, N. J., P. Schofield, M. Gierliński, C. Cole, A. Sherstnev *et al.*, 2016 *How many biological replicates are needed in an RNA-seq ex-periment and which differential expression tool should you use? RNA* 22(6): 839–851. https://doi.org/10.1261/rna.053959.115

Shoemaker, C., M. Ramsey, J. Queen, and D. Crews, 2007 *Expression of Sox9, Mis, and Dmrt1 in the gonad of a species with temperature-dependent sex determination*. *Dev. Dyn.* 236(4): 1055–1063. https://doi.org/10.1002/dvdy.21096

Spanier, K. I., M. Jansen, E. Decaestecker, G. Hulsmans, D. Becker *et al.*, 2017 *Conserved Transcription Factors Steer Growth-Related Genomic Programs in Daphnia. Genome Biol. Evol.* 9(6): 1821–1842. https://doi.org/10.1093/gbe/evx127

Supek, F., M. Bošnjak, N. Škunca, T. Šmuc, and C. O'Donovan, 2011 *REVIGO summarizes and visualizes long lists of gene ontology terms*. *PLoS One* 6: e21800. https://doi.org/10.1371/journal.pone.0021800

Svendsen, N., C. M. O. Reisser, M. Duki, V. Thuillier, A. Segard *et al.*, 2015 *Uncovering Cryptic Asexuality in Daphnia magna by RAD Sequencing*. *Genetics* 201(3): 1143–1155. https://doi.org/10.1534/genetics.115.179879

Torres Maldonado, L., A. Landa Piedra, N. Moreno Mendoza, A. Marmolejo Valencia, A. Meza Martínez *et al.*, 2002 *Expression profiles of Dax1, Dmrt1, and Sox9 during temperature sex determination in gonads of the sea turtle Lepidochelys olivacea*. *Gen. Comp. Endocrinol.* 129(1): 20–26. https://doi.org/10.1016/S0016-6480(02)00511-7

Toyota, K., T. D. Williams, T. Sato, N. Tatarazako, and T. Iguchi, 2017 *Comparative ovarian microarray analysis of juvenile hormone-responsive genes in water flea Daphnia magna: potential targets for toxicity*. *J. Appl. Toxicol.* 37(3): 374–381. https://doi.org/10.1002/jat.3368

Voigt, C., G. Ball, and J. Balthazart, 2009 *Sex differences in the expression of sex steroid receptor mRNA in the quail brain*. *J. Neuroendocrinol.* 21(12): 1045–1062. https://doi.org/10.1111/j.1365-2826.2009.01929.x

Wijchers, P. J., and R. J. Festenstein, 2011 *Epigenetic regulation of autosomal gene expression by sex chromosomes*. *Trends Genet.* 27(4): 132–140. https://doi.org/10.1016/j.tig.2011.01.004

Wright, A. E., M. Fumagalli, C. R. Cooney, N. I. Bloch, F. G. Vieira *et al.*, 2018 *Male-biased gene expression resolves sexual conflict through the evolution of sex-specific genetic architecture*. *Evolution Letters*.

Xu, S., M. S. Ackerman, H. Long, L. Bright, K. Spitze *et al.*, 2015 *A Male-Specific Genetic Map of the Microcrustacean Daphnia pulex Based on Single-Sperm Whole-Genome Sequencing*. *Genetics* 201(1): 31–38. https://doi.org/10.1534/genetics.115.179028

Yang, X., E. E. Schadt, S. Wang, H. Wang, A. P. Arnold *et al.*, 2006 *Tissue-specific expression and regulation of sexually dimorphic genes in mice*. *Genome Res.* 16(8): 995–1004. https://doi.org/10.1101/gr.5217506

Yatsu, R., S. Miyagawa, S. Kohno, B. B. Parrott, K. Yamaguchi *et al.*, 2016 *RNA-seq analysis of the gonadal transcriptome during Alligator mississippiensis temperature-dependent sex determination and differentiation*. *BMC Genomics* 17(1): 77. https://doi.org/10.1186/s12864-016-2396-9

Ye, Z., S. Xu, K. Spitze, J. Asselman, and X. Jiang, 2017 *A new reference genome assembly for the microcrustacean Daphnia pulex*. *G3: Genes, Genomes*. *Genetics* 7(5): 1405–1416. https://doi.org/10.1534/g3.116.038638

Communicating editor: B. Andrews