Skip to main content

A comprehensive epigenomic analysis of phenotypically distinguishable, genetically identical female and male Daphnia pulex



Daphnia species reproduce by cyclic parthenogenesis involving both sexual and asexual reproduction. The sex of the offspring is environmentally determined and mediated via endocrine signalling by the mother. Interestingly, male and female Daphnia can be genetically identical, yet display large differences in behaviour, morphology, lifespan and metabolic activity. Our goal was to integrate multiple omics datasets, including gene expression, splicing, histone modification and DNA methylation data generated from genetically identical female and male Daphnia pulex under controlled laboratory settings with the aim of achieving a better understanding of the underlying epigenetic factors that may contribute to the phenotypic differences observed between the two genders.


In this study we demonstrate that gene expression level is positively correlated with increased DNA methylation, and histone H3 trimethylation at lysine 4 (H3K4me3) at predicted promoter regions. Conversely, elevated histone H3 trimethylation at lysine 27 (H3K27me3), distributed across the entire transcript length, is negatively correlated with gene expression level. Interestingly, male Daphnia are dominated with epigenetic modifications that globally promote elevated gene expression, while female Daphnia are dominated with epigenetic modifications that reduce gene expression globally. For examples, CpG methylation (positively correlated with gene expression level) is significantly higher in almost all differentially methylated sites in male compared to female Daphnia. Furthermore, H3K4me3 modifications are higher in male compared to female Daphnia in more than 3/4 of the differentially regulated promoters. On the other hand, H3K27me3 is higher in female compared to male Daphnia in more than 5/6 of differentially modified sites. However, both sexes demonstrate roughly equal number of genes that are up-regulated in one gender compared to the other sex. Since, gene expression analyses typically assume that most genes are expressed at equal level among samples and different conditions, and thus cannot detect global changes affecting most genes.


The epigenetic differences between male and female in Daphnia pulex are vast and dominated by changes that promote elevated gene expression in male Daphnia. Furthermore, the differences observed in both gene expression changes and epigenetic modifications between the genders relate to pathways that are physiologically relevant to the observed phenotypic differences.


Daphnia (Crustacea: Cladocera) are fresh-water branchiopods, recognized as a model organisms by the U.S. National Institutes of Health [1]. Daphnia are used as a model organism in various fields of research, including ecotoxicology, ecology, population genetics and molecular studies [2,3,4,5]. Species of Daphnia typically reproduce by cyclical parthenogenesis. During the asexual phase female Daphnia produce genetically identical offspring [6]. When environmental conditions deteriorate (due to crowding, shortage of food or change in day-light cycle and temperature), Daphnia can switch to sexual reproduction, where female Daphnia produce both male and female offspring [7,8,9,10,11]. The female Daphnia produce haploid eggs which are fertilized by the male during mating to form diapausing resting eggs contained in an ephippium. These resting eggs can lay dormant in the sediment for prolonged periods of time, and hatch when environmental conditions improve [12,13,14].

The male and female offspring produced during the sexual reproduction are genetically identical in Daphnia [6], with sex being determined entirely by environmental factors, a system known as environmental sex determination (ESD). Daphnia offers unique opportunities in studying ESD, because the parthenogenetic female Daphnia can be maintained indefinitely in laboratory conditions via ameiotic reproduction to form clonal lineages and subjected to experimental manipulation [1]. The switch to male production can be manipulated either by altering the environment [11] or by administering methyl farnesoate (MF) or some other juvenile hormone analog [15, 16].

The genetically identical male and female Daphnia have a variety of morphological and behavioural differences, including lipid metabolism, mortality, and body size [17,18,19,20,21,22,23]. Previous studies have investigated gene expression differences between female and male Daphnia in several species [1, 24,25,26]. Despite differences in analysis techniques and quality of reference genomes these studies have identified substantial overlap in genes with sex-biased expression [26]. In this study, our aim was to further expand our understanding of the molecular differences between genetically identical female and male Daphnia that show clear phenotypic differences. Epigenetic factors, are known to contribute to phenotypic diversity in the absent of genetic differences [27, 28]. Therefore, we compared whole genome bisulphite sequencing (WGBS) data, histone modification data (H3K4me3 and H3K27me3) from chromatin immunoprecipitation sequencing, splicing and gene expression data collected from female and male Daphnia pulex under laboratory conditions.

Previous research of DNA methylation has shown that CpG-methylation is conserved among Daphnia species [29, 30]. We have also shown that in Daphnia and other arthropods high levels of DNA methylation within gene bodies as significantly correlated with elevated gene expression levels [30]. Since all of the previous studies on DNA methylation were carried out in female Daphnia, we wanted to see if DNA methylation was also conserved in male Daphnia, or if sex-specific differences could be observed, with correlated changes in gene expression and possibly alternative splicing. The application of ChIP-seq to study histone modifications (H3K4me3 and H3K27me3) is novel for Daphnia, but immunological studies have demonstrated that histone modifications do occur non-uniformly in Daphnia and are altered during development [31, 32].

This is the first comprehensive study that combines multiple epigenomics data with the aim of achieving a comprehensive understanding of the epigenetic differences between female and male Daphnia with environmental sex determination. Our data provides strong evidence that epigenetic markers are differently distributed between the two genders. Furthermore, it provides evidence in support of the hypothesis that epigenetic modifications may contribute towards an overall higher expression of majority of genes in male Daphnia compared to female Daphnia and this higher overall expression of genes in male Daphnia may contribute and explain some of the phenotypic differences observed between the two genders.


A multiomics approach was used to characterise the molecular profile of genetically identical female and male Daphnia pulex Eloise Butler strain. The aim of this study was to achieve a better understanding of sex-dependent molecular differences between genetically identical female and male D. pulex. To achieve this goal, the omics data (gene expression, ChIP-seq, DNA methylation and splicing data) were analysed both individually and in association with each other. This study provides the first insight into epigenetic and transcriptional differences between genetically identical genders of the model organism Daphnia that have evolved distinct morphological, physiological and behavioural differences.

Gene expression changes between male and female D. pulex

We analysed expression differences between male and female Daphnia pulex at the transcriptome and gene level. A significant expression difference (with Posterior Probability of Equivalent Expression: PPEE< 0.05) was observed in 11.2% (12,266/109,840) of the transcripts, which originate from 23.6% (7830/33,139) of the genes. The expression differences are symmetrically distributed, except for a slight excess of transcripts (55% higher in female Daphnia compared to 45% higher in male Daphnia) with higher expression in female Daphnia (Fig. 1a; Additional file 1: Table S1A).

Fig. 1

Differentially expressed transcripts between male and female Daphnia pulex (EB45) a) Volcano plot of differentially expressed transcripts. The transcripts marked with color are significantly different (Posterior Probablity of Equivalent Expression; PPEE < 0.05) between the sexes, (red = higher expression in female, blue = higher expression in male, pink = only expressed in female, light blue = only expressed in male). b) Reactome enrichment analysis for differentially (PPEE< 0.05) expressed transcripts. The enrichment analysis is carried out separately for the transcripts that have higher expression in male or female as well as for transcripts that are uniquely expression in one gender

The transcripts with higher expression in female Daphnia are enriched for RNA processing pathways (in particular rRNA) and translation, while the transcripts with higher expression in male Daphnia are enriched for muscle contraction, cardiac conduction, neuronal systems and cell signalling (Additional file 2: Table S2A). A small subset (13%) of transcripts (1614 transcripts in 1313 genes) are exclusively expressed in one gender. Half of these (805 transcript) are male specific (not expressed in female Daphnia), and half are female specific (809 transcripts; Fig. 1a). The transcripts that are uniquely expressed in female are not significantly enriched, and the male specific transcripts are enriched for the same pathways identified for the full set of differentially expressed transcripts (Additional file 2: Table S2A; Fig. 1b).

Most of the genes with differentially expressed transcripts were also differentially expressed when analysed at the gene level (71%; 5553/7830; Additional file 1: Table S1B), while a small subset of genes were differentially expressed only at the transcript level (either alternative splicing, alternative start or stop site usage) (Additional file 1: Table S1A; Additional file 1: Table S1B). The genes with only transcript level differences were enriched for the same pathways identified for the full set of differentially expressed transcripts (including RNA processing, muscle contraction and cell-cell communication; Additional file 2: Table S2A – S2C).

We detected 3291 potential splicing events using KisSplice (Additional file 1: Table S1C). The most common splicing event was intron retention (1244), followed by alternative acceptor and/or donator site usage (1142), with exon skipping being the third most common type (524). Very few splicing events (284) were significantly (FDR < 0.05) altered between male and female Daphnia. The splicing types were the same for the sex-specific events and all detected splicing events (chi-squared = 80, p value = 0.24), and they occurred mostly in the same genes that were already identified as having differentially expressed transcripts (80%; 226/284). The genes with detected sex-specific splicing changes were not significantly enriched for Reactome pathways (Additional file 2: Table S2D).

DNA methylation changes between male and female D. pulex

We performed whole genome bisulfite sequencing (WGBS) of Daphnia pulex Eloise Butler strain (genotypes EB31 and EB45). We quantified the methylation level of individual CpG sites (the ratio of methylated reads to read coverage at each site). The majority of CpG sites in Daphnia are unmethylated or have extremely low methylation level [29, 30, 33]. The high methylated CpGs (with median methylation level > 50%) are located mostly within exons (83%; 10,599/12,790 CpGs). Almost all of these (94.5%) are within the first four exons (with 1803, 4278, 2901 and 1074 CpGs in exons 1–4 respectively) in the primary transcripts, with exon 2 having the highest occurrence (40.4%) of high methylated CpGs. The primary transcripts containing high methylated CpGs (within exons 2–4) also have substantially higher expression level compared to the transcripts with only low methylated CpGs (Fig. 2).

Fig. 2

Density plot of mean expression levels (log2 FPKM) for genes that contain high methylated CpGs (> 50% median methylation; 2747 genes) and low methylated CpGs (< 50% median methylation; 33,139 genes) within exons 2–4 in the primary transcript in Daphnia pulex Eloise Butler (EB45)

After filtering out CpG sites with no methylated reads in more than half of the samples, only 18,951 sites remained for further analysis. The variation among the samples in the filtered CpG sites could be primarily attributed to differences between genotypes (EB45 vs EB31; PC1: 47% of variation) and sexes (female vs male; PC2: 41% of variation) (Fig. 3a). The CpG methylation level in male samples is globally higher than in female samples, with more than 70% of all CpGs having higher methylation level in male compared to female samples (Fig. 3b). A statistically significant difference in the methylation levels in CpG sites (FDR < 0.05) was observed for 1841 CpGs (9.71%) between male and female Daphnia (Additional file 1: Table S1D). The differentially methylated CpGs (DMC) are located within gene bodies (97.56%; 1796/1841), and in particular within the first four exons (78.67%; 1413/1796). Very few DMCs are located outside of known genes (2.4%; 45/1841) (Additional file 1: Table S1D) and almost all of the DMCs have higher methylation level in male Daphnia (96.46%, 1776/1841 DMCs) compared to female Daphnia (Fig. 3b).

Fig. 3

DNA methylation differences between male and female Daphnia pulex in Eloise Butler strain (genotypes EB31 and EB45), using a filtered dataset; CpGs not covered in all samples and with no methylated reads in more than half of the samples were excluded. a) Principal component analysis (PCA) of DNA methylation (CpG) levels. Samples are represented by points along PC1 (x-axis) and PC2 (y-axis), which account for the majority of variance in the data. Genotypes separate by PC1 which accounts for 47% of the variance in methylation and the sexes separate by PC2 which accounts for 41% of variance. b) Volcano plot of DNA methylation (CpG) differences between male and female. The differentially methylated CpGs (DMCs; FDR < 0.05) are indicated in red

The DMCs with higher methylation in male Daphnia are not significantly enriched for any known pathways (Additional file 2: Table S2E). This potentially indicates that the higher methylation of genes in male Daphnia compared to female Daphnia is non-specific and global. The few genes with lower methylation level in male Daphnia compared to female Daphnia are however enriched for specific cellular functions, including cellular senescence, interleukin-17 signalling and negative regulation of FGFR signalling (Additional file 2: Table S2E). The transcripts containing DMC with decreased methylation in male Daphnia also demonstrate a reduced expression compared to female Daphnia for ~ 80% of the transcripts (Fig. 4), while the DMCs with increased methylation in male Daphnia have no association with expression level at the transcript level.

Fig. 4

Heatmap of expression and DNA methylation for the transcripts that contain differentially methylated CpGs, where the methylation is significantly lower in male compared to female Daphnia. The expression and methylation levels were scaled from 0 to 1, with red indicating high expression or high methylation and blue low expression or low methylation. The sidebar shows the average direction of expression change, with red indicating increased expression and blue decreased expression in female compared to male Daphnia

Histone modification changes between male and female D. pulex

The initial ChIP-peaks identified with MACS2 are substantially smaller for histone H3 trimethylated at lysine 27 (H3K27me3; with a median size of 318 bp) compared to histone H3 trimethylated at lysine 4 (H3K4me3; 800 bp). During the differential peak analysis (DiffBind) overlapping peaks are merged resulting in slightly larger peaks (488 bp for H3K27me3 and 968 bp for H3K4me3). The H3K4me3 peaks are more sparsely located in the genome with a 3089 bp median distance between peaks, compared to 430 bp for H3K27me3 (with long stretches of nearby peaks). The peak intensity (ChIP compared to input controls) for H3K4me3 is higher than for H3K27me3, with a median fold change of 5.15 vs 2.02 in the initial peak discovery, and 7.08 vs 4.95 in the differential peak analysis for H3K4me3 and H3K27me3, respectively (Additional file 1: Table S1E; Additional file 1: Table S1F).

We identified 10,092 H3K4me3 peaks, 95% (9602) of which are consistently found (FDR < 0.05) in all samples (n = 6) compared to input controls (n = 2) (Additional file 1: Table S1E). Almost all (97%; 9365) of these peaks are within 200 bp of known genes (10,968 genes, with some peaks overlapping with more than one gene), and enriched at the start of the gene, with 90% (8438) overlapping with exon 1. About 10% (1061) of the H3K4me3 peaks have sex-specific differences in intensity (FDR < 0.05), with 78% (833) of the sex-specific peaks having higher intensity in male Daphnia (in 1068 genes) and 22% (228) having higher intensity in female Daphnia (in 275 genes) (Fig. 5a). The genes with higher H3K4me3 intensity in female Daphnia compared to male Daphnia are enriched for multiple Reactome pathways, including collagen formation, lipid metabolism, heme biosynthesis, extracellular matrix organization and cell motility via c-Met signalling pathway. Whereas the genes with higher H3K4me3 intensity in male Daphnia are only marginally enriched for cardiac conduction and related pathways (Fig. 5c; Additional file 2: Table S2F).

Fig. 5

Differentially regulated histone modifications between male and female Daphnia pulex. A) Volcano plot for H3K4me3, B) Volcano plot for H3K27me3 between male and female. The differentially modified histone peaks (FDR < 0.05) are indicated in red. C) Reactome enrichment analysis of differential histone modifications analysed separately for transcripts that have higher peak intensity in male or female

We identified almost three times as many peaks (29,162) for H3K27me3, compared to H3K4me3. Similar to H3K4me3, most of the peaks (97%) are consistently found (28,372/29,162; FDR < 0.05) in all samples compared to the input controls, and are associated (99%; 28,284 peaks) with known genes (12,901 genes; Additional file 1: Table S1F). Overall, 41% (12,102) of the H3K27me3 peaks (in 7329 genes) had different intensities (FDR < 0.05) between male and female Daphnia. In contrast to the gene expression promoting H3K4me3 histone modification, the expression suppressing H3K27me3 modification occurred predominantly (> 86%; 10,356) in female Daphnia (in 6123 genes), while only 14% (1753) of the H3K27me3 peaks had higher intensity in male Daphnia (in 1296 genes) (Fig. 5b). The genes with higher H3K27me3 intensity in female compared to male Daphnia are enriched for multiple Reactome pathways including GPCR signalling, transport of small molecules, G-Protein alpha-i signalling, digestion, muscle contraction and neuronal systems. Whereas the genes with higher H3K27me3 intensity in male Daphnia are not significantly enriched for any Reactome pathways (Fig. 5c; Additional file 2: Table S2G).

The histone modifications have significant association with gene expression. Genes with H3K4me3 modifications have two times higher mean expression (FPKM 31.97 vs 15.95), compared to genes without H3K4me3 modifications (Fig. 6a). The opposite pattern is observed for H3K27me3 modifications. Genes with H3K27me3 modifications have two times lower mean expression (FPKM 14.20 vs 24.28), compared to genes without H3K27me3 modifications (Fig. 6b). While genes containing both modifications have an intermediate expression level (Fig. 6c).

Fig. 6

Expression densities for genes with or without histone modifications. a H3K4me3, b) H3K27me3, c) both H3K4me3 and H3K27me3. The expression level (FPKM) is averaged across all samples and log2-tranformed.

Integrative analysis: covariation and co-occurrence

The DNA methylation and histone modifications affect gene expression in an additive manner (Fig. 7a). DNA methylation (in exons) increases gene expression (from mean FPKM 18.17 to 32.21) regardless of histone modifications. The presence of H3K4me3 in methylated genes additionally increase expression (to FPKM 40.25), while H3K27me3 decreases expression (to FPKM 11.62). The highest expression is observed in genes that have both DNA methylation, contain H3K4me3 and are absent of H3K27me3 modifications (mean FPKM 41.59). While the lowest expressed genes are absent of all modifications. The very low expressed genes undoubtedly contain genes with mapping problems (highly variable or partial genes), which could result in reduced detection in all datasets.

Fig. 7

Combined comparison of DNA methylation, histone modifications and gene expression. a Violin plot of gene expression separated by presence/absence of DNA methylation and histone modifications: H3K4me3 and H3K27me3. The average gene expression in all samples, Fragments Per Kilobase of transcript per Million (FPKM) is log2 scaled. b Venn diagram of genes with DNA methylation and histone modifications, for all genes with detectable modifications above the filtering cutoffs specified in the methods. c Venn diagram for genes, with significant differences (FDR < 0.05) between male and female Daphnia pulex for the modifications. d Heatmap of ranked values for gene expression (FPKM), histone modifications (H3K4me3 and H3K27me3), and DNA methylation (CpG). Red indicates high level of expression or modification, blue indicates low level of expression or modifications. Genes separate into 5 main clusters by omics profile. e Enrichment results for the most significant Reactome pathways in the main clusters from the heatmap (1–5)

The majority of the genes containing DNA methylation (69.19%) also contain H3K4me3 histone modifications (chi-squared = 7615.5, p value = 2.9e− 1656), which is more than twice the value one would expect by chance (5346 genes observed with both modifications compared to 2281 genes expected by chance). While the overlap among genes with H3K27me3 and DNA methylation (obs: 2543 vs exp.: 2759; chi-squared = 34.1, p value = 5.2e− 09) or H3K27me3 and H3K4me3 (obs: 2181 vs exp.: 3493; chi-squared = 1087.1, p value = 2.1e− 238) is significantly underrepresented considering the huge number of genes containing these modifications (Fig. 7b).

While the overlap is substantially smaller for the genes where these modifications are different between male and female Daphnia (Fig. 7c), the overlap is still significantly different than one would expect by chance. The overlap between DNA methylation and H3K4me3 is significantly enriched (111 genes with both modifications compared to 41 expected by chance; chi-squared = 123.7, p value = 1.0e− 28) as is the overlap between H3K4me3 and H3K27me3 (obs: 326 vs exp.: 188; chi-squared = 128.9, p value = 7.0e− 30). The overlap between DNA methylation and H3K27me3 is significantly underrepresented (obs: 211 vs exp.: 271; chi-squared = 16.8, p value = 4.1e− 05).

Most genes with sex-specific differences contain a single modification, especially when contrasted to the global background of DNA methylation and histone methylations, where the overlap is substantial. The few sex-specific genes that contain multiple modifications are not significantly enriched to any particular pathway after multiple testing corrections (Additional file 2: Table S2H). Granted each of these omics methods has their own unique strengths and weakness (unique error profiles, statistical power and heterogeneity among samples), which could partially explain the low level of overlap. On the other hand, the global background of co-occurring DNA methylation and histone modifications (Fig. 7d) are significantly enriched for multiple pathways (Fig. 7e; Additional file 2: Table S2I). Many of the same pathways are enriched in the individual analyses for DNA methylation, histone modification and gene expression for the sex-specific differences (Table 1).

Table 1 Summary of enriched Reactome pathways across multiple omics datasets comparing male and female Daphnia pulex

The sex-specific changes in gene expression, DNA methylation, histone modifications and alternative splicing are evenly distributed across the genome (scaffolds assigned to chromosomes according to Ye et al. 2017) (Fig. 8), with slight excess from an expected distribution in chromosomes 9 and 11 for DNA methylation, H3K4me3 and H3K27me3 and chromosome 4 for gene expression and alternative splicing.

Fig. 8

a) Density plot of the epigenetic modification. Showing the relative locations of histone modifications (H3K4me3 and H3K27me3) and CpG methylations (in different colour), scaled to the start and end of the primary transcript of each gene. Extremely short and long genes are excluded (transcript length below 1000 bp or above 10,000 bp). The modifications are mapped to the nearest gene, taking the relative distance to the start of the gene from the peak maximum separately for each sample (genders are indicated with line type). b Circos-plot of sex-specific differences in the multiple omics datasets, distributed across the genome. The scaffolds assignment to chromosomes is based on [34]. The direction of change is indicated with color; blue = higher in male, pink = higher in female. The differentially spliced genes are indicated in purple and the alternatively spliced genes that do not display sex-specific difference are indicated in green


Sex determination, a fundamental biological process, impact the development of most organs and causes sex-specific differences in behaviour, physiology and morphology [35]. Sex determination in majority of organisms is underpinned and regulated either by genetic factors (GSD: genetic sex determination) or environmental factors (ESD: environmental sex determination). The latter is initiated by environmental cues, such as temperature, nutrition, population density and photoperiod. ESD is observed in range of species across the animal taxa, such as rotifers, nematodes, crustaceans, insects, fishes and reptiles [35]. The crustacean Daphnia are also subject to environmental sex determination [36], where by the same genotype may develop into phenotypically distinct male and female Daphnia depending on the environmental cues [26, 37]. The genetically identical male and female Daphnia demonstrate differences in their phenotype and life-history traits, including metabolic activity, mortality, morphology (antenna and carapace) and body size [17, 18]. In particular the female Daphnia are larger, almost twice the size of male Daphnia, while the male Daphnia have a higher metabolic rate and shorter lifespan compared to female Daphnia [21,22,23]. Once the sex is determined it is maintained throughout the life of the organism even in the absence of the initial environmental cue [11, 16]. The maintenance of the acquired gender throughout the life of an organism can be caused by early developmental changes which result in a cascade of differences including structural alterations. It can also include regulatory factors such as hormones that need to be constantly maintained at specific levels. Such regulatory factors can also include epigenetic factors that help to maintain the acquired phenotype, leading to creating a sex-specific molecular signature. Our goal in this study was to achieve a better characterisation and understanding of the sex-specific differences (signature) at a molecular level with a specific focus on contribution of epigenetic factors (histone modifications and DNA methylation). To achieve this objective, we generated omics data at multiple levels to create a molecular signature for female and male Daphnia.

Previous studies have investigated differences in the expression levels of genes between female and male Daphnia (in D. pulex, D. magna and D. galeata) [1, 24,25,26]. Our study differs from previous published work as in addition to investigate differences in expression at the gene level, we also investigated changes in expression at the transcriptome level identifying variation in alternative splicing and usage of alternative start and stop sites. Our data indicated that the genes with the same basal expression level containing differentially expressed alternative isoforms between female and male Daphnia were enriched for RNA processing pathways and translation regulation. However, the genes with detected splicing variations were not significantly enriched for Reactome pathways. The alternative splice variants that are differentially regulated between the sexes may represent a diverse set of tissue specific changes, in line with the morphological differences between the sexes.

Our results, similar to previous findings, demonstrate that a large portion of genes display significant differences in the expression between male and female Daphnia, affecting more than 1/5 of all annotated genes. We further compared our list of sex-specific genes to D. magna [26]. The two species D. magna and D. pulex are among the most distantly related Daphnia species and span the entire phylogeny of the genus, having diverged more than 200 million years ago [38]. In D. magna, 42% of the genes are reported as differentially expressed between male and female [26], which is substantially higher than what we detected in D. pulex (~ 20%). Out of the 11,197 differentially expressed genes in D. magna, we could find a reliable ortholog in D. pulex for 7920 genes (using blastp with e-value <1e− 20). The agreement between D. pulex and D. magna for the identified 7920 sex-specific genes is substantial. Without filtering the data based on significance, in D. pulex > 73% of the genes have the same direction of expression change as in D. magna. When selecting only the genes that we detected as significantly differentially expressed (3093/7920 genes) the agreement increases to > 86%. Not only is the direction of change the same but also the magnitude of expression change is highly correlated (R2 = 0.55, p value < 2.2e− 16), especially for the genes with higher expression in female Daphnia (Additional file 3: Fig. S1). This potentially indicates that sex-specific genes and the enriched pathways (e.g. RNA metabolisms, signalling and development) are conserved between the two species and are essential for maintaining sex-specific characteristics.

It is worth highlighting that these conserved genes included known sex-determination factors. For example, in Daphnia there are several orthologs for the Drosophila doublesex (dsx) gene, which are not alternatively spliced as in insects, but regulate sex determination by expression level [39]. In Daphnia magna two of these genes (DapmaDsx1: APZ42_027481, DapmaDsx2: APZ42_027480) have elevated expression in male Daphnia, with DapmaDsx1 being capable of regulating male morphology when ectopically applied and female traits when knocked-down during embryogenesis [35]. The Daphnia pulex orthologs of DapmaDsx1 (Daplx7pEVm013292) and DapmaDsx2 (Daplx7pEVm013921) both have higher expression in male Daphnia (log2FC = − 4.02 and log2FC = − 6.18, respectively, with PPEE< 2.2e− 16 for both; Additional file 3: Fig. S1), and also contain significant differences in the modification H3K4me3, with higher level in male Daphnia (log2FC = − 8.25, FDR = 1.20e− 25 and log2FC = − 4.98, FDR = 7.87e− 06), whereas female Daphnia had higher level of H3K27me3 modification in both genes (log2FC = 12.40, FDR = 1.79e− 40 and log2FC = 13.34, FDR = 1.87e− 54; Additional file 1: Table S1).

Histone modifications can rapidly regulate the expression of genes [40, 41]. In this study, we analysed two histone modifications, H3K4me3 and H3K27me3, known to regulate the expression of genes in a variety of species [42, 43]. H3K4me3 modification is a hallmark of actively transcribed genes and it is commonly associated with transcription start sites (TSS) and promoter regions [44], whereas H3K27me3 peaks at the TSS and promoter region, it is more spread out along the length of the affected genes than the H3K4me3 modification. Furthermore, H3K27me3 is strongly associated with down-regulation of nearby genes via the formation of heterochromatic regions [45]. Both active and inactive modifications can be found in Daphnia in the expected locations (Fig. 8a). The H3K4me3 modifications were concentrated at start of the genes, with 97% of the detected peaks within 200 bp of the known transcription start site. While H3K27me3 modifications occurred throughout the gene body and intergenic regions. The majority of the histone modification peaks were observed in both male and female Daphnia. The effect of the histone modifications on gene expression level was clear and in line with the expectations (Fig. 6) with H3K4me3 modification promoting higher expression level and H3K27me3 modification generally suppressing the expression of the genes. Most interestingly, we observed that the majority of sex-specific H3K4me3 peak are higher in male Daphnia (78%), while female Daphnia are dominated by higher H3K27me3 peaks (86%). This difference can potentially indicate a higher basal level of global expression in male compared to female Daphnia. We also detected a relatively small number of genes where both modifications were present (Fig. 7b). This resulted in an intermediate expression level (Fig.6c) potentially creating genes in a poised status ready to be either expressed or suppressed (higher expression compared to genes with only H3K27me3 and lower than genes with only H3K4me3) [46,47,48]. However, the latter category requires further investigation to remove the possibility of mix peak signal due to presence of multiple cell populations.

In addition to histone modifications, we investigated the differences in CpG methylation between the two genders. Similar to our previous findings, the majority of the methylated CpG sites in both genders were located within the gene body and mostly concentrated at exons 2–4 region [30]. Genes with high levels of CpG methylation (> 50%) in both genders demonstrated an elevated expression level compared to rest of the genes (Fig. 2; similar to Kvist et al., 2018). Furthermore, based on our data, the two epigenetic modifications of CpG methylation and H3K4me3 demonstrated a complementary and additive effect on gene expression. As demonstrated in Fig. 7a, genes with both modifications had a significantly higher expression level compared to the rest of the genes. Most interestingly, CpG methylation levels are overall significantly higher (96% of all DMC) in male compared to female Daphnia. This observed non-specific global higher level of methylation in male Daphnia coupled with higher H3K4me3 peaks in male compared to female Daphnia could further suggest a potential basal global higher gene expression in males. However, at the gene expression level there is no obvious bias in male Daphnia demonstrating a higher expression for majority of the genes compared to female Daphnia. In fact, there are slightly more genes (5% more) with higher expression in female compared to male Daphnia. Although our data does not support a male biased higher gene expression level, the existence of such bias in gene expression cannot be entirely ruled out at this stage as methods used for normalising the data, library preparation and RNA-seq analysis can mask global biases [49]. In order to evaluate if a global bias in gene expression truly exists between male and female Daphnia, one would need to use external spike-in references during sample preparation, which would tie the cell counts to mRNA yields, and permit absolute quantification of gene expression. The traditional normalization methods used (in this study and all other Daphnia gene expression studies) assume that most genes are expressed at the same level among samples, and cannot detect global biases that affect all or most genes [49]. An alternative explanation is that the lack of male biased gene expression level, which is observed at histone modification and CpG methylation levels, could be real. It is feasible that there are compensatory changes in female Daphnia (besides the ones studied here) that balance, and slightly (5% of genes) increase, the level of gene expression between female and male Daphnia. For example in mouse lymphocytes, elevated expression of a single transcription factor (c-myc) can result in a global transcriptional amplification of all actively transcribed genes [50]. The Daphnia pulex ortholog of c-myc (Daplx7pEVm006187) was indeed elevated in female Daphnia pulex in this study (log2FC = 1.39 higher in female compare to male, PPEE< 2.2e− 16). As well as in D. magna (APZ42_014785) in another study (log2FC = 0.64 higher expression in female compared to male, adjusted p value = 5.3e− 05) [26].

Enrichment analysis demonstrated that genes with higher CpG methylation and histone modifications in male Daphnia were not enriched for specific pathways and were mostly randomly distributed across the genome. In contrast, genes containing higher CpG methylation levels in female Daphnia were enriched for partially linked pathways related to immune response (Toll like receptor cascades, Interleukin-17 signalling, Class I MHC mediated antigen processing & presentation, and TRAF6 mediated induction of NFkB and MAP kinases upon TLR7/8 or 9 activation) and ageing (Cellular senescence, Senescence-Associated Secretory Phenotype, MAP kinase activation, and Negative regulation of FGFR signalling). Enrichment of these particular pathways in female Daphnia may be related to the fact that female Daphnia typically have a longer lifespan compared to male Daphnia [21,22,23], although few male strains maintained under specific conditions have shown to outlive females [51]. The enriched pathways could explain some of the phenotypic differences observed between female and male Daphnia. For examples, the heat shock response protects the cells against a plethora of external and internal damage, including elevated temperature, oxidative damage, metal stress and also ageing related protein misfolding and aggregation [52, 53]. Heat shock proteins (HSPs) can also activate innate immune system [54]. HSPs are differentially expressed between sexes in Daphnia, with most HSPs having higher expression in female Daphnia. Also HSPs react more strongly to heat stress in female Daphnia [55]. In comparisons among Daphnia species elevated HSP expression is associated with longer lifespan [56]. We observed 80% of the differentially expressed heat shock proteins (11/14 genes) having higher expression in female compared to male Daphnia, including heat shock transcription factor 1 (HSF1; Daplx7pEVm005655, log2FC = 0.52), despite HSF1 having (9.43%) higher methylation level in male Daphnia.

Male Daphnia grow more slowly compared to female Daphnia and reach a smaller body size [17, 18]. Female Daphnia accumulate lipids they acquire from their food [19], which are used for producing eggs (sexual and asexual) [57, 58]. These morphological differences are in line with the enrichment results for the relatively few genes that had higher H3K4me3 levels in female Daphnia (Metabolism of lipids, Biological oxidations and Heme biosynthesis). Male Daphnia are typically smaller than female Daphnia, are more active, and faster swimmers [20], have faster heartbeat rate [22] and in general have higher metabolic activity compared to female Daphnia. These differences are reflected in the patterns of gene expression with enriched pathways for muscle activity (Ion homeostasis, Muscle contraction and Cardiac conduction) for genes with higher expression in male compared to female Daphnia (Additional file 2: Table S2F).


Overall, our study indicates that genetically identical female and male Daphnia have evolved distinct DNA methylation, histone modification and gene expression patterns which could explain the differences in morphology, physiology and behaviour between male and female Daphnia. As discussed, some of the changes observed at the gene (doublesex genes and HSP genes) and pathway (cellular senescence pathway and immune response) levels support this hypothesis. Furthermore, this is the first multi-omics study that provides insight into interactions between histone modifications (H3K4me3 and H3K27me3), DNA methylation and gene expression in any Daphnia species. We demonstrate the impact of the two histone modifications and DNA methylation individually, and more interestingly when they co-occur, on gene expression. Finally, this study provides further evidence in support of use of Daphnia as a model organism for research into epigenetic regulation of traits and phenotypic plasticity.


Daphnia pulex maintenance and induction of males

Cultures of Daphnia pulex Eloise Butler strain (genotypes EB31 and EB45, originally sampled from Eloise Butler pond in Minnesota, [59] were maintained in standard COMBO as previously described [30, 60, 61]. To induce male Daphnia, sexually mature individual female Daphnia were treated with the crustacean reproductive hormone, methyl (2E, 6E)-farnesoate (MF) at a final concentration of 400 nM. This concentration is sufficient to induce male Daphnia at 100% efficiency [16]. Due to the instability of MF, medium was changed daily to ensure consistent exposure. The first brood was discarded, and male neonates were collected from 2nd – 3rd broods. Female Daphnia used in the ‘omics studies were not exposed to MF. Similar to the male samples, neonates from 2nd-3rd broods were collected and used in this study. Female and male cultures were maintained separately.

DNA and RNA extraction and sequencing

Genomic DNA and RNA were extracted from a pool of samples with a mixture of different ages (3, 8 and 15 days old) using MasterPure DNA purification kit (Epicentre, USA) and RNeasy Micro Kit (Qiagen Ltd., UK), respectively as described by Athanasio et al. 2016 and 2018 [61, 62]. DNA for the whole genome bisulfite sequencing (WGBS) was extracted from both genotypes (EB31 & EB45), from 3 female and 3 male Daphnia pools from each genotype. The ChIP-seq and RNA-seq samples were prepared from only one genotype (EB45). DNA for the ChIP-seq was extracted from 3 female, 3 male and 2 input control pools. RNA for the gene expression and splicing analysis was extracted from 2 female and 3 male Daphnia pools. The whole genome bisulfite sequencing (WGBS) libraries and the RNA sequencing libraries (RNA-seq) were prepared as described in our previous publication [30]. Briefly, the EpiGenome Methyl-Seq kit (Epicentre, USA) was used to prepare the WGBS libraries and sequenced (2x80bp) using Illumina NextSeq 500 platform at the Centre for Genomics and Bioinformatics, Indiana University. The RNA-seq libraries were prepared using the Illumina TruSeq standard mRNA kit and sequenced (1x85bp) using Illumina NextSeq 500 platform at the Centre for Genomics and Bioinformatics, Indiana University. The chromatin immunoprecipitation sequencing libraries (ChIP-seq) were prepared using the iDeal-seq kit, H3K4me3 (C15410003–50, 1 μg/reaction), H3K27me3 (C15410195, 1 μg/reaction) antibodies and sequenced using Illumina HiSeq 2500 (1 × 50 bp) as part of a service provided by Diagenode (Belgium). Briefly, Daphnia samples (30 mg wet tissue per sample) were homogenised in 1 ml of PBS/1%formaldehyde using Dounce homogenizer. The collected cells were lysed and the nuclei were collected and sonicated to a final size of 80–400 bp. The mentioned antibodies were used to prepare test samples according to the manual for the iDeal ChIP-seq kit. The IP samples and input samples were quantified using the Qubit dsDNA HS kit. Library preparation was performed on the IP and input samples using the MicroPLEX library preparation protocol on 500 pg of DNA. The amplified libraries (13 PCR cycles) were purified using AMPure beads, quantified using the Qubit ds DNA HS kit and analysed on Bioanalyzer. The prepared libraries were then sequenced on HiSeq 2500. This project has been deposited at NCBI GEO under accession GSE12442.

Pre-processing, mapping, preliminary analysis

lllumina adapters (using core sequence: AGATCGGAAGAGC) and nucleotides with low quality (Phred score < 20) were removed with cutadapt (v.1.11) [63]. The filtered reads were mapped to the reference genome of Daphnia pulex PA42 (GCA_900092285.1) [34] using BWA Meth (v.0.10) [64] for bisulfite-treated DNA samples, BWA-MEM (v.0.7.15-r1140) [65] for the non-bisulfite treated DNA samples (ChIP-seq and reference DNA), and with RSEM (v.1.3.0) [66] using STAR aligner (v.2.5.3a) [67] for the RNA-seq samples, with default settings. The Daphnia pulex gene models used in the analysis are from November 2017 obtained from the arthropod database in eugenes (Genomic Information for Eukaryotic Organisms; produced by Don Gilbert using EvidentialGene [68].

Analysis of gene expression and splicing data

Expression changes were analysed at gene and transcript levels using EBSeq (v.1.20.0) [69], with default settings. Genes and transcripts with significant expression difference between male and female Daphnia (with posterior probability of differential expression < 0.05) were analysed further. An additional alternative splicing analysis was conducted on the same filtered reads used for the expression analysis, using the de novo splicing predictor, KisSplice (v2.4.0-p1) [70] with default settings. The potential splicing events detect by KisSplice (type_1) were mapped back to the reference genome (GCA_900092285.1) with STAR aligner (v2.5.2a) [67], using default settings. The mapping results were analysed with KisSplice2RefGenome (v.1.0.0) [71] to identify the types of splicing events that occurred in the samples. Alternative splicing events were analysed for sex induced (male vs female) differential changes with kissDE (v1.5.0) [71]. Splicing events that did not map to known genes or mapped to multiple locations as well as events that were low coverage were excluded. Splicing events that were insertions, deletions or SNPs according to the genomic mapping were also removed.

Analysis of DNA methylation data

Differential methylation analysis was done using methylKit (v.1.3.0) [72]. CpG sites with abnormally high (> 98 percentile) coverage were removed, as well as sites that were not covered in all samples or had zero methylated reads in more than half of the samples (n = 6/12). Logistic regression was used to analyse differential CpG methylation between male and female, using genotype (EB31 and EB45) as a co-variable. The Q-values were adjusted using the SLIM method [73].

Analysis of chromatin immunoprecipitation sequencing data

The DNA fragments containing histone modification (H3K4me3 and H3K27me3) were purified, sequenced and aligned to the genome. The ChIP-seq reads were filtered by mapping quality (MAPQ > 30) to reduce background noise from unspecific mapping. The genomic locations where the DNA fragments were concentrated (peaks) were identified. The peaks corresponding to histone modifications (H3K4me3 and H3K27me3) were called with MACS2 (v. [74], separately for each sample without sifting model building using 132Mbp as an estimate of the mappable genome size and predicted fragment sizes 134 bp (for H3K4me3) and 144 bp (for H3K27me3) as estimated from the data. Differential analysis of histone peaks (narrowPeak) were achieved using DiffBind (v.2.8.0) [75], by comparing the male and female samples against each other (n = 3 for both sexes and histone modifications) and against the input controls (n = 2). The peaks for H3K27me3 were mapped to the nearest transcript, and the peaks for H3K4me3 were mapped against the nearest exon 1. Differential peaks (FDR < 0.05) within 200 bp of known transcripts (H3K27me3) or exon 1 (H3K4me3) were retained for further analysis.

Enrichment analysis

The differentially regulated (FDR < 0.05) genes (containing CpG methylation, modified histones, expression or splicing changes) were analysed for enrichment in Reactome pathways [76] with ClusterProfiler (v.3.8.1) [77] and ReactomePA (v.1.24.0) [78]. Since Daphnia pulex genes are not annotated in Reactome, we used protein blast (with e-value <1e− 20) to identify orthologous genes in humans. The reference genes (universe) for the enrichment analysis were limited to only those human genes that were identified by blast and had NCBI gene IDs (9992 Daphnia pulex genes, matching to 6013 unique genes). 40% (4014) of these genes were annotated in the Reactome database.

Availability of data and materials

This project has been deposited at NCBI GEO under accession GSE12442. The reference genome and chromosomal assignment of scaffolds for Daphnia pulex is based on Ye et al. 2017 (DOI: The Daphnia pulex gene models are from the arthropod database in eugenes (Genomic Information for Eukaryotic Organisms) produced by Don Gilbert using EvidentialGene (DOI: Expression data for Daphnia magna sex-biased genes are from Molinier et al. 2018 (DOI:



Chromatin Immunoprecipitation


Differentially methylated CpGs


Doublesex gene


Environmental Sex Determination


False discovery rate


Fragments per kilobase of transcript per million mapped reads


Histone H3 trimethylated at lysine 27


Histone H3 trimethylated at lysine 4


Heat Shock Proteins


Methyl Farnesoate


Whole genome bisulfite sequencing


  1. 1.

    Colbourne JK, Pfrender ME, Gilbert D, Thomas WK, Tucker A, Oakley TH, et al. The ecoresponsive genome of Daphnia pulex. Science. 2011;331:555–61.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  2. 2.

    Orsini L, Gilbert D, Podicheti R, Jansen M, Brown JB, Solari OS, et al. Daphnia magna transcriptome by RNA-Seq across 12 environmental stressors. Sci Data. 2016;3:160030.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  3. 3.

    Harris KDM, Bartlett NJ, Lloyd VK. Daphnia as an emerging epigenetic model organism. Genet Res Int. 2012;2012:147892.

    Article  Google Scholar 

  4. 4.

    Lynch M, Gutenkunst R, Ackerman M, Spitze K, Ye Z, Maruki T, et al. Population Genomics of Daphnia pulex. Genetics. 2017;206:315–32.

    CAS  Article  PubMed  Google Scholar 

  5. 5.

    Shaw JR, Pfrender ME, Eads BD, Klaper R, Callaghan A, Sibly RM, et al. Daphnia as an emerging model for toxicological genomics. Adv Exp Biol. 2008;2.

    Chapter  Google Scholar 

  6. 6.

    Hebert PDN, Ward RD. Inheritance during parthenogenesis in Daphnia magna. Genetics. 1972;71:639–42

    CAS  PubMed  PubMed Central  Google Scholar 

  7. 7.

    Brown LA, Banta AM. Control of sex in Cladocera: VIII. Food and crowding as factors in male production. Physiol Zool. 1935;8:138–55.

    Article  Google Scholar 

  8. 8.

    Stross RG, Hill JC. Photoperiod control of winter diapause in the fresh-water crustacean. Daphnia Biol Bull. 1968;134:176–98.

    Article  Google Scholar 

  9. 9.

    Stuart CA, Cooper HJ, Miller R. A critical analysis of excretory products as sex controlling agents in cladocera. Wilhelm Roux Arch Entwickl Mech Org. 1932;126:325–47.

    CAS  Article  PubMed  Google Scholar 

  10. 10.

    Korpelainen H. The effects of temperature and photoperiod on life history parameters of Daphnia magna (Crustacea: Cladocera). Freshw Biol. 1986;16:615–20.

    Article  Google Scholar 

  11. 11.

    Hobæk A, Larsson P. Sex determination in Daphnia magna. Ecology. 1990;71:2255–68.

    Article  Google Scholar 

  12. 12.

    Schwartz S, Hebert P. Methods for the activation of the resting eggs of Daphnia. Freshw Biol. 1987;17:373–9.

    Article  Google Scholar 

  13. 13.

    Orsini L, Schwenk K, De Meester L, Colbourne JK, Pfrender ME, Weider LJ. The evolutionary time machine: using dormant propagules to forecast how populations can adapt to changing environments. Trends Ecol Evol. 2013;28:274–82.

    Article  PubMed  PubMed Central  Google Scholar 

  14. 14.

    Cáceres CE. Interspecific variation in the abundance, production, and emergence of Daphnia diapausing eggs. Ecology. 1998;79:1699–710.[1699:IVITAP]2.0.CO;2.

    Article  Google Scholar 

  15. 15.

    Oda S, Tatarazako N, Watanabe H, Morita M, Iguchi T. Production of male neonates in four cladoceran species exposed to a juvenile hormone analog, fenoxycarb. Chemosphere. 2005;60:74–8.

    CAS  Article  PubMed  Google Scholar 

  16. 16.

    Olmstead AW, Leblanc GA. Juvenoid hormone methyl farnesoate is a sex determinant in the crustacean Daphnia magna. J Exp Zool. 2002;293:736–9.

    CAS  Article  PubMed  Google Scholar 

  17. 17.

    MacArthur JW, Baillie WHT. Metabolic activity and duration of life. I. Influence of temperature on longevity in Daphnia magna. J Exp Zool. 1929;53:221–42.

    Article  Google Scholar 

  18. 18.

    Munro IG, White RWG. Comparison of the influence of temperature on the egg development and growth of Daphnia longispina O.F. Müller (Crustacea: Cladocera) from two habitats in southern England. Oecologia. 1975;20:157–65.

    CAS  Article  PubMed  Google Scholar 

  19. 19.

    Taipale SJ, Kainz MJ, Brett MT. Diet-switching experiments show rapid accumulation and preferential retention of highly unsaturated fatty acids in Daphnia. Oikos. 2011;120:1674–82.

    Article  Google Scholar 

  20. 20.

    Brewer CM. Mating behaviours of Daphnia pulicaria, a cyclic parthenogen: comparisons with copepods. Philos Trans R Soc London Ser B Biol Sci. 1998;353:805–15.

    Article  Google Scholar 

  21. 21.

    Euent S, Menzel R, Steinberg CEW. Gender-specific lifespan modulation in Daphnia magna by a dissolved humic substances preparation. Ann Environ Sci. 2008;2.

  22. 22.

    MacArthur JW, Baillie WHT. Sex differences in mortality and metabolic activity in Daphnia magna. Science (80- ). 1926;64:229 LP – 230. doi:

    CAS  Article  Google Scholar 

  23. 23.

    Schwarzenberger A, Christjani M, Wacker A. Longevity of Daphnia and the attenuation of stress responses by melatonin. BMC Physiol. 2014;14:8.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  24. 24.

    Huylmans AK, López Ezquerra A, Parsch J, Cordellier M. De novo Transcriptome assembly and sex-biased gene expression in the cyclical Parthenogenetic Daphnia galeata. Genome Biol Evol. 2016;8:3120–39.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  25. 25.

    Raborn RT, Spitze K, Brendel VP, Lynch M. Promoter architecture and sex-specific gene expression in Daphnia pulex. Genetics. 2016;204:593–612.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  26. 26.

    Molinier C, Reisser CMO, Fields P, Ségard A, Galimov Y, Haag CR. Identification of General Patterns of Sex-Biased Expression in Daphnia, a Genus with Environmental Sex Determination. G3 (Bethesda). 2018;8:1523–33.

    CAS  Article  Google Scholar 

  27. 27.

    Norouzitallab P, Baruah K, Vanrompay D, Bossier P. Can epigenetics translate environmental cues into phenotypes? Sci Total Environ. 2019;647:1281–93.

    CAS  Article  PubMed  Google Scholar 

  28. 28.

    Kilvitis HJ, Hanson H, Schrey AW, Martin LB. Epigenetic potential as a mechanism of phenotypic plasticity in vertebrate range expansions. Integr Comp Biol. 2017;57:385–95.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  29. 29.

    Asselman J, De Coninck DIM, Pfrender ME, De Schamphelaere KAC. Gene body methylation patterns in Daphnia are associated with gene family size. Genome Biol Evol. 2016;8:1185–96.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  30. 30.

    Kvist J, Gonçalves Athanàsio C, Shams Solari O, Brown JB, Colbourne JK, Pfrender ME, et al. Pattern of DNA methylation in Daphnia: evolutionary perspective. Genome Biol Evol. 2018;10:1988–2007.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  31. 31.

    Robichaud NF, Sassine J, Beaton MJ, Lloyd VK. The epigenetic repertoire of Daphnia magna includes modified histones. Genet Res Int. 2012;2012:174860.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  32. 32.

    Gómez R, Van Damme K, Gosálvez J, Morán ES, Colbourne JK. Male meiosis in Crustacea: synapsis, recombination, epigenetics and fertility in Daphnia magna. Chromosoma. 2016;125:769–87.

    Article  PubMed  Google Scholar 

  33. 33.

    Asselman J, De Coninck DIM, Beert E, Janssen CR, Orsini L, Pfrender ME, et al. Bisulfite sequencing with Daphnia highlights a role for epigenetics in regulating stress response to Microcystis through preferential differential methylation of serine and threonine amino acids. Environ Sci Technol. 2017;51:924–31.

    CAS  Article  PubMed  Google Scholar 

  34. 34.

    Ye Z, Xu S, Spitze K, Asselman J, Jiang X, Ackerman MS, et al. A New Reference Genome Assembly for the Microcrustacean <em>Daphnia pulex</em>. G3 Genes|Genomes|Genetics. 2017;7:1405–16

    CAS  Article  Google Scholar 

  35. 35.

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

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  36. 36.

    Camp AA, Haeba MH, LeBlanc GA. Complementary roles of photoperiod and temperature in environmental sex determination in <em>Daphnia</em> spp. J Exp Biol. 2019;222.

  37. 37.

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

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  38. 38.

    Colbourne JK, Hebert PDN. The systematics of north American Daphnia (Crustacea: Anomopoda): a molecular phylogenetic approach. Philos Trans R Soc London B Biol Sci. 1996;351:349–60

    CAS  Article  Google Scholar 

  39. 39.

    Price DC, Egizi A, Fonseca DM. The ubiquity and ancestry of insect doublesex. Sci Rep. 2015;5:13068.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  40. 40.

    Dong X, Weng Z. The correlation between histone modifications and gene expression. Epigenomics. 2013;5:113–6.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  41. 41.

    Stillman B. Histone modifications: insights into their influence on gene expression. Cell. 2018;175:6–9.

    CAS  Article  PubMed  Google Scholar 

  42. 42.

    Ha M, Ng D, Li W-H, Chen ZJ. Coordinated histone modifications are associated with gene expression variation within and between species. Genome Res. 2011.

    CAS  Article  Google Scholar 

  43. 43.

    Akkers RC, van Heeringen SJ, Jacobi UG, Janssen-Megens EM, Françoijs K-J, Stunnenberg HG, et al. A hierarchy of H3K4me3 and H3K27me3 Acquisition in Spatial Gene Regulation in Xenopus embryos. Dev Cell. 2009;17:425–34.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  44. 44.

    Guenther MG, Levine SS, Boyer LA, Jaenisch R, Young RA. A chromatin landmark and transcription initiation at Most promoters in human cells. Cell. 2007;130:77–88.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  45. 45.

    Ferrari KJ, Scelfo A, Jammula S, Cuomo A, Barozzi I, Stützer A, et al. Polycomb-dependent H3K27me1 and H3K27me2 regulate active transcription and enhancer Fidelity. Mol Cell. 2014;53:49–62.

    CAS  Article  PubMed  Google Scholar 

  46. 46.

    Lesch BJ, Page DC. Poised chromatin in the mammalian germ line. Development. 2014;141:3619–26.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  47. 47.

    Bernhart SH, Kretzmer H, Holdt LM, Jühling F, Ammerpohl O, Bergmann AK, et al. Changes of bivalent chromatin coincide with increased expression of developmental genes in cancer. Sci Rep. 2016;6:37393.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  48. 48.

    Voigt P, Tee W-W, Reinberg D. A double take on bivalent promoters. Genes Dev. 2013;27:1318–38.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  49. 49.

    Lovén J, Orlando DA, Sigova AA, Lin CY, Rahl PB, Burge CB, et al. Revisiting global gene expression analysis. Cell. 2012;151:476–82.

    Article  PubMed  PubMed Central  Google Scholar 

  50. 50.

    Nie Z, Hu G, Wei G, Cui K, Yamane A, Resch W, et al. c-Myc Is a Universal Amplifier of Expressed Genes in Lymphocytes and Embryonic Stem Cells. Cell. 2012;151:68–79.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  51. 51.

    Pietrzak B, Grzesiuk M, Dorosz J, Mikulski A. When males outlive females: sex-specific effects of temperature on lifespan in a cyclic parthenogen. Ecol Evol. 2018;8:9880–8.

    Article  PubMed  PubMed Central  Google Scholar 

  52. 52.

    Kim YE, Hipp MS, Bracher A, Hayer-Hartl M, Ulrich HF. Molecular chaperone functions in protein folding and Proteostasis. Annu Rev Biochem. 2013;82:323–55.

    CAS  Article  PubMed  Google Scholar 

  53. 53.

    Calderwood SK, Murshid A, Prince T. The shock of aging: molecular chaperones and the heat shock response in longevity and aging a mini-review. Gerontology. 2009;55:550–8.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  54. 54.

    Wallin RPA, Lundqvist A, Moré SH, von Bonin A, Kiessling R, Ljunggren H-G. Heat-shock proteins as activators of the innate immune system. Trends Immunol. 2002;23:130–5.

    CAS  Article  PubMed  Google Scholar 

  55. 55.

    Mikulski A, Bernatowicz P, Grzesiuk M, Kloc M, Pijanowska J. Differential levels of stress proteins (HSPs) in male and female Daphnia magna in response to thermal stress: a consequence of sex-related behavioral differences? J Chem Ecol. 2011;37:670–6.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  56. 56.

    Schumpert C, Handy I, Dudycha JL, Patel RC. Relationship between heat shock protein 70 expression and life span in Daphnia. Mech Ageing Dev. 2014;139:1–10.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  57. 57.

    Wacker A, Martin-Creuzburg D. Allocation of essential lipids in Daphnia magna during exposure to poor food quality. Funct Ecol. 2007;21:738–47.

    Article  Google Scholar 

  58. 58.

    Abrusán G, Fink P, Lampert W. Biochemical limitation of resting egg production in Daphnia. Limnol Oceanogr. 2007;52:1724–8.

    Article  Google Scholar 

  59. 59.

    Yampolsky LY, Zeng E, Lopez J, Williams PJ, Dick KB, Colbourne JK, et al. Functional genomics of acclimation and adaptation in response to thermal stress in Daphnia. BMC Genomics. 2014;15:859.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  60. 60.

    Kilham SS, Kreeger DA, Lynn SG, Goulden CE, Herrera L. COMBO: a defined freshwater culture medium for algae and zooplankton. Hydrobiologia. 1998;377.

    CAS  Article  Google Scholar 

  61. 61.

    Athanasio CG, Chipman JK, Viant MR, Mirbahai L. Optimisation of DNA extraction from the crustacean Daphnia. PeerJ. 2016;4:e2004.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  62. 62.

    Athanasio CG, Sommer U, Viant MR, Chipman JK, Mirbahai L. Use of 5-azacytidine in a proof-of-concept study to evaluate the impact of pre-natal and post-natal exposures, as well as within generation persistent DNA methylation changes in Daphnia. Ecotoxicology. 2018;27:556–68.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  63. 63.

    Martin M. Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnet J. 2011;17:10

    Article  Google Scholar 

  64. 64.

    Pedersen BS, Eyring K, De S, Yang I V., Schwartz DA. Fast and accurate alignment of long bisulfite-seq reads. arXiv Prepr arXiv14011129. 2014.

  65. 65.

    Li H, Durbin R. Fast and accurate short read alignment with burrows–wheeler transform. Bioinformatics. 2009;25:1754–60.

    CAS  Article  Google Scholar 

  66. 66.

    Li B, Dewey CN. RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinformatics. 2011;12:323.

    CAS  Article  Google Scholar 

  67. 67.

    Dobin A, Davis CA, Schlesinger F, Drenkow J, Zaleski C, Jha S, et al. STAR: ultrafast universal RNA-seq aligner. Bioinformatics. 2013;29:15–21

    CAS  PubMed  PubMed Central  Google Scholar 

  68. 68.

    Gilbert DG. Longest protein, longest transcript or most expression, for accurate gene reconstruction of transcriptomes? bioRxiv. 2019. doi:

  69. 69.

    Leng N, Dawson JA, Thomson JA, Ruotti V, Rissman AI, Smits BMG, et al. EBSeq: an empirical Bayes hierarchical model for inference in RNA-seq experiments. Bioinformatics. 2013;29:1035–43

    CAS  Article  Google Scholar 

  70. 70.

    Sacomoto GAT, Kielbassa J, Chikhi R, Uricaru R, Antoniou P, Sagot M-F, et al. KIS SPLICE: de-novo calling alternative splicing events from RNA-seq data. BMC Bioinformatics. 2012;13:1–12.

    Article  Google Scholar 

  71. 71.

    Benoit-Pilven C, Marchet C, Chautard E, Lima L, Lambert M-P, Sacomoto G, et al. Annotation and differential analysis of alternative splicing using de novo assembly of RNAseq data bioRxiv 2016. doi:

  72. 72.

    Akalin A, Kormaksson M, Li S, Garrett-Bakelman FE, Figueroa ME, Melnick A, et al. methylKit: a comprehensive R package for the analysis of genome-wide DNA methylation profiles. Genome Biol. 2012;13:R87.

    Article  PubMed  PubMed Central  Google Scholar 

  73. 73.

    Wang H-Q, Tuominen LK, Tsai C-J. SLIM: a sliding linear model for estimating the proportion of true null hypotheses in datasets with dependence structures. Bioinformatics. 2011;27:225–31.

    CAS  Article  PubMed  Google Scholar 

  74. 74.

    Zhang Y, Liu T, Meyer CA, Eeckhoute J, Johnson DS, Bernstein BE, et al. Model-based analysis of ChIP-Seq (MACS). Genome Biol. 2008;9:R137.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  75. 75.

    Stark R, Brown G. DiffBind: differential binding analysis of ChIP-Seq peak data. R Packag version. 2011;100:3–4.

    Google Scholar 

  76. 76.

    Fabregat A, Sidiropoulos K, Garapati P, Gillespie M, Hausmann K, Haw R, et al. The Reactome pathway knowledgebase. Nucleic Acids Res. 2016;44:D481–7.

    CAS  Article  Google Scholar 

  77. 77.

    Yu G, Wang L-G, Han Y, He Q-Y. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS. 2012;16:284–7.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  78. 78.

    Yu G, He Q-Y. ReactomePA: an R/bioconductor package for reactome pathway analysis and visualization. Mol BioSyst. 2016;12:477–9.

    CAS  Article  PubMed  Google Scholar 

Download references


We thank the Environmental Omics Facility, University of Birmingham, and the Centre for Genomics and Bioinformatics, Indiana University for sequencing the RNA-seq and whole genome bisulfite sequencing libraries.


This research was supported by the Wellcome Trust ISSF funding, joint University of Birmingham and Brazil [Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (CAPES/BEX8572-12-7)] scholarship program, and the University of Birmingham Fellowship Award. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Author information




JK performed the bulk of the computation analysis and manuscript preparation, JBB contributed to data analysis. CGA performed all the experimental work. MEP provided the Daphnia pulex EB clones and contributed to manuscript preparation. JKC contributed to the manuscript writing. LM developed the experimental design, contributed to the experimental work and manuscript preparation. All authors read and approved the final manuscript.

Corresponding authors

Correspondence to Jouni Kvist or Leda Mirbahai.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

JBB is an editorial member for BMC Genomics journal.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary information

Additional file 1: Table S1.

Differentially expressed or regulated events between male and female Daphnia pulex. Results are in separate tabs; A) Differentially expressed transcripts, B) Differentially expressed genes, C) Differentially expressed splice variants, D) Differentially methylated CpGs, E) Differentially modified histone H3 trimethylated at lysine 4 (H3K4me3), F) Differentially modified histone H3 trimethylated at lysine 27 (H3K27me3).

Additional file 2: Table S2

. Enriched Reactome pathways in differentially expressed or regulated events between male and female Daphnia pulex. A) Enrichment of differentially expressed transcripts between male and female D. pulex, B) Enrichment of differentially expressed genes, C) Enrichment of differentially expressed transcripts not contained in differentially expressed genes, D) Enrichment of differentially regulated spice variants between male and female, E) Enrichment of genes with differentially methylated CpGs, F) Enrichment of genes with differentially regulated histone H3 trimethylation at lysine 4 (H3K4me3), G) Enrichment of genes with differentially regulated histone H3 trimethylation at lysine 4 (H3K4me3), H) Enrichment analysis of genes with co-occurring modifications (CpG methylation, H3K4me3 and H3K27me3), I) Enrichment analysis of main clusters of genes with co-varying expression and modifications from heatmap in Fig. 7.d.

Additional file 3: Fig. S1

. Comparison of sex-biased expression between D. magna and D. pulex. The orthologs were identified with blastp and limited to single gene matches with e-values <1e− 20. Sex-biased genes in D. magna are based on [26]. Genes marked in red are significantly different (PPEE< 0.05) in D. pulex, based on this study. Doublesex genes (DapmaDsx1 and DapmaDsx2) are highlighted in purple.

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (, which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Kvist, J., Athanàsio, C.G., Pfrender, M.E. et al. A comprehensive epigenomic analysis of phenotypically distinguishable, genetically identical female and male Daphnia pulex. BMC Genomics 21, 17 (2020).

Download citation


  • Epigenetics
  • Gene expression
  • Evolution
  • Non-conventional model organisms