# Evolution of Gene Expression during a Transition from Environmental to Genetic Sex Determination

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

► To cite this version:

Cécile Molinier, Céline Reisser, Peter Fields, Adeline Ségard, Yan Galimov, et al.. Evolution of Gene Expression during a Transition from Environmental to Genetic Sex Determination. Molecular Biology and Evolution, 2019, 36 (7), pp.1551-1564. 10.1093/molbev/msz123. hal-02356483

HAL Id: hal-02356483

https://hal.science/hal-02356483

Submitted on 16 Jun 2021

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
chromosomes (Peichel et al. 2001; Lindholm and Breden 2002; Khil et al. 2004; Bachtrog et al. 2008). Yet, many details, especially regarding the early stages of sex chromosome evolution, are still unknown or controversial. The molecular mechanisms and pathways underlying the different forms of sex determination are just starting to be elucidated (Bachtrog et al. 2014; Herpin and Schartl 2015). To investigate these issues, many studies have focused on intermediate systems, which are rather common between hermaphroditism and dioecy (Charlesworth 2006; Barrett 2010; Weeks et al. 2014) but rare and understudied for the ESD to GSD transition.

Daphnia magna is a freshwater crustacean which offers unique opportunities to study the evolutionary transition from ESD to GSD because it contains mixed populations with some individuals being genetically determined females and others having ancestral ESD (Galimov et al. 2011). Females are determined by a single nuclear chromosome region (Galimov et al. 2011; Reisser et al. 2017). The *D. magna* system thus closely fits the scenario outlined above for the consensus model. Like most other *Daphnia* species, *D. magna* reproduces by cyclical parthenogenesis, a partly sexual and partly asexual life cycle (Hebert 1978; Ebert 2005). During asexual reproduction, females produce clonal, live-born offspring, whereas sexual reproduction leads to the production of diapause eggs which are required to survive periods of freezing or drought (Cáceres 1998; Ebert 2005). In ESD individuals, the sex of the clonally produced offspring is determined by the environment (Kleiven et al. 1992). Specifically, male production is induced by a juvenile hormone emitted by the mother in response to specific conditions, such as shortened photoperiod and/or increased population density (Olmstead and Leblanc 2002; Roulin et al. 2013). Male production can also be experimentally induced by adding the crustacean juvenile hormone methyl farnesoate (MF) to the culture medium at a precise moment of the ovarian cycle (Olmstead and Leblanc 2002). The individuals with ESD are also called “male producers” (“MP”) for their ability to produce males during clonal reproduction. In contrast to these MP clones, some other clones cannot produce males, neither in nature nor under artificial hormone exposure (Galimov et al. 2011), and not even at very low frequencies (Svendsen et al. 2015). These clones, called “non-male producers” (“NMP”), are genetically determined females who participate in sexual reproduction only as females. Diapause eggs are still produced sexually (NMP are thus not obligate parthenogens, which also occur in some *Daphnia* species [Lynch et al. 2008]). Hence, to be able to undergo successful diapause egg production, they need to coexist with MP clones who will provide the males necessary for sexual reproduction.

Previous studies have shown that the MP/NMP polymorphism in *D. magna* is determined by a large (~2 Mb), non-recombining chromosome region on linkage group 3 (LG 3) (Reisser et al. 2017). This “NMP chromosome region” segregates as a single Mendelian locus with a dominant female-determining allele on an incipient W chromosome. Heterozygous individuals (genotype WZ) are genetically determined females (NMP), whereas homozygotes (ZZ) are MP

phenotypes with ESD, and crosses between NMP females and MP males result in 50% NMP and 50% MP offspring (Galimov et al. 2011; Reisser et al. 2017). Hence, MP and NMP do not differ genetically, except for the NMP region, for which $F_{ST}$ is 0.18 between MP and NMP individuals (Reisser et al. 2017). Although the inheritance mode and general genetic architecture of the NMP trait have been elucidated, we still do not know which of the over 600 genes in the NMP chromosome region are involved in producing the NMP phenotype, nor do we know the identity of downstream genes and molecular networks involved in the difference between the MP and NMP phenotypes.

To investigate the genes and networks potentially involved in the expression of the NMP phenotype, we sequenced mRNA extracted from MP (ESD) and NMP (GSD) females, both under control conditions (normal culturing conditions) and after hormone exposure. Whole adult females including the maturing oocytes in their ovaries (but not the offspring in their brood pouches) were sampled at the moment of greatest hormone-sensitivity for sex determination, as it is unclear whether the primary response to the hormone occurs in oocytes or in other maternal tissue. Under control conditions (absence of MF hormone), both MP and NMP females produce female parthenogenetic offspring. Hence, any genes that are differentially expressed (DE) under control conditions represent constitutive gene expression differences between the two adult phenotypes and may be related to the potential of producing males (MP females) or absence thereof (NMP females). Exposure to MF leads to a change in offspring sex (male) in MP females, but not in NMP females (Innes and Dunbrack 1993; Innes 1997; Olmstead and Leblanc 2002; Olmstead and LeBlanc 2003; Tatarazako et al. 2003). It can be hypothesized that hormone exposure during the hormone-sensitive period may lead to a larger number and stronger gene-expression changes in MP than in NMP females. Specifically, if NMP is controlled by a simple loss of hormone-sensitivity mutation, the genes of the male sex-determination cascade that intervene upstream of this mutation can be hypothesized to show differential expression upon hormone exposure in both phenotypes, but downstream genes only in MP. Previous work has shown that expression of *DapmaDsx1*, one of the *doublesex* (*dsx*) genes of *D. magna*, is sufficient to induce male production (Kato et al. 2011). It is therefore likely that the NMP mutation interferes in the male sex-determination cascade upstream of *DapmaDsx1*. In addition, because MF is a hormone involved in various biological pathways (Homola and Chang 1997), additional genes, independent of the sex-determination pathway, may show a reaction to MF in both phenotypes.

Our experimental design comprised four distinct MP clones and four distinct NMP clones (they are referred to as “clones” because each of the eight clones is a genetically distinct isolate from a natural population, cultivated exclusively by parthenogenetic reproduction in the laboratory). The main analysis of differential gene expression was carried out for four contrasts: The first two tested differences between MP and NMP under the two experimental conditions (control: contrast “MP vs. NMP” and hormone exposure:
"MP-MF vs. NMP-MF"), the two others tested for changes in gene expression upon hormone exposure in the two phenotypes (MP: "MP vs. MP-MF" and NMP: "NMP vs. NMP-MF"). Statistical tests for the first two contrasts were carried out by comparing four MP clones with four NMP clones. In these tests, any among-clone variation in gene expression (Huylmans et al. 2016; Orsini et al. 2016; Tams et al. unpublished data) will increase the error variance and hence affect statistical power but not the likelihood of false positives. Tests for the other two contrasts were done in a pairwise design, thus taking into account clone identity, which was possible because each clone was tested in both conditions. In the presence of significant among-clone variation in gene expression, these tests should have improved statistical power. For some analyses, a fifth contrast was added: for a hierarchical clustering analysis (see below), as well as to assess whether DE genes showed sex-biased expression, we used gene expression data obtained from males of the same four MP clones used in the current study (data from Molinier et al. [2018]); males were reared at the same time as the present experiment, so as to reduce batch effects. For all four main contrasts, we investigated how many genes were DE, the direction of their expression, and assessed the genomic location of these genes with respect to the NMP chromosome region. Finally, we also investigated the identities of DE genes, specifically with respect to their possible involvement in sex-determination pathways.

## Results

### Data Generated and Basic Features of Data Analysis

RNA sequencing of the 16 libraries (eight clones, two culture conditions each) resulted in a total of 1.1 billion raw, Illumina paired-end reads. An average of 99.1% of raw reads passed quality control. After end-trimming, an average of 87.1% aligned to the reference genome, resulting in an average of 54 million aligned reads per library (see supplementary table S1, Supplementary Material online, for the percentages of reads retained at each step and sample). These reads were used to generate raw read count data for each of the 26,646 genes of the current *D. magna* genome annotation, which were then analyzed using state of the art software (mainly DESeq2, Love et al. 2014). Genes that were DE (multiple comparison-adjusted *P* value [*P*-adj] <0.05, |log₂ FC| > 1) in at least one of the four contrasts are listed in supplementary table S2, Supplementary Material online, along with *P*-adj and fold-change (FC) values for each contrast.

### Differential Gene Expression between MP and NMP Females

We found 126 DE genes between MP and NMP females under control conditions, of which 80 had a 2-fold-change or higher (table 1 and fig. 1). Of these 80 genes, 32.5% were overexpressed in MP females, and 67.5% in NMP females (table 2). The heatmap clearly distinguishes gene expression patterns between MP and NMP females (fig. 1). Without correction for multiple comparisons, 1,542 genes (out of 20,352 for which a *P* value was obtained) had *P* < 0.05, which is clearly higher

than the 1,018 genes (i.e., 5% of 20,352 genes) expected to be found by chance if all differences between the two phenotypes were purely random. The discrepancy between these figures (1,542 - 1,018 = 524 is substantially higher than 126) suggests that there may exist a class of DE genes that were not detected, likely due to too low statistical power.

When comparing MP females with NMP females when both were exposed to MF (MP-MF vs. NMP-MF), 265 genes were significantly DE (with multiple test correction, as in all subsequent tests), of which 163 had a 2-fold-change or more. Of these 163 genes, 48.5% were overexpressed in MP-MF, and 51.5% were overexpressed in NMP-MF (tables 1 and 2 and fig. 1). However, among the strongly DE genes (|log₂ FC| > 2, supplementary table S3, Supplementary Material online), there was a clear bias toward a higher number of genes being overexpressed in NMP compared with MP, both in control conditions (78.6%) and under hormone exposure (88.9%).

### Changes in Gene Expression upon Hormone Exposure

In MP females, 686 genes significantly changed their expression upon hormone exposure (table 1), of which 139 had a 2-fold-change or more (table 1). Of these 139 genes, 57.6% increased their expression upon exposure, and 42.4% reduced expression (table 2 and fig. 1). In NMP females, 1,563 genes significantly changed their expression upon hormone exposure, of which 566 had a 2-fold or more (table 1). Of these 566 genes, 99% increased their expression upon hormone exposure, and only 1% decreased expression (table 2 and fig. 1). Among the strongly DE genes (|log₂ FC| > 2), the bias toward a higher number of genes being upregulated rather than downregulated upon hormone exposure was even clearer as all of these genes were upregulated upon hormone exposure (both phenotypes), with the exception of a single gene which was downregulated upon hormone exposure in NMP (supplementary table S3, Supplementary Material online).

In the analysis of DE genes upon hormone exposure in MP (i.e., contrast MP vs. MP-MF), the clustering algorithm in DESeq2 clustered one of the control samples inside the hormone-exposed samples (fig. 1). However, as can be seen in the figure, many genes showed consistent differential expression between control conditions and hormone exposure in all four MP clones.

### Hierarchical Clustering Analysis

To globally identify the main drivers of gene expression differences, we carried out a hierarchical clustering analysis using the software WGCNA ("Weighted gene coexpression analysis") based on filtered and normalized read counts of 18,252 genes for the 20 samples (including the four males from Molinier et al. [2018]). We found that sex was by far the most important driver of variation in expression patterns (fig. 2). However, clone identity also seemed to have an important effect and apparently affected the general patterns of gene expression more strongly than the difference between control conditions and hormone exposure: In females, the two replicates of each clone (one exposed to MF, one not) always clustered together (fig. 2). In addition, MP clone 4 had the most distant expression pattern compared with the three
<table><caption>**Table 1.** Numbers of Significantly (*P*-adj < 0.05) DE Genes and Their Degree of Expression Bias in the Four Contrasts.</caption><thead><tr><th>Contrast</th><th>All</th><th>&lt;2-Fold</th><th>2- to 5-Fold</th><th>5- to 10-Fold</th><th>&gt;10-Fold</th><th>&gt;2-Fold</th></tr></thead><tbody><tr><td>MP vs. NMP</td><td>126</td><td>46</td><td>61</td><td>11</td><td>8</td><td>80</td></tr><tr><td>MP-MF vs. NMP-MF</td><td>265</td><td>102</td><td>143</td><td>10</td><td>10</td><td>163</td></tr><tr><td>MP vs. MP-MF</td><td>686</td><td>547</td><td>135</td><td>4</td><td>0</td><td>139</td></tr><tr><td>NMP vs. NMP-MF</td><td>1,563</td><td>997</td><td>552</td><td>13</td><td>1</td><td>566</td></tr></tbody></table>

NOTE.—Fold-changes are given without considering the direction of the bias, for each gene comparing the more strongly expressed condition to the more weakly expressed one (i.e., all fold-changes are >1). The direction of bias is listed in table 2.

<figure><img src="image_2.png" /><figcaption>**FIG. 1.** Heatmaps of differential gene expression in each of the four contrasts: (A) MP vs. NMP, (B) MP-MF vs. NMP-MF, (C) MP vs. MP-MF, and (D) NMP vs. NMP-MF. Each row represents a gene, each column a biological replicate. Relative expression intensities among biological replicates are shown with a color code, varying from dark blue (strongly underexpressed) to dark red (strongly overexpressed), as shown on the right-hand side of the heatmaps. The dendrograms on the top show the results of the hierarchical clustering as implemented in DESeq2, depicting gene expression differences among replicates (based on DE genes only).</figcaption></figure>
other MP clones in males as well as in females. Finally, among females, the four NMP clones clustered together, though the separation of this node is not very distinct. The presence of a relatively strong clone effect on overall gene expression patterns suggests that, in the main analysis of differential gene expression, the two within-clone contrasts (MP vs. MP-MF and NMP vs. NMP-MF) had higher statistical power to detect DE genes than the two between-clone contrasts (MP vs. NMP and MP-MF vs. NMP-MF).

## DE Genes in the NMP Chromosome Region

The physical position of 16,111 (out of 26,646) genes could be identified on the reference genetic map of *D. magna* (Dukić et al. 2016). The NMP chromosome region was enriched for genes that were DE between MP and NMP under control conditions: 15.6% of these genes (>2-fold-change) occurred in this region compared with only 2.8% expected by chance (P = 0.0003; table 3). The NMP chromosome region was also slightly enriched for genes that were DE between MP-MF and NMP-MF (P = 0.0041; table 3). In contrast, genes that changed expression upon hormone exposure (MP vs. MP-MF and NMP vs. NMP-MF) were slightly underrepresented in the NMP chromosome region, compared with chance expectations (table 3).

<table><caption>**Table 2.** Direction of Expression Bias of DE Genes and Test for Deviation from 50:50 Ratios for Each of the Four Contrasts.</caption><thead><tr><th>Contrast</th><th>Number of Genes Overexpressed in (1)</th><th>Number of Genes Overexpressed in (2)</th><th>Test for Deviation from 50:50, P Value</th></tr></thead><tbody><tr><td>MP (1) vs. NMP (2)</td><td>26</td><td>54</td><td>0.0023</td></tr><tr><td>MP-MF (1) vs. NMP-MF (2)</td><td>79</td><td>84</td><td>0.75</td></tr><tr><td>MP (1) vs. MP-MF (2)</td><td>59</td><td>80</td><td>0.09</td></tr><tr><td>NMP (1) vs. NMP-MF (2)</td><td>6</td><td>560</td><td>&lt;0.0001</td></tr></tbody></table>

## Genes that Are DE in More than One Contrast

Of the 80 DE genes between MP and NMP under control conditions (MP vs. NMP), 38 were also DE between MP and NMP when exposed to the MF (MP-MF vs. NMP-MF, fig. 3), which is significantly more than expected by chance, (P < 0.0001). All these 38 genes were biased in the same direction: 30 genes were overexpressed in NMP in both conditions (control and hormone exposure), and 8 genes were overexpressed in MP in both conditions (fig. 3). Of the genes that changed expression between control conditions and hormone exposure, 34 did so in MP as well as in NMP. The large majority of these expression changes showed the same directionality in the two phenotypes: 28 genes were overexpressed under hormone exposure and 2 genes were overexpressed under control conditions in both phenotypes. The only exception are four genes that were upregulated upon hormone exposure in NMP but downregulated upon hormone exposure in MP (fig. 3). A Venn-diagram showing the number of genes that were DE in any combination of all four contrasts is given in supplementary figure S1, Supplementary Material online. For the interpretation of this diagram, it is important to remember that the statistical analysis controls for the overall false discovery rate, but not for statistical power. Hence, the fact that a gene was not found to be DE in a specific contrast should be interpreted with caution. For all genes that were DE (P-adj < 0.05, |log<sub>2</sub> FC| > 1) in at least one contrast, the relative expression levels and significance tests for all four contrasts (even nonsignificant ones) are listed in supplementary table S2, Supplementary Material online. An analysis of these data reveals a strong positive correlation in the degree and direction of expression difference between MP and NMP when the expression difference (log<sub>2</sub> FC) under control conditions was plotted against the expression difference under hormone exposure (r = 0.78, N = 205, P < 0.005, fig. 4). In other words, genes that were significantly DE between MP and NMP in one of the two conditions (control condition or hormone exposure) were likely also DE (with an expression bias toward the same phenotype) in the other

<figure><img src="image_3.png" /><figcaption>**FIG. 2.** Dendrogram showing the results of the hierarchical clustering analysis of global gene expression patterns (all genes) and including, for comparison, also a male of each of the four MP clones.</figcaption></figure>
condition, especially for genes with a strong expression bias. A similar correlation also exists for genes that were significantly DE between control and hormone exposure: Their expression change (log₂ FC between control and hormone exposure) was positively correlated between the two phenotypes (r = 0.35, N = 671, P < 0.05, fig. 4). In other words, that a substantial proportion of genes that changed expression upon hormone exposure in a certain direction in one

phenotype did so also in the other one. Yet, the correlation was weaker compared with DE genes between MP and NMP. Indeed, closer inspection of figure 4 reveals the existence of a substantial number of genes that showed a reaction to hormone exposure that was specific to one of the two phenotypes. These include genes whose expression responded to hormone exposure in one phenotype but not in the other, as well as the four above-mentioned genes whose expression changed in opposite direction upon hormone exposure in MP compared with NMP.

<table><caption>**Table 3.** Number and Proportions of the DE Genes Inside and Outside the NMP Region.</caption><thead><tr><th>Contrast</th><th>Number (proportion) of DE Genes in the NMP Region</th><th>Number of DE Genes Outside the NMP Region</th><th>P Value (Fisher's exact test)</th></tr></thead><tbody><tr><td>MP vs. NMP</td><td>7 (15.6%)</td><td>38</td><td>0.0003</td></tr><tr><td>MP-MF vs. NMP-MF</td><td>9 (8.18%)</td><td>101</td><td>0.0041</td></tr><tr><td>MP vs. MP-MF</td><td>0</td><td>89</td><td>0.19</td></tr><tr><td>NMP vs. NMP-MF</td><td>3 (0.7%)</td><td>424</td><td>0.0029</td></tr></tbody></table>

NOTE.—The expected proportion of DE genes in the NMP region is 2.8% (458 of 16,111 mapped genes). Significant P values indicate significant overrepresentation (>2.8%) or significant underrepresentation (<2.8%) of DE genes in the NMP region.

## Sex-Biased Expression of DE Genes in the Four Contrasts

We also compared the genes being DE in the four contrasts with a list of sex-biased genes (Molinier et al. 2018, fig. 3). Interestingly, among the genes that were DE between control conditions and hormone exposure, we found a higher proportion of sex-biased genes in MP (i.e., in the contrast MP vs. MP-MF) than in NMP (contrast NMP vs. NMP-MF; 59%, N = 139 compared with 33%, N = 566, Fisher's exact test, P < 0.0001). Furthermore, taking as a reference the overall proportion of genes with male-biased expression (57.7%) among all sex-biased genes (N = 8,384, Molinier et al. 2018), we found that, in MP, hormone exposure led to a significant

A

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

C

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

E

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

B

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

D

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

F

<figure><img src="image_9.png" /><figcaption>**FIG. 3.** Number of genes found to be DE in multiple contrasts. (A) Number of genes being DE between control and hormone exposure in NMP (left circle), MP (right circle), or both (overlap). (B) Number of genes being DE between MP and NMP under hormone exposure (left circle), control conditions (right circle), or both (overlap). (C–F) Number of genes found to be DE in each contrast compared with number of genes with sex-biased expression (from Molinier et al. [2018]). The tables next to the Venn-diagrams show the direction of the expression bias for genes that were found to be DE in only one contrast (margins) and genes being DE in both contrasts (center).</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>FIG. 4. Degree of differential expression for genes being DE (*P*-adj < 0.05, ||log₂FC| > 1) in at least one contrast. (A) Genes changing their expression upon hormone exposure. (B) Genes being DE between MP and NMP. Positive values relate to genes being downregulated under hormone exposure (A) or overexpressed in NMP relative to MP (B).</figcaption></figure>

upregulation of male-biased genes (Fisher's exact test, *P* = 0.0014) and to a significant downregulation of female-biased genes (*P* = 0.0005). In contrast, in NMP, hormone exposure led to a significant upregulation of female-biased genes (*P* = 0.016).

Of the genes that were DE between MP and NMP, about 50% were also sex-biased, both in control conditions and under hormone exposure. Under control conditions, there was no deviation from expectations in the proportion of female-biased versus male-biased genes, neither among genes overexpressed in MP nor among genes overexpressed in NMP (*P* = 0.29 and *P* = 0.16, respectively, fig. 3). However, under hormone exposure (MP-MF vs. NMP-MF), the genes that were overexpressed in MP showed a strong tendency to be male biased (54 male-biased genes compared with only three female-biased ones, *P* < 0.0001, fig. 3). There was no such tendency among genes that were overexpressed in NMP (*P* = 0.8345).

## Gene Ontology (GO) Enrichment

The GO-enrichment analysis was carried out on all DE genes (*P*-adj < 0.05, no restriction on ||log₂FC||), in order to have a sufficient number of genes for the analysis of each of the four contrasts. The results are reported in supplementary figure S2, Supplementary Material online. Among the genes that are DE between MP and NMP under control conditions, "catabolism activities of sugar" were overrepresented in the category "biological processes" and "peptidase activity" in the category "molecular functions". Under hormone exposure, several other biological processes, mainly related to DNA replication, were overrepresented among the DE genes between MP and NMP. Several molecular functions linked to transcription and regulation (e.g., "nucleic acid binding transcription factor activity," "DNA binding," and "sequence-specific DNA binding") were overrepresented (supplementary fig. S2, Supplementary Material online). The GO-enrichment

analyses of genes being DE between control conditions and hormone exposure showed an overrepresentation of many different GO-terms in both MP and NMP, with some common terms linked to carbohydrates and lipids, such as "lipid transport" (supplementary fig. S2, Supplementary Material online).

## Functional Annotation of the Genes with the Strongest Expression Differences

We first concentrated on genes that showed particularly strong expression differences (||log₂FC| > 2, i.e., being DE by more than a 4-fold-change). In total, 90 such genes were found in at least one of the contrasts (supplementary table S3, Supplementary Material online). From the General Feature Format (GFF) file or the BLAST2GO analysis, we were able to identify a functional annotation for a bit <50% of these genes (supplementary table S3, Supplementary Material online). Several of these genes involved in carbohydrate-related processes, which have already earlier been identified as important components of MF pathways and/or sex-determination pathways (Toyota et al. 2017). Other notable gene families included hemoglobin and serine protease, which also had been identified as components of the MF and/or sex-determination pathway (Toyota et al. 2014, 2015, 2017; Abe et al. 2015). The genes with the highest expression differences in reaction to hormone exposure were a cell wall-associated hydrolase (in MP) and a cytosolic sulfotransferase 1B family member (in MP). In addition, "tectonin beta-propeller repeat-containing 2," which may be involved in autophagy and cell cycle regulation (Alexander et al. unpublished data), was strongly (almost 100 times) upregulated in NMP compared with MP under both control conditions and hormone exposure.

Included in the list of genes with strong expression differences are also the four genes that showed an opposite reaction to hormone exposure in MP compared with NMP, being
upregulated upon hormone exposure in NMP but downregulated upon hormone exposure in MP (fig. 3). Although the individual FCs are lower than for the above genes, the changes were in opposite direction in the two phenotypes, which was the rationale for including them. Three of them had a functional annotation (supplementary table S3, Supplementary Material online), yet with unclear meaning with respect to the differential reaction of the two phenotypes to the hormone. Nonetheless, all four genes are interesting candidates for being involved, as key components, in the phenotype-specific effects of hormone exposure.

## Functional Annotation of DE Genes in the NMP Chromosome Region

Across all four contrasts, 14 DE genes ($P$-adj < 0.05 and $|log_2 FC| > 1$) were found in the NMP chromosome region, but only five with an annotated function (supplementary table S4, Supplementary Material online). Most notably, one of the genes that was DE between MP and NMP under control conditions, “ecdysone 20-monooxygenase isoform X1” (Weirich et al. 1984), is a key enzyme in the formation of the moulting hormone which belongs to the cytochrome P450 family known to be involved in sex determination in other organisms (Doctor 1985; Verma 1996).

## Expression Patterns of Genes with Known Functions

We investigated whether genes that had earlier been identified as potential candidate genes involved in the NMP phenotype or in sex determination/sex differentiation are DE in any of the four contrasts. First, we searched for the most probable *D. magna* homolog of gene 8960, the gene primarily responsible for the NMP phenotype in *Daphnia pulex* (Ye Z, Molinier C, Zhao C, Haag CR, Lynch M, in review). We obtained just a single hit, in the *D. magna* assembly, as well as in the RNA-sequencing data by Orsini et al. (2016): Gene APZ42_021088 (no functional annotation) shows 51.2% amino acid similarity to gene 8960. It is located between positions 633974 and 635100 on scaffold 1036 of the 2.4 *D. magna* assembly, on linkage group 10 of the reference genetic map (Dukić et al. 2016). Because it is outside the NMP chromosome region on LG 3, it is unlikely that gene APZ42_021088 carries the causal mutation responsible for NMP in *D. magna*. Furthermore, in contrast to *D. pulex*, where gene 8960 is upregulated upon hormone exposure in MP but not in NMP phenotypes (Ye Z, Molinier C, Zhao C, Haag CR, Lynch M, in review), gene APZ42_021088 was not DE in any of the contrasts except for NMP vs. NMP-MF, where it was slightly upregulated under hormone exposure (supplementary table S5, Supplementary Material online). Second, we compiled a list based on published literature of 34 candidate genes potentially involved in sex determination in *D. magna* and checked if they were DE in any of the contrasts. Three genes were significantly DE (supplementary table S5, Supplementary Material online): *Doublesex1* was overexpressed in NMP compared with MP (both experimental conditions) but did not show an expression change upon hormone exposure. Furthermore, *FTZ-F1* was downregulated upon hormone exposure in MP but did not show a change in

expression upon hormone exposure in NMP, and *Vrille* was upregulated upon hormone exposure in MP and NMP (though in NMP the FC was <2, supplementary table S5, Supplementary Material online).

## Discussion

### Constitutive Differences in Gene Expression between MP and NMP Phenotypes

We found a bit more than 100 genes to be constitutively DE between MP and NMP even in the absence of hormone. This is far less than, for instance, the number of DE genes between sexes (Molinier et al. 2018). However, under control conditions, phenotypic differences between NMP and MP individuals are invisible: both phenotypes are females producing female parthenogenetic offspring. Thus, the gene expression differences between MP and NMP females are likely linked to the constitutive ability (or the absence thereof) of these females to produce males in case of an environmental change (or under artificial hormone exposure). Furthermore, the expression differences between phenotypes were strongly correlated between control conditions and hormone exposure, suggesting that the constitutive expression differences are largely condition-independent.

The constitutively DE genes were particularly concentrated in the NMP chromosome region. This large, nonrecombining region likely contains several genetic variants involved in the NMP phenotype (Reisser et al. 2017). Among them is the primary variant causing a loss of male-producing function, perhaps a loss-of-function mutation in a gene essential for making males, or in the environment-dependent switch, from female to male offspring production. However, the region likely contains several additional mutations that are functionally involved in the NMP phenotype (Reisser et al. 2017). These functions may be related to secondary fine-tuning and integration of the effects of the primary mutation. Notwithstanding these points, it is possible that some of the mutations are located in regulatory regions, leading to differential gene expression, and thus potentially explaining the enrichment of the region for DE genes.

### Gene Expression Differences in Response to the Hormone in MP and NMP Phenotypes

Hormone exposure induced changes in gene expression in a large number (several 100s to over 1000s) of genes in MP and NMP phenotypes, respectively. A subset of these genes showed parallel changes in the two phenotypes (i.e., changes in the same direction upon hormone exposure). These may be genes involved in hormone pathways that are either unrelated to the induction of the male production cascade, or in the part of the cascade that is common to both phenotypes. Yet, a substantial fraction of genes reacted differently to hormone exposure in the two phenotypes: They either showed changes in opposite direction upon hormone exposure or showed strong expression changes in one phenotype while reacting only weakly or not at all in the other.

Genes showing MP-specific expression changes upon hormone exposure might be involved in the cascades that lead to
male production exclusively in MP. In line with this idea, hormone exposure led to a "masculinization" of gene expression in MP females, that is, to an upregulation of male-biased genes and to a downregulation of female-biased genes. Note, however, that the developing oocytes, whose sex was determined during hormone exposure, were still inside their mothers' ovaries at the time of sampling (Hiruta et al. 2010). At this stage, they make up only a tiny fraction of the whole tissue sampled, and their transcriptional activity is unclear. Most of the hormone-induced gene-expression changes likely occurred in the mothers. This suggests that male production in MP via ESD is under maternal control, initiated by hormone-dependent, maternal gene-expression changes taking place much before the termination of oogenesis.

Inferring the role of the genes with NMP-specific expression changes upon hormone exposure is less straightforward, due to a lack of a visible phenotype. It is possible that some of these changes are due to unspecific reactions to the hormone, unrelated to the phenotype-specific cascades that lead to continued female production in NMP. Yet, unspecific changes are difficult to reconcile with the larger number of genes that changed expression upon hormone exposure in NMP than in MP and with the observation of a slight but significant "feminization" of gene expression (i.e., preferential upregulation of genes with female-biased expression) in NMP. It is thus possible that, despite the constitutive expression differences between MP and NMP, extensive gene regulation changes are required for homeostasis of the female-producing pathway under hormone exposure in NMP.

## Functional Analysis of Genes Involved in the NMP Phenotype Determination

The functional annotation of DE genes revealed a large proportion of genes with unknown functions. Nonetheless, for each contrast, a list of genes with known function could be identified, and these lists contain several candidates of key genes involved in constitutive differences between MP and NMP, as well as in their differential reaction to hormone exposure. Among the most promising candidates are genes located in the NMP chromosomal region, genes with large FCs, as well as genes that are DE in several contrasts, such as those with opposite responses to hormone exposure in the two phenotypes. However, even the genes with annotated functions are involved in a multitude of biological processes, suggesting that the molecular construction of the phenotypic differences between MP and NMP, as well as their differential reaction to hormone, is multifaceted.

Some of the genes identified in this paper as reacting to hormone exposure had already been identified before (Eads et al. 2008; Hannas et al. 2011; Toyota et al. 2014, 2015, 2017; Abe et al. 2015), including some genes known to be involved in sex determination. In particular, the expression of *Doublesex1* (*Dsx1*) is known to be sensitive to hormone treatment and sufficient to trigger male production in MP (*Xu et al. 2014; LeBlanc and Medlock 2015*). Indeed, *Dsx* plays a key role in sex differentiation in many organisms (*Burtis and Baker 1989; Verhulst et al. 2010*). In *D. magna*, *Dsx1* is overexpressed in early male development compared with females

(Kato et al. 2011; Nong et al. 2017), but the difference appears only post-ovulation (Nong et al. 2017). The likely reason for why we did not find *Dsx1* to be DE between control conditions and hormone exposure in MP is thus that the hormone-triggered change of *Dsx1* expression reported in the literature occurs downstream of the expression changes observed in our study. We did, however, observe a clear constitutive expression difference of *Dsx1* between MP and NMP females in both culturing conditions, with *Dsx1* being overexpressed in NMP. *Dsx1* may thus play a role in the constitutive differences between phenotypes.

Gene APZ42_021088, the putative homolog of *gene 8960* in *D. pulex*, clearly does not have the same central function in determining the NMP phenotype as in *D. pulex* (Ye Z, Molinier C, Zhao C, Haag CR, Lynch M, in review). Neither is it located in the NMP chromosome region, nor was it found to be DE in MP upon hormone exposure. However, it was upregulated upon hormone exposure in NMP, which suggests that it may still be involved in the same pathway as in *D. pulex*, though in a different role. These observations, together with the fact that no other putative homolog of *gene 8960* was found in *D. magna*, suggest parallel evolution of the NMP phenotype in *D. magna* and *D. pulex*, which is probably not so surprising given their divergence time of ~150 My (Kotov and Taylor 2011).

## General and Evolutionary Implications of the Results

Exposure to MF leads to a change in phenotype in MP but not in NMP females. Our initial hypothesis was that hormone exposure should lead to more pronounced changes in gene expression in MP than in NMP females. Furthermore, if NMP was entirely controlled by a simple loss of hormone sensitivity, expression differences between the two phenotypes may occur only in presence of the hormone (due to the hormone-sensitive genes being differentially regulated in MP). Our results do not conform to these expectations: A substantial number of genes showed constitutive expression differences between MP and NMP under control conditions, and hormone exposure led to expression changes in a larger number of genes in NMP than in MP. Taken together, our findings suggest that the homeostasis of female production is regulated by a rather complex mechanism, involving differential expression of many genes (constitutive and hormone induced). It seems likely that some of these changes are due to secondary modifications, which occurred after the establishment of the initial female-determining mutation. Overall, the complex nature of the gene expression patterns underlying the maintenance of female production suggests that the evolution of the NMP phenotype is not a very recent event. Concomitantly, this also suggests that the sex-chromosome-like region (the NMP chromosome region) is not very young and likely has incurred secondary changes, including possibly sex-antagonistic mutations and additional recombination suppression. In other words, the NMP chromosome region appears to be at an advanced stage of evolution of an incipient W chromosome. More generally, our study illustrates that the evolution of a genetically determined sex from ESD
may be complex, involving modifications of multiple genes and pathways.

## Materials and Methods

### Daphnia Clones Used in the Study

All *D. magna* clones used in this study originated from a single population in Moscow (55.763514°N, 37.581667°E). The clones have been used in another study (Reisser et al. 2017), where their phenotype was verified by hormone tests, and genetic markers revealed that each of the clones was genetically distinct.

### Preparation and Treatments for RNA Sequencing

We carried out RNA sequencing on adult NMP females and adult MP females, kept under control conditions (i.e., standard culturing conditions under which the females typically produce female parthenogenetic offspring) or exposed to hormone prior to sampling. Each phenotypic class and experimental condition was replicated four times by using four NMP clones and four MP clones ("biological replicates"). One library was prepared per biological replicate, resulting in a total of 16 libraries. Furthermore, each library was based on eight technical replicates, that is, eight replicate individuals of the same genotype, phenotype and treatment. Hence, a total of 128 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. The eight individuals per biological replicate were pooled just prior to RNA extraction.

### Experimentation

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 °C under a 16:8 h light–dark photoperiod. Each technical replicate was reared under these standard conditions during two pre-experimental clonal generations to remove maternal effects due to different culturing conditions (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 generation. Third-generation offspring were used for RNA sequencing. The experimental procedure ensured that these individuals were derived from germinal cells that had started their differentiation under standard culturing conditions. Throughout the experiment culture medium was exchanged daily. Just before the moult during which third-generation females released their first clutch into the water column, all individuals were transferred individually to a 1.5-ml well on a culture plate, where they were kept for about three days. Standard medium (controls) or medium containing 400 nM of MF (MF treatment, MF was obtained from Echelon Biosciences, catalog number S-0153) was exchanged

daily, and the juveniles released during the moult were removed (the new clutch that was deposited into the brood pouch shortly after this moult was removed after sampling). All individuals were sampled 60 h past moult, which corresponds to the moment of highest sensibility of the maturing oocytes in the ovaries to MF for male production (Olmstead and Leblanc 2002). 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 induction of starvation-dependent gene expression. The efficacy of the hormone batch was successfully tested on the second generation females: Using MF, we produced brothers of the MP females used in the experiment (Molinier et al. 2018).

### Sampling

To remove as much culture medium as possible, the individuals were blotted with absorbing paper (previously sterilized with ultraviolet for 30 min), and then transferred to a 1.5-ml tube that was directly immersed in liquid nitrogen. Directly after flash-freezing, three volumes of RNA later ICE solution were added to preserve RNA, and samples were placed at -80 °C. Treatment with RNA later ICE is advantageous when samples have to be dissected prior to RNA extraction, as it prevents RNA degradation during thawing and dissection. Here, prior to RNA extraction, all eggs were removed from the brood pouch of females to avoid noise from developmental genes not induced under treatment conditions (developing oocytes, whose sex was induced under treatment conditions, were still in the ovaries at the time of sampling; these eggs would have been deposited in the brood pouch only during the next moult). The technical replicates (see above) were subsequently pooled.

### RNA Extraction, Library Preparation, and RNA Sequencing

Total RNA extraction and purification of the 16 samples was carried out following the protocol of the *Daphnia* genomic consortium (DGC; DGC, Indiana University, October 11, 2007), using Trizol Reagents and the Qiagen RNEasy Mini Kit. The extracted and purified RNA samples were then put at -80 °C and shipped on dry ice to the BSSE Genomic Sequencing Facility, University of Basel, Switzerland. Each of the 16 samples was labeled using the TruSeq preparation kit. All libraries were sequenced using a single flow cell on an Illumina NextSeq 500 sequencer with 76 cycles in paired-end (strand information kept).

### Quality Control

The software FastQC v.0.10.1 (http://www.bioinformatics.babraham.ac.uk/projects/fastqc, last accessed April 21, 2018) was used to analyze read quality. The paired-end sequences were subjected to adapter trimming and quality filtering using trimmomatic v.0.36 (Bolger et al. 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 as long as their average quality scores were below 15.

## Mapping and Counting

Filtered reads were mapped to the *D. magna* genome assembly (v2.4; GenBank assembly accession: GCA_001632505.1, including a genome annotation GFF-file with 26,646 genes) using the RNA-Seq aligner STAR (Dobin et al. 2013) with default settings. The raw counts (number of mapped reads per transcript per sample) were obtained with the software program featureCounts (Liao et al. 2014). Counts were summarized at the gene level using the GFF annotation file.

## Differential Gene Expression

Differential expression analysis was carried out with DESeq2 (version 1.10.1) implemented in R (Love et al. 2014). In the following text, "NMP" or "MP" refer to individuals that were reared under normal culturing conditions whereas "NMP-MF" and "MP-MF" refer to those treated with MF. Four pairwise comparisons ("contrasts") were carried out. Because the four MP clones were different from the four NMP clones, a one-factor analysis was used for the comparisons MP vs. NMP and MP-MF vs. NMP-MF. However, as the same clones were used in both treatments, the contrasts MP vs. MP-MF and NMP vs. NMP-MF were analyzed using a two-factor design, as implemented in DESeq2, thereby taking into account clone identity. This assures that tests are carried out according to the pairwise design. All *P* values were adjusted for multiple testing with the Benjamini-Hochberg method as implemented in DESeq2. Genes were considered DE if they had an adjusted *P* value <0.05 (false discovery rate = 5%). The magnitude of differential expression was determined by the FC differences between the phenotypes or treatments.

## Hierarchical Clustering Analysis

After a classical analysis of DE genes, we used the package WGCNA version 1.50 (Langfelder and Horvath 2008), implemented in R. WGCNA is a hierarchical clustering-based method that creates networks of genes whose expression is similar among each other in the sample set, and then correlates these networks (called modules) with particular traits and conditions of interest.

Raw read counts were first filtered by removing genes with counts fewer than three in more than 75% of samples and then normalized using the variance stabilizing normalization available in DESeq2 (see below). These filtered and normalized counts led to 18,252 genes to be used as input in the WGCNA package. The analysis was performed on an extended data set with the 20 biological replicates, including the four males from a previous study (Molinier et al. 2018) because they have a strong pattern of differential expression.

We used this method as a broad-scale method to explore the main factors driving global variance in gene expression patterns. Thus, WGCNA provides an opportunity to assess all potential drivers (phenotype, treatment, and clone) simultaneously, rather than using the contrasts of the DESeq2 analysis. The classical use of the method (identification of

modules of genes with similar pattern of expression changes) was not employed here, due to a lack of a high enough number of DE genes in most of our contrasts.

## Location of DE Genes on the Genetic Map

To investigate the location of DE genes (*P* < 0.05; |log₂ FC| > 1) with respect to the NMP chromosome region, we used a high-density genetic map of *D. magna* (Dukić et al. 2016). For the purpose of this study, we defined the NMP chromosome region as the region between cM-positions 69 and 95 cM on LG 3, that is, slightly larger than the one previously used by Reisser et al. (2017) because closely linked genes could contribute to sex chromosome evolution. Among the 16,111 mapped genes, 2.8% were within the NMP chromosome region. For each of the four contrasts, we then estimated the proportion of DE genes inside and outside the NMP chromosome region and tested for deviations from the expected proportions using two-tailed Fisher's exact tests.

## Gene Ontology Enrichment Analysis

For each of the four contrasts, we performed a gene ontology enrichment analysis to test for overrepresentation of GO-terms in biological processes and molecular functions among the DE genes (*P* < 0.05). This was done using the GOatools Python script (https://github.com/tanghaibao/goatools), which performs a Fisher's exact test of overrepresentation GO-terms of DE genes compared with non-DE genes. GO-terms with *P* < 0.01 were considered significantly enriched. We used the software REVIGO (Supek et al. 2011) to summarize enriched GO-terms (by a reduction of GO-term complexity and levels) and to visualize them.

## Annotation of DE Genes

In order to functionally annotate the DE genes found in each of the four contrasts, we extracted the annotation from the UniprotKB database, corresponding to the genome annotation (GFF file used). Additionally, and in order to potentially further complete this annotation, we performed a Blast2GO annotation (version 4.0.7, Conesa et al. 2005), using the NCBI nr database, allowing for 20 output alignments per query sequence with an *e*-value threshold of 0.001. The subsequent mapping and annotation steps implemented in BLAST2GO were run with default settings. Additionally, InterPro IDs from InterProScan were merged to the annotation for further accuracy. Using these annotations, we first concentrated on the genes that showed the strongest expression differences in any of the four contrasts (we used an arbitrary cut-off of |log₂ FC| > 2). Second, we inspected the function of DE genes in the NMP chromosome region. Third, we searched for genes potentially involved in sex determination or sex differentiation in the entire list of *D. magna* genes to assess whether any of these are DE. Based on the literature (Caudy et al. 1988; Murre et al. 1994; Zelzer et al. 1997; Heinrichs et al. 1998; Tokishita et al. 2006; Hasselmann et al. 2008; Kato et al. 2008, 2010, 2011; Verhulst et al. 2010; LeBlanc et al. 2013; Guo et al. 2015; Herpin and Schartl 2015; LeBlanc and Medlock 2015; Mohamad Ishak et al. 2016; Toyota et al. 2016, 2017; Mohamad Ishak et al. 2017; Nong et al. 2017; Ye Z, Molinier C, Zhao C, Haag CR,
Lynch M, in review), we compiled a list of candidate genes potentially involved in sex determination or sex differentiation with a specific emphasis on crustaceans and insects. If no direct information on *D. magna* was available, we tried to identify the most probable homolog of these candidates in *D. magna* using the best hit from a blast of the protein sequence using BLAST2GO.

## Supplementary Material

Supplementary data are available at *Molecular Biology and Evolution* online.

## Acknowledgments

We thank D. Frey for help with indoors cultures. We thank C. Beisel and I. Nissen, Department of Biosystems Science and Engineering, ETH Zurich for Illumina sequencing, and gratefully acknowledge support by M.-P. Dubois, the 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éen Environnement Biodiversité (CeMEB). We thank the University of Fribourg, the Montpellier Bioinformatic Biodiversity platform, and the Labex CeMEB for access to high-performance computing. We also thank M. Cordellier and an anonymous reviewer for their constructive and helpful comments. The sequence data of the *D. magna* genome V2.4, available through GenBank assembly accession GCA_001632505.1, were produced by the Center for Genomics and Bioinformatics, Indiana University, and distributed via wFleaBase in collaboration with the *Daphnia* Genomics Consortium (National Institutes of Health grant 5R24GM078274-02 "Daphnia Functional Genomics Resources"). This work was supported by the Swiss National Science Foundation (grant no. 31003A_138203 to C.R.H.), the Russian Foundation for Basic Research (grant 16-04-01579 to Y.G.), and the European Union (Marie Curie Career Integration grant PCIG13-GA-2013-618961, DamaNMP to C.R.H.). Raw Illumina sequence reads from this study have been deposited in the NCBI SRA database under the BioProject accession number PRJNA533017. The sequence of *D. pulex* gene 8960 is available through the *D. pulex* genome assembly PA42 v3.0, GenBank assembly accession GCA_900092285.1, with the annotation file S3 from Ye et al. (2017).

## References

* Abe R, Toyota K, Miyakawa H, Watanabe H, Oka T, Miyagawa S, Nishide H, Uchiyama I, Tollefsen KE, Iguchi T, et al. 2015. Diofenolan induces male offspring production through binding to the juvenile hormone receptor in *Daphnia magna*. *Aquat Toxicol.* 159:44–51.

* Alexander J, Stroebel T, Georgitsi M, Schuster M, Penz T, Bock C, Honigschnabl S, Reiner A, Fischer P, Paschou P, Kovacs GG, unpublished data, www.biorxiv.org/content/10.1101/157693v1, last accessed May 22, 2019.

* Ashman T. 2002. The role of herbivores in the evolution of separate sexes. *Ecology* 83(5):1175–1184.

* Bachtrog D, Hom E, Wong KM, Maside X, de Jong P. 2008. Genomic degradation of a young Y chromosome in *Drosophila miranda*. *Genome Biol.* 9(2):R30.

* Bachtrog D, Perrin N, Ming R, Valenzuela N, Mayrose I, Peichel CL, Hahn MW, Ashman T-L, Vamosi JC, Ross L, et al. 2014. Sex determination: why so many ways of doing it? *PLoS Biol.* 12(7):e1001899.

* Barrett S. 2010. Understanding plant reproductive diversity. *Philos Trans R Soc B* 365(1537):99–109.

* Bergero R, Charlesworth D. 2009. The evolution of restricted recombination in sex chromosomes. *Trends Ecol Evol.* 24(2):94–102.

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

* Bolger AM, Lohse M, Usadel B. 2014. Trimmomatic: a flexible trimmer for Illumina sequence data. *Bioinformatics* 30(15):2114–2120.

* Burtis KC, Baker BS. 1989. Drosophila *doublesex* gene controls somatic sexual differentiation by producing alternatively spliced mRNAs encoding related sex-specific polypeptides. *Cell* 56(6):997–1010.

* Cáceres CE. 1998. Interspecific variation in the abundance, production, and emergence of *Daphnia* diapausing eggs. *Ecology* 79(5):1699–1710.

* Caudy M, Vässin H, Brand M, Tuma R, Jah LY, Jan YN. 1988. *daughterless*, a *Drosophila* gene essential for both neurogenesis and sex determination, has sequence similarities to myc and the achaete-scutte complex. *Cell* 55(6):1061–1067.

* Charlesworth B. 1991. The evolution of sex chromosomes. *Science* 251(4997):1030–1033.

* Charlesworth B, Charlesworth D. 2000. The degeneration of Y chromosomes. *Philos Trans R Soc B Biol Sci.* 355(1403):1563–1572.

* Charlesworth D. 2006. Evolution of plant breeding systems. *Curr Biol.* 16(17):R726–R735.

* Conesa A, Gotz S, Garcia-Gomez JM, Terol J, Talon M, Robles M. 2005. Blast2GO: a universal tool for annotation, visualization and analysis in functional genomics research. *Bioinformatics* 21(18):3674–3676.

* Dobin A, Davis CA, Schlesinger F, Drenkow J, Zaleski C, Jha S, Batut P, Chaisson M, Gingeras TR. 2013. STAR: ultrafast universal RNA-seq aligner. *Bioinformatics* 29(1):15–21.

* Doctor J. 1985. The pupal cuticle of *Drosophila*: biphasic synthesis of pupal cuticle proteins in vivo and in vitro in response to 20-hydroxyecdysone. *J Cell Biol.* 101(1):189–200.

* Dorken ME, Barrett S. 2004. Sex determination and the evolution of dioecy from monoecy in *Sagittaria latifolia* (Alismataceae). *Proc R Soc B Biol Sci.* 271(1535):213.

* Dukić M, Berner D, Roesti M, Haag CR, Ebert D. 2016. A high-density genetic map reveals variation in recombination rate across the genome of *Daphnia magna*. *BMC Genet.* 17(1):137.

* Eads BD, Andrews J, Colbourne JK. 2008. Ecological genomics in *Daphnia*: stress responses and environmental sex determination. *Heredity* (Edinb). 100(2):184–190.

* Ebert D. 2005. Ecology, epidemiology and evolution of parasitism in *Daphnia*. Bethesda (MD): National Center for Biotechnology Information.

* Galimov Y, Walser B, Haag CR. 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.

* Gorbi G, Moroni F, Sei S, Rossi V. 2011. Anticipatory maternal effects in two different clones of *Daphnia magna* in response to food shortage. *J Limnol.* 70(2):222–230.

* Guo X, Xu S, Yan X, Zhou W, Dai X, Zou X, Wang C, Wang D, Zhao Y. 2015. Cloning and expression of a *Chk1* gene in *Daphnia pulex* during different modes of reproduction. *Genes Genom.* 37(9):797–807.

* Hannas BR, Wang YH, Thomson S, Kwon G, Li H, LeBlanc GA. 2011. Regulation and dysregulation of vitellogenin mRNA accumulation in daphnids (*Daphnia magna*). *Aquat Toxicol.* 101(2):351–357.

* Hasselmann M, Gempe T, Schiott M, Nunes-Silva CG, Otte M, Beye M. 2008. Evidence for the evolutionary nascence of a novel sex determination pathway in honeybees. *Nature* 454(7203):519–522.

* Hebert P. 1978. The population biology of *Daphnia* (Crustacea, Daphnidae). *Biol Rev.* 53(3):387–426.
Heinrichs V, Ryner LC, Baker BS. 1998. Regulation of sex-specific selection of fruitless 5' splice sites by transformer and transformer-2. *Mol Cell Biol.* 18(1):450-458.

Herpin A, Scharl M. 2015. Plasticity of gene-regulatory networks controlling sex determination of masters, slaves, usual suspects, newcomers, and usurpersors. *EMBO Rep.* 16(10):1260-1274.

Hiruta C, Nishida C, Tochinai S. 2010. Abortive meiosis in the oogenesis of parthenogenetic *Daphnia pulex*. *Chromosome Res.* 18(7):833-840.

Homola E, Chang ES. 1997. Methyl farnesoate: crustacean juvenile hormone in search of functions. *Comp Biochem Physiol B Biochem Mol Biol.* 117:347-356.

Huylmans AK, Ezquerra AL, Parsch J, Cordellier M. 2016. De novo transcriptome assembly and sex-biased gene expression in the cyclical parthenogenetic *Daphnia galeata*. *Genome Biol Evol.* 8(10):3120-3139.

Innes DJ. 1997. Sexual reproduction of *Daphnia pulex* in a temporary habitat. *Oecologia* 111(1):53-60.

Innes DJ, Dunbrack RL. 1993. Sex allocation variation in *Daphnia pulex*. *J Evolution Biol.* 6(4):559-575.

Janzen FJ, Paukstis GL. 1991. Environmental sex determination in reptiles: ecology, evolution, and experimental design. *Q Rev Biol.* 66(2):149-179.

Kato Y, Kobayashi K, Oda S, Colbourn JK, Tatarazako N, Watanabe H, Iguchi T. 2008. Molecular cloning and sexually dimorphic expression of DM-domain genes in *Daphnia magna*. *Genomics.* 91(1):94-101.

Kato Y, Kobayashi K, Oda S, Tatarazako N, Watanabe H, Iguchi T. 2010. Sequence divergence and expression of a transformer gene in the branchiopod crustacean, *Daphnia magna*. *Genomics* 95(3):160-165.

Kato Y, Kobayashi K, Watanabe H, Iguchi T. 2011. Environmental sex determination in the branchiopod crustacean *Daphnia magna*: deep conservation of a Doublesex gene in the sex-determining pathway. *PLoS Genet.* 7(3):e1001345.

Khil PP, Smirnova NA, Romanienko PJ, Camerini-Otero RD. 2004. The mouse X chromosome is enriched for sex-biased genes not subject to selection by meiotic sex chromosome inactivation. *Nat Genet.* 36(6):642-646.

Kleiven OT, Larsson P, Hobaek A, Hobaek A. 1992. Sexual reproduction in *Daphnia magna* requires three stimuli. *Oikos* 65(2):197.

Klüttgen B, Dülmer U, Engels M, Ratte HT. 1994. ADaM, an artificial freshwater for the culture of zooplankton. *Water Res.* 28(3):743-746.

Kotov AA, Taylor DJ. 2011. Mesozoic fossils (>145 Mya) suggest the antiquity of the subgenera of *Daphnia* and their coevolution with chaoborid predators. *BMC Evol Biol.* 11:129.

Langfelder P, Horvath S. 2008. WGCNA: an R package for weighted correlation network analysis. *BMC Bioinformatics* 9(1):559.

LeBlanc GA, Medlock EK. 2015. Males on demand: the environmental-neuro-endocrine control of male sex determination in daphnids. *FEBS J.* 282(21):4080-4093.

LeBlanc GA, Wang YH, Holmes CN, Kwon G, Medlock EK. 2013. A transgenerational endocrine signaling pathway in Crustacea. *PLoS One* 8(4):e61715.

Liao Y, Smyth GK, Shi W. 2014. FeatureCounts: an efficient general purpose program for assigning sequence reads to genomic features. *Bioinformatics* 30(7):923-930.

Lindholm A, Breden F. 2002. Sex chromosomes and sexual selection in poeciliid fishes. *Am Nat.* 160(S6):S214.

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

Lynch M, Seyfert A, Eads B, Williams E. 2008. Localization of the genetic determinants of meiosis suppression in *Daphnia pulex*. *Genetics* 180(1):317-327.

Mohamad Ishak NS, Kato Y, Matsuura T, Watanabe H. 2016. Sequence conservation and sexually dimorphic expression of the Ftz-F1 gene in the crustacean *Daphnia magna*. *PLoS One* 11(5): e0154636.

Mohamad Ishak NS, Nong QD, Matsuura T, Kato Y, Watanabe H. 2017. Co-option of the bZIP transcription factor Vrille as the activator of

*Daphnia magna*. PLoS Genet. 13(11): e1006953.

Molinier C, Reisser CMO, Fields P, Haag CR, Galimov Y, Ségard A. 2018. Identification of general patterns of sex-biased expression in *Daphnia*, a genus with environmental sex determination. *G3 (Bethesda)* 8(5):1523-1533.

Murre C, Bain G, van Dijk MA, Engel I, Furnari BA, Massari ME, Matthews JR, Quong MW, Rivera RR, Stuiver MH. 1994. Structure and function of helix-loop-helix proteins. *BBA Gene Struct Expr.* 1218(2):129-135.

Nong QD, Mohamad Ishak NS, Matsuura T, Kato Y, Watanabe H. 2017. Mapping the expression of the sex determining factor Doublesex1 in *Daphnia magna* using a knock-in reporter. *Sci Rep.* 7(1):13521.

Olmstead AW, Leblanc GA. 2002. Juvenile hormone methyl farnesoate is a sex determinant in the crustacean *Daphnia magna*. *J Exp Zool.* 293(7):736-739.

Olmstead AW, LeBlanc GA. 2003. Insecticidal juvenile hormone analogs stimulate the production of male offspring in the crustacean *Daphnia magna*. *Environ Health Perspect.* 111(7):919-924.

Orsini L, Gilbert D, Podicheti R, Jansen M, Brown JB, Solari OS, Spanier KI, Colbourne JK, Rush D, Decaestecker E, et al. 2016. *Daphnia magna* transcriptome by RNA-Seq across 12 environmental stressors. *Sci Data* 3:160030.

Peichel CL, Nereng KS, Ohgi KA, Cole BLE, Colosimo PF, Buerkle CA, Schlüter D, Kingsley DM. 2001. The genetic architecture of divergence between threespine stickleback species. *Nature* 414(6866):901-905.

Reisser CMO, Fasel D, Hürlimann E, Dukić M, Haag-Liautard C, Thuillier V, Galimov Y, Haag CR. 2017. Transition from environmental to partial genetic sex determination in *Daphnia* through the evolution of a female-determining incipient W chromosome. *Mol Biol Evol.* 34(3):575-588.

Roulin AC, Routtu J, Hall MD, Janicke T, Colson I, Haag CR, Ebert D. 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.

Supek F, Bošnjak M, Škunca N, Šmuc T. 2011. Revigo summarizes and visualizes long lists of gene ontology terms. *PLoS One* 6(7):e21800.

Svendsen N, Reisser CMO, Dukić M, Thuillier V, Ségard A, Liautard-Haag C, Fasel D, Hürlimann E, Lenormand T, Galimov Y, et al. 2015. Uncovering cryptic asexuality in *Daphnia magna* by RAD sequencing. *Genetics* 201(3):1143-1155.

Tams V, Nickel J, Ehring A, Cordellier M, unpublished data, www.biorxiv.org/content/10.1101/503904v1, last accessed May 22, 2019.

Tatarazako N, Oda S, Watanabe H, Morita M, Iguchi T. 2003. Juvenile hormone agonists affect the occurrence of male *Daphnia*. *Chemosphere* 53(8):827-833.

Tokishita S-I, Kato Y, Kobayashi T, Nakamura S, Ohta T, Yamagata H. 2006. Organization and repression by juvenile hormone of a vitellogenin gene cluster in the crustacean, *Daphnia magna*. *Biochem Biophys Res Commun.* 345(1):362-370.

Toyota K, Gavin A, Miyagawa S, Viant MR, Iguchi T. 2016. Metabolomics reveals an involvement of pantothenate for male production responding to the short-day stimulus in the water flea, *Daphnia pulex*. *Sci Rep.* 6(April):25125.

Toyota K, Kato Y, Miyakawa H, Yatsu R, Mizutani T, Ogino Y, Miyagawa S, Watanabe H, Nishide H, Uchiyama I, et al. 2014. Molecular impact of juvenile hormone agonists on neonatal *Daphnia magna*. *J Appl Toxicol.* 34(5):537-544.

Toyota K, Miyakawa H, Hiruta C, Furuta K, Ogino Y, Shinoda T, Tatarazako N, Miyagawa S, Shaw JR, Iguchi T. 2015. Methyl farnesoate synthesis is necessary for the environmental sex determination in the water flea *Daphnia pulex*. *J Insect Physiol.* 80:22-30.

Toyota K, Miyakawa H, Yamaguchi K, Shigenobu S, Ogino Y, Tatarazako N, Miyagawa S, Iguchi T. 2015. NMDA receptor activation upstream of methyl farnesoate signaling for short day-induced male offspring production in the water flea, *Daphnia pulex*. *BMC Genomics.* 16(1):186.
Toyota K, Williams TD, Sato T, Tatarazako N, Iguchi T. 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.

Verhulst EC, van de Zande L, Beukeboom LW. 2010. Insect sex determination: it all evolves around transformer. *Curr Opin Genet Dev.* 20(4):376-383.

Verma R. 1996. Genetics of sex determination. Vol. 4. Elsevier, USA.

Weeks SC, Brantner JS, Astrop TI, Ott DW, Rabet N. 2014. The evolution of hermaphroditism from dioecy in Crustaceans: selfing hermaphroditism described in a fourth spinicaudatan genus. *Evol Biol.* 41(2):251-261.

Weirich GF, Svoboda JA, Thompson MJ. 1984. Ecdysone 20-monooxygenases. Berlin/Heidelberg (Germany): Springer. p. 227-233.

Xu SL, Zhou W, Chen P, Zhou JK, Zou X, Wang CL, Wang DL, Zhao YL. 2014. Identification and expression analysis of a *doublesex1* gene in

*Daphnia pulex* during different reproductive stages. *Dev Genes Evol.* 224(3):147-157.

Yatsu R, Miyagawa S, Kohno S, Parrott BB, Yamaguchi K, Ogino Y, Miyakawa H, Lowers RH, Shigenobu S, Guillette LJ, et al. 2016. RNA-seq analysis of the gonadal transcriptome during *Alligator mississippiensis* temperature-dependent sex determination and differentiation. *BMC Genomics*. 17(1):77.

Ye Z, Xu S, Spitze K, Asselman J, Jiang X, Ackerman MS, Lopez J, Harker B, Raborn RT, Thomas WK, et al. 2017. A new reference genome assembly for the microcrustacean *Daphnia pulex*. *G3 (Bethesda)* 7(5):1405-1416.

Zelzer E, Wappner P, Shilo BZ. 1997. The PAS domain confers target gene specificity of *Drosophila bHLH/PAS* proteins. *Genes Dev.* 11(16):2079-2089.