Epigenetic mechanisms such as DNA methylation are a key component of dosage compensation on sex chromosomes and have been proposed as an important source of phenotypic variation inﬂuencing plasticity and adaptive evolutionary pro- cesses, yet little is known about the role of DNA methylation in an ecological or evolutionary context in vertebrates. The threespine stickleback (Gasterosteus aculeatus) is an ecological and evolutionary model system that has been used to study mechanisms involved in the evolution of adaptive phenotypes in novel environments as well as the evolution heteromorphic sex chromosomes and dosage compensation in vertebrates. Using whole genome bisulﬁte sequencing, we compared genome-wide DNA methylation patterns between threespine stickleback males and females and between stickleback reared at different environmental salinities. Apparent hypermethylation of the younger evolutionary stratum of the stick- leback X chromosome in females relative to males suggests a potential role of DNA methylation in the evolution of hetero- morphic sex chromosomes. We also demonstrate that rearing salinity has genome-wide effects on DNA methylation levels, which has the potential to lead to the accumulation of epigenetic variation between natural populations in different environments. Key words: epigenetics, Gasterosteus aculeatus, bisulﬁte sequencing, WGBS, BS-Seq. Introduction environment (Boyko et al. 2010; Kucharski et al. 2008; Epigenetic variation has the potential to impact ecological Cooney et al. 2002). VariationinDNA methylationlevels and evolutionary processes, and thus affect species distri- have therefore been hypothesized to play a key role in butions and evolutionary trajectories (Bossdorf et al. mediating phenotypic responses to environmental change 2008; Flores et al. 2013; Jablonka and Raz 2009; (Bossdorf et al. 2008; Hofmann 2017; Flores et al. 2013), Varriale 2014; Franks and Hoffmann 2012). Currently, and may represent a dynamic source of heritable variation one of the best-studied mechanisms underlying epige- that can respond to changes in the environment and in- netic variation is DNA methylation, a heritable epigenetic ﬂuence phenotypic variation over multiple time-scales modiﬁcation in which a methyl group is added to position (Richards 2006). 5 of the pyrimidine ring on a cytosine (5mC), most com- In addition to its potential role in regulating gene expres- monly found on cytosine-phosphate-guanine (CpG) dinu- sion in response to environmental change, DNA methylation cleotides in vertebrates (Heard and Martienssen 2014). is also critical in regulating gene expression in dosage com- Changes in DNA methylation can have profound effects pensation systems that have evolved to minimize the unequal on chromatinstructure, which caninturnalter gene ex- expression of genes on heteromorphic sex chromosomes pression (Klose and Bird 2006; Jaenisch and Bird 2003). (Graves 2016). In older XY sex chromosome systems, such The addition or removal of these methyl groups can be as those found in most mammalian species, DNA methylation dynamically regulated in response to changes in the is involved in the global silencing of one of the two X The Author(s) 2018. Published by Oxford University Press on behalf of the Society for Molecular Biology and Evolution. This is an Open Access article distributed under the terms of the Creative Commons Attribution Non-Commercial License (http://creativecommons.org/licenses/by-nc/4.0/), which permits non-commercial re-use, distribution, and reproduction in any medium, provided the original work is properly cited. For commercial re-use, please contact firstname.lastname@example.org Genome Biol. Evol. 10(3):775–785. doi:10.1093/gbe/evy034 Advance Access publication February 6, 2018 775 Downloaded from https://academic.oup.com/gbe/article-abstract/10/3/775/4840697 by Ed 'DeepDyve' Gillespie user on 16 March 2018 Metzger and Schulte GBE chromosome in females (Graves 2016). In ZW sex chromo- Materials and Methods some systems (which have a female-speciﬁc W chromosome) Fish Collection such as those in birds and some reptiles and ﬁshes, DNA All animal experimentation was conducted according to methylation is involved in gene-speciﬁc dosage compen- University of British Columbia approved animal care protocols sation via the activation or suppression of particular dos- (A10-0285 and A11-0372). Adult threespine stickleback (G. age sensitive genes (Graves 2016). While epigenetic aculeatus) of the fully plated “marine” ecotype were collected silencing has been well-established as a mechanism in- at the beginning of their natural spawning season in May volved in dosage compensation of older heteromorphic 2013 from Oyster Lagoon, British Columbia in Canada sex chromosome systems, patterns of DNA methylation (GPS: 49.6121, -124.0314). Fish were separated into six in young sex chromosome systems are less well-under- 110-l glass tanks (20 ﬁsh per tank) and acclimated to 21 stood, but have been hypothesized to play a key role sex ppt salt water (dechlorinated Vancouver municipal tap water chromosome evolution (Gorelick 2003). supplemented with Instant Ocean Sea Salt), 18 Cand The threespine stickleback (Gasterosteus aculeatus)has 14:10 h light:dark photoperiod. These conditions mimic the been extensively used to investigate the genetic basis of adap- natural environmental conditions at the collection location at tive evolution to novel environments (Jones et al. 2012a,b). the time of collection. Fish were fed daily to satiation with Following the last glaciation, ancestral marine populations of Hakari Bio-Pure frozen Mysis Shrimp and were acclimated to stickleback colonized and adapted to newly available fresh- laboratory conditions for four weeks. water habitats in the north-temperate zone (Bell and Foster 1994). Adaptation to these novel environments drove the rapid parallel evolution of divergent phenotypes including changes in body shape, armor plate number, gene expression Fertilization and Rearing Procedure levels, and gene expression plasticity (Jones et al. 2012a; Gibbons et al. 2017; Morris et al. 2014; McCairns and To determine the impact of salinity on fertilization and hatch- Bernatchez 2009; Ishikawa et al. 2017). Several studies have ing, eggs were collected from gravid female stickleback and used reduced representation approaches to characterize var- immediately divided into six different petri dishes containing iation in DNA methylation patterns between stickleback that 5 ml of 2, 7, 14, 21, 28, or 35 ppt saltwater. Testes were vary in their lateral plate morphology, and have suggested collected from males displaying sexually mature characteristics that variation in DNA methylation patterns may contribute and individually macerated in a 1.75 ml microcentrifuge tube to the phenotypic divergence observed between marine and containing 300 ml Ginzberg’s ﬁsh Ringer’s solution. Eggs from freshwater populations (Smith et al. 2015; Trucchi et al. 2016; a single clutch were fertilized with sperm solution from a sin- Artemov et al. 2017). Threespine stickleback also have a rel- gle male across all salinities (50 ml of sperm solution for each atively young XY sex chromosome pair that has evolved since petri dish at each different salinity). Following fertilization, an the species ﬁrst arose 13–16 Ma (Bell et al. 2009; Kawahara additional 10 ml of water at the appropriate salinity was et al. 2009; Ross et al. 2009), and this species has become a added to each petri dish. This process was repeated ten times powerful model system to explore the evolution of hetero- creating a total of ten different families, each fertilized at all morphic sex chromosomes and dosage compensation mech- salinities. Petri dishes were partially covered to prevent water anisms (Schultheiß et al. 2015; White et al. 2015). Thus, the loss via evaporation and to allow for surface gas exchange. threespine stickleback is an ideal model in which to investigate Eggs were monitored twice daily during which time any the complementary roles of DNA methylation in both envi- unfertilized eggs were removed and 10 ml of water was ronmental adaptation and the evolution of sex chromosome changed with sterilized water of the appropriate salinity to systems. prevent mold growth. Percent fertilization and percent hatch In this study, we present the ﬁrst high-resolution analysis of were recorded. Percent hatch is recorded as the proportion of DNA methylation patterns in the stickleback genome using fertilized embryos that hatched. The effect of salinity on fer- whole genome bisulﬁte sequencing (WGBS). This approach tilization and hatching success was analyzed using a logistic allowed us to characterize prominent features in the DNA regression (Warton and Hui 2011) in the R v3.3.2 base stats methylation landscape of stickleback, including variation in package. Tukey post hoc analysis was performed using the DNA methylation patterns between males and females along glht() function in the multicomp v1.4-6 R package. the entire stickleback sex chromosome, which provides insight After all ﬁsh in a petri dish had hatched and the yolks had into the relationship between epigenetic mechanisms and sex been absorbed (15 days post fertilization), larvae were trans- chromosome evolution. By rearing putatively ancestral marine ferred to hanging net boxes (Aquaclear) in 110 L glass aquaria stickleback at both low and high salinity, we also describe the containing water at the fertilization salinity. Sponge ﬁlters effects of environmental salinity on genomic DNA methyla- were used for ﬁltration and aeration. Each family was kept tion patterns, and highlight potential salinity responsive genes separate throughout the experiment. At one-month post that may be differentially regulated by DNA methylation. hatch whole animals were snap frozen in liquid nitrogen. 776 Genome Biol. Evol. 10(3):775–785 doi:10.1093/gbe/evy034 Advance Access publication February 6, 2018 Downloaded from https://academic.oup.com/gbe/article-abstract/10/3/775/4840697 by Ed 'DeepDyve' Gillespie user on 16 March 2018 DNA Methylation Landscape of Stickleback GBE WGBS 10 kb windows of the genome. All ﬁgures were generated in R. Genomic DNA was isolated from one-month old whole ﬁsh To obtain nearest neighboring gene annotations, the gene samples from the 2 ppt and the 21 ppt salinity treatments coordinates in the annotation ﬁle (.gtf) provided by Glazer using a Qiagen DNeasy Blood and Tissue Kit following the et al. (2015) in the dryad digital repository (which contains manufacturer’s recommended protocol for RNA-free geno- gene coordinates that correspond to the stickleback genome mic DNA using RNAase A. The sex of each sample was iden- available in Ensembl) were converted to the coordinates in the tiﬁed by PCR analysis using primers designed to idh, gasm6, updated assembly using the R script convertCoordinate.R that and stn190 following previously described methods (Toli et al. is provided by the authors of the revised assembly. Distances 2016). Genomic DNA samples from three males and three to nearest neighboring genes were calculated using the females from 2 ppt and 21 ppt (twelve samples total) were annotatePeakInBatch() function in the R package sent to the McGill University and Genome Quebec Innovation ChIPpeakAnno v3.6.5. Center for DNA quality assessment, library preparation, bisul- CpG islands for the revised stickleback genome assembly ﬁte treatment, and 150 base pair paired-end sequencing us- were calculated using python scripts (https://github.com/ ing an Illumina HiSeqX. The 12 samples were split evenly lucasnell/TaJoCGI) that apply an algorithm based on the across three sequencing lanes (4 samples/lane) such that methods described by Takai and Jones (2002). The observed one male and female sample from each of the salinity treat- distribution of DMCs was compared with the distribution of ments were represented on each sequencing lane. Average CpGs across the genome using a Chi-square test. sequencing library size was 102,011,555 reads (6 13,147,873 SD). Candidate Gene Analysis WGBS Data Analysis Previous RNA-seq studies have identiﬁed many genes that respond to changes in salinity in stickleback (2,771 in gill tis- Reads were mapped to a revised assembly of the stickleback sue, Gibbons et al. 2017 and 1,844 in kidney tissue, Wang genome (Glazer et al. 2015) obtained from the Dryad Digital et al. 2014). To determine whether DNA methylation could be Repository (http://datadryad.org/resource/doi: 10.5061/dryad. involved in the differential regulation of these candidate q018v) and DNA methylation levels were calculated using the genes we compiled a list of salinity responsive genes from bisulﬁte sequencing plugin v1.2 in CLC Genomics Workbench previous studies (4,615 candidate genes) and compared v10.0. Average mapping efﬁciency was 89.5% (6 1% SD). them to genes within 2 kb of DMCs in stickleback reared at DNA methylation data were exported from CLC and analyzed different salinities. using in R v3.3.2 using methylKit package v3.5 (Akalin et al. Similarly, several studies have also characterized sex-biased 2012) following previously recommended guidelines for bisul- gene expression patterns in stickleback. We therefore com- ﬁte sequence analysis (Ziller et al. 2015; Wreczycka et al. piled a list of 2,282 genes that display sex-biased gene expres- 2017). Sequenced CpG loci were ﬁltered so that only sites sion patterns in brain (1,255 genes, Metzger and Schulte with at least 10 reads in each of the 12 samples were retained. 2016) and liver tissue (1,268 genes, Leder et al. 2010), and Sites that were in the 99.9th percentile of coverage were also compared these to genes within 2 kb of DMCs that were removed from the analysis to account for potential PCR bias. identiﬁed between male and female stickleback. Hierarchical cluster analysis was conducted using Ward’s We also compared genes within 2 kb of DMCs between method and was implemented using the clusterSamples() individuals reared a different salinities to genes associated function. Pairwise comparisons between groups were per- with single nucleotide polymorphisms (SNPs) under positive formed using a logistic regression model with a correction selection in threespine stickleback from freshwater environ- for overdispersion using the calculateDiffMeth() function ments compared with marine environments (Jones et al. followed by a Chi-square test to identify signiﬁcantly differ- 2012a). entially methylated cytosines (DMCs) between groups. The P- values for DMCs were false discovery rate (FDR) corrected using the sliding linear model method (SLIM) with a maximum Results and Discussion q-value threshold of 0.05 and a minimum change in percent Characterization of the Stickleback Methylome methylation of 10%. For the comparison between males and females, salinity and family were included as covariates. For We performed WGBS on ﬁsh from a marine population of the comparison between salinity rearing treatments, sex and threespine stickleback reared from fertilization to the age of family were included as covariates. For comparisons between 1 month at a salinity of either 2 ppt or 21 ppt. These salinities families, sex and salinity were included as covariates. To cal- represent the widest range that still allows good fertilization culate mean methylation levels across 10 kilobase (kb) geno- and hatching success in this population (supplementary ﬁg. mic regions, the tileMethylCounts() function in methylKit was S1, Supplementary Material online). In this study, we utilized a used to calculate DNA methylation values across sequential balanced design with WGBS performed on one male and one Genome Biol. Evol. 10(3):775–785 doi:10.1093/gbe/evy034 Advance Access publication February 6, 2018 777 Downloaded from https://academic.oup.com/gbe/article-abstract/10/3/775/4840697 by Ed 'DeepDyve' Gillespie user on 16 March 2018 Metzger and Schulte GBE female from each of three families and each of the salinity dimorphic methylation may be highly variable among teleosts, rearing treatments. This design was chosen because genetic consistent with the wide range of sex-determining mecha- variation has been shown to have substantial effects on the nisms in this group (Devlin and Nagahama 2002). divergence of DNA methylation patterns among individuals in The distribution of DMCs across genomic features (e.g., both plants and animals (Gertz et al. 2011; McRae et al. 2014; promoters, exons, CpG islands) did not differ from the relative Vidalis et al. 2016). Consistent with this observation, we proportions of these features in the genome (supplementary detected strong effects of family on DNA methylation (sup- table S2, Supplementary Material online). No DMCs were plementary ﬁgs. S2 and S3, Supplementary Material online). identiﬁed between males and females in the mitochondrial However, family-level variation in DNA methylation levels genome. The majority of DMCs (90%; 16,626 DMCs) be- could also be indicative of transgenerational environmental tween males and females showed a bias towards higher or maternal effects (Jablonka and Raz 2009). methylation in females relative to males suggesting that fe- CpG loci in the stickleback genome had an average meth- male stickleback genome is hypermethylated relative to male ylation level of 70.3%, which is consistent with whole ge- stickleback genome. nome assessments of methylation in other ﬁsh species (Feng The most striking pattern in these data is the apparent et al. 2010; Shao et al. 2014; Zemach et al. 2010). However, hypermethylation of chr19 (the threespine stickleback sex there were several hypomethylated regions (<40% methyla- chromosome) in females relative to males, which is where tion) that are indicative of DNA methylation canyons or valleys 65% of the putative DMCs identiﬁed between males and (Jeong et al. 2014; Xie et al. 2013) with the most prominent females are located (12,112 DMCs). Chr19 also had the high- of these located on chromosomes 4, 10, 11, 12, and 16 est proportion of DMCs relative to the number of CpG sites (ﬁg. 1, supplementary ﬁg. S4, Supplementary Material online). on the chromosome (5%) compared with the rest of the While the factors that determine the size of DNA methylation autosomes where 0.07% of CpG loci were differentially canyons remains unknown, larger hypomethylated canyons methylated between the sexes (ﬁg. 2). such as those described here in stickleback have been shown Three distinct regions (strata) have been characterized on to be under strong transcriptional suppression due to in- chr19 based on the extent of divergence in these regions creased abundance of repressive histone H3 lysine 27 meth- between the X and Y chromosome (Ross and Peichel 2008; ylation that interacts with hypomethylated DNA (Nakamura White et al. 2015), and two of these strata no longer recom- et al. 2014). This mechanisms of transcriptional repression is bine between the X and Y: Stratum two (the younger evolu- thought to maintain these regions in a “poised” transcrip- tionary stratum located between 2.5 Mb and 12 Mb), and tional state to allow rapid activation of gene transcription in stratum one (the older evolutionary stratum located from these regions at speciﬁc times during embryonic development 12 Mb to the end of the chromosome). There is also a pseu- (Nakamura et al. 2014), but whether these hypomethylated doautosomal region (PAR) located in the ﬁrst 2.5 Mb of chr19 canyons play a functional role in adults is unknown. However, that is thought to still recombine between the X and Y chromo- genes that are essential for proper development typically somes (White et al. 2015). To assess whether these evolutionary dominate these regions (Jeong et al. 2014; Nakamura et al. strata are also associated with sex-speciﬁc DNA methylation, we 2014; Xie et al. 2013). Consistent with this pattern, genes divided chr19 into 10 kb consecutive nonoverlapping bins and located in the hypomethylated canyons in the stickleback ge- calculated the frequency of CpG loci that were putatively iden- nome include protocadherins on chromosome 4 and homeo- tiﬁed as hypermethylated or hypomethylated in female stick- box genes on chromosomes 10, 11, 12, and 16 leback relative to male stickleback (ﬁg. 3). This analysis revealed (supplementary table S1, Supplementary Material online), clear patterns that correspond to the evolutionary strata of suggesting a conserved role of hypomethylated canyons chr19 for loci that were hypermethylated in females relative across vertebrates. to males, with the greatest apparent hypermethylation in stra- tum two and the least in the PAR (ﬁg. 3B). Sex-Biased DNA Methylation Patterns Identifying Differential Methylation on the Sex We identiﬁed a total of 18,564 DMCs between males and Chromosome females (ﬁg. 2, supplementary data set S1, Supplementary Material online). Although relatively few studies have exam- One of the challenges for unambiguously determining levels ined differential methylation patterns between males and of sex-speciﬁc DNA methylation on stickleback sex chromo- females at the whole genome level in ﬁshes, a study in tilapia somes is that there is currently no publically available se- detected a similar number of DMCs between males and quence for the Y chromosome and the published sequence females in muscle tissue (17,112 DMCs; Wan et al. 2016), for chr19 is predominately derived from X chromosome se- whereas a study of sex-speciﬁc differential methylation in quence. Because DNA methylation is detected as sequence zebraﬁsh brain detected only 914 DMCs (Chatterjee et al. differences between bisulﬁte-treated DNA and the reference 2016). These data suggest that the extent of sexually sequence at CpG sites, both divergence between the X and Y 778 Genome Biol. Evol. 10(3):775–785 doi:10.1093/gbe/evy034 Advance Access publication February 6, 2018 Downloaded from https://academic.oup.com/gbe/article-abstract/10/3/775/4840697 by Ed 'DeepDyve' Gillespie user on 16 March 2018 DNA Methylation Landscape of Stickleback GBE FIG.1.—Mean CpG methylation level across chromosomes 4, 10, 11, 12, and 16. Each point represents the mean methylation level across all twelve individuals for a single 10 kb window. The solid line represents the smoothed spline ﬁt to these data. Position along the x-axis represents the base position along the chromosome. The y-axis is the average DNA methylation level. Genome Biol. Evol. 10(3):775–785 doi:10.1093/gbe/evy034 Advance Access publication February 6, 2018 779 Downloaded from https://academic.oup.com/gbe/article-abstract/10/3/775/4840697 by Ed 'DeepDyve' Gillespie user on 16 March 2018 Metzger and Schulte GBE FIG.2.—Differentially methylated CpG loci between male and female stickleback. (A) Differentially methylated CpG (DMC) loci between male and female threespine stickleback. Each point represents an individual DMC. The y-axis indicates the percent difference in methylation between males and females. A positive value on the y-axis indicates a DMC that is hypermethylated in females relative to males. A negative value on the y-axis indicates a DMC that is hypomethylated in females relative to males. The x-axis indicates the position of the DMC in the stickleback genome. Chromosome boundaries are represented by vertical dashed lines. Only DMCs for which a change in methylation of >10% are presented. (B) The percentage of CpG loci on a given chromosome that were differentially methylated between male and female threespine stickleback. The light shading represents DMCs that are hyper- methylated and the dark shading represents DMCs that are hypomethylated in female stickleback compared with males. chromosome at CpG sites and differential methylation can cannot unambiguously determine whether differential meth- result in the same patterns in the sequence data. Thus, the ylation or X Y polymorphism is the cause of the apparent signal of differential methylation that we observe could be hypermethylation in females in stratum two. attributed to one of three possible mechanisms: 1) identical Stratum one is thought to be the most divergent region of sequences but differential methylation between the X and Y the sex chromosome, and thus we predicted that few reads chromosome, 2) differential methylation between X chromo- from the Y chromosome would map to the chr19 reference, somes in males and females, or 3) sequence divergence be- resulting in a ratio of 0.5 for the number of reads mapped in tween the X and Y chromosomes resulting in alteration of males and females. Again, the results from the read coverage CpG sites. analysis mostly matched this prediction (ﬁg. 3D), although To address this issue, we again divided chr19 in to 10 kb there were speciﬁc regions where the read count ratio was consecutive nonoverlapping segments and compared the close to one. These regions in stratum one may correspond to number of reads that mapped to chr19 in males and females regions that are thought to be under purifying selection to for each 10 kb segment. If chr19 reads map uniquely to the maintain dosage sensitive genes from being lost on the Y published X chromosome then we would expect half the chromosome (White et al. 2015). We also detected apparent number of reads to map to chr19 in males compared with hypermethylation in females in these regions (ﬁg. 3B and D). females. Given that the PAR is known to recombine between Thus, we cannot conclusively determine whether this appar- the X and Y chromosomes, suggesting low levels of diver- ent hypermethylation of chr19 in females is due to differential gence in this region, we predicted that sequencing reads de- methylation or from the accumulation of TG polymorphisms rived from both the X and Y chromosomes would map to the on the Y that are being interpreted as unmethylated loci. reference sequence, resulting in a ratio of one for the number While we are unable to unambiguously determine the ul- of reads mapped in males and females in the PAR. The results timate cause of the apparent DNA methylation differences from the read coverage analysis matched this prediction between males and females on chr19, whether the patterns (ﬁg. 3D). we observe are the result of sequence polymorphism be- For the younger, less diverged stratum (stratum two), we tween the X and Y chromosome that alters CpG sites, or predicted a read count ratio between 0.5 and 1 because se- are due to differential methylation of conserved sequences quence similarity between the X and Y chromosome would between males and females in chr19, the ultimate effect result in reads from both chromosomes mapping to the X would be differences in methylation between the sex chro- chromosome reference sequence. The results from the read mosomes. Thus, taken together, the patterns of putative dif- coverage analysis matched this prediction (ﬁg. 3D). Thus, we ferential methylation that we observe suggest that divergence 780 Genome Biol. Evol. 10(3):775–785 doi:10.1093/gbe/evy034 Advance Access publication February 6, 2018 Downloaded from https://academic.oup.com/gbe/article-abstract/10/3/775/4840697 by Ed 'DeepDyve' Gillespie user on 16 March 2018 DNA Methylation Landscape of Stickleback GBE in DNA methylation patterns between males and females on the stickleback sex chromosome are closely associated with PAR Stratum Two Stratum One the known evolutionary history of this chromosome. 60 DNA Methylation and Sex Chromosome Evolution Sex chromosome evolution from autosomes is thought to in- volve recombination suppression in sex determining regions, followed by the accumulation of deleterious mutations and the degeneration of the sex-speciﬁc (e.g., Y) chromosome (Graves 2016). Degradation of the sex-speciﬁc heteromorphic 0 500 1000 1500 2000 chr 19 (10kb) sex chromosome following recombination suppression has the potential to cause imbalances in gene expression. Many 0.20 taxa with heteromorphic sex chromosome pairs have evolved PAR Stratum Two Stratum One dosage compensation mechanisms to resolve this effect, but 0.15 the nature and extent of these dosage compensation mech- anisms varies greatly among taxa (Graves 2016). DNA meth- 0.10 ylation has been proposed as a key mechanism responsible for regulating every step of the evolution of sex chromosomes 0.05 from recombination suppression in the early stages of sex chromosome evolution to dosage compensation in more de- rived sex chromosome systems (Gorelick 2003); however, 0.00 0 500 1000 1500 2000 there has been little empirical evidence to test this hypothesis. chr 19 (10kb) Taxa with relatively “young” heteromorphic sex chromo- somes, such as the threespine stickleback, provide an oppor- 0.20 PAR Stratum Two Stratum One tunity to investigate the potential role of DNA methylation in regulating sex chromosome recombination and the evolution 0.15 of dosage compensation mechanisms. In the following sec- tion, we discuss the apparent differential methylation be- 0.10 tween males and females on chr19 in the context of the different stages of sex chromosome evolution in stickleback. 0.05 DNA methylation promotes the formation of heterochro- matin (Mirouze et al. 2012; Yelina et al. 2015; Melamed- 0.00 Bessudo and Levy 2012), and it is thought to play a role in 0 500 1000 1500 2000 chr 19 (10kb) suppressing recombination of sex chromosomes in plants (Zhang et al. 2008). In stickleback we observed apparent hypermethylation of the younger evolutionary stratum in 2.5 PAR Stratum Two Stratum One females (stratum two), and relatively less differential methyl- 2.0 ation between males and females along the older evolution- ary stratum (stratum one) and the PAR. The apparent 1.5 hypermethylation of stratum two on the X chromosome in 1.0 females (hypomethylated in males) corresponds to the region hypothesized to have undergone the ﬁrst chromosomal inver- 0.5 sion during the evolution of the Y chromosome (Ross and Peichel 2008). The apparent differential methylation in 0.0 0 500 1000 1500 2000 chr 19 (10kb) FIG.3.—Differential methylation between sexes on chromosome 19 FIG.3.—Continued (chr19). (A) Mean DNA methylation levels for CpG loci along chromosome total number of DMCs in a 10 kb window relative to the number of CpG chr19. Each point represents the mean DNA methylation level in a 10 kb loci in that same 10 kb window. (D) Ratio of mapped read counts for males window for six individual stickleback that were either male (blue) or female relative to females along chr19. Each point represents the ratio of mean (red). Solid lines represent the smooth spline ﬁt for the DNA methylation counts for a 10 kb window in males compared with females. The solid levels in males (blue) and females (red). (B, C) Proportion of DMCs along black line is the smooth spline ﬁt. Vertical dashed lines represent the chr 19 that are hypermethylated (B) or hypomethylated (C) in female stick- boundaries between the three evolutionary strata on chr19: The pseu- leback compared with male stickleback. Values on the y-axis represent the doautosomal region (PAR), stratum two, and stratum one. Genome Biol. Evol. 10(3):775–785 doi:10.1093/gbe/evy034 Advance Access publication February 6, 2018 781 Downloaded from https://academic.oup.com/gbe/article-abstract/10/3/775/4840697 by Ed 'DeepDyve' Gillespie user on 16 March 2018 Male/Female read counts hypomethylated CpG/ total CpG hypermethylated CpG/ total CpG Percent methylation Metzger and Schulte GBE stratum two between males and females (hypomethylated in gene-speciﬁc regulation, which is not consistent with a gener- malesand hypermethylatedinfemales) couldhave playeda alized hypertranscription of genes in stratum one in females role in suppressing recombination between male and female (Schultheiß et al. 2015). Because the male to female coverage sex chromosomes and in establishing the boundaries in which ratio is similar in these localized regions in stratum one, this this inversion ﬁrst occurred. localized pattern may be more consistent with the potential The next stage in sex chromosome evolution, following preservation of dosage sensitive genes in these regions (White recombination suppression, is thought to be the accumulation et al. 2015); however, it is also possible that the differential of genetic variation and degradation in the nonrecombining methylation in these regions is caused by the accumulation region(s). If methylated cytosines on the female X chromo- of C to T polymorphisms which is less consistent with purifying some correspond to thymines on the male Y chromosome as selection acting in these regions (White et al. 2015). previously discussed, then it is possible that the accelerated Taken together, the apparent differential methylation be- mutation rate of methylated cytosines, which can be deami- tween male and female stickleback described in this study is nated to become thymines (Coulondre et al. 1978; Shen et al. consistent with the proposed role of DNA methylation in the 1994), could play an important role in the divergence be- evolution of sex chromosomes (Gorelick 2003). Thus differen- tween X and Y chromosomes. Alternatively, instead of being tial DNA methylation could be playing a role in the evolution a direct result of C to T polymorphisms, the observed increase of this “young” heteromorphic sex chromosome system, ei- in hypermethylated loci in females could be closely linked to ther through inﬂuencing patterns of recombination or poten- the accumulation of genetic polymorphisms. The frequency at tially through mediating the early stages of the development which SNPs occur in genomes has been shown to be higher of dosage compensation. near methylated CpG loci (Qu et al. 2012). The CGCG motif has been identiﬁed as a candidate cis-element associated with DNA Methylation and Sex-Biased Gene Expression this observation and is enriched in hypomethylated regions (Qu et al. 2012). Therefore, individuals with higher DNA To determine whether the apparent variation in DNA meth- methylation levels at particular loci would be predicted to ylation patterns between males and females could be have a higher degree of sequence divergence near those inﬂuencing previously described sex-biased gene expression loci relative to individuals with lower methylation levels. The patterns in stickleback, we compared the list of genes near PAR had the highest frequency of the CGCG motif (1.98/1 kb) DMCs between males and females (supplementary data set whereas stratum one and stratum two had a lower frequency S2, Supplementary Material online) to previously identiﬁed of the CGCG motif (1.22/1 kb and 1.18/1 kb respectively). genes that exhibit sex-biased gene expression patterns These data suggest that DNA methylation may also be asso- (Metzger and Schulte 2016; Leder et al. 2010). Of the ciated with the accumulation of genetic variation between 2,282 genes that have been shown to exhibit sex biased the nonrecombining regions of the male and female sex gene expression patterns from these studies, 490 overlapped chromosomes. with genes near DMCs in our study of which 269 are on chr19 We next explored whether the apparent differential meth- (supplementary data set S3, Supplementary Material online) ylation patterns between males and females on chr19 are including genes located in the region considered to be tightly consistent with the regulation of dosage sensitive genes. In linked to sex determination in stickleback such as sema4ba stickleback, there are two conﬂicting hypotheses regarding and idh2 (Peichel et al. 2004). This pattern is consistent with the existence of a dosage compensation system. One hypoth- differential methylation between males and females playing a esis is that there is locally conﬁned partial dosage compensa- role in regulating sex-biased patterns of gene expression. tion in stratum one in males that is also associated with a hypertranscription of genes in stratum one in females Effects of Environmental Salinity on DNA Methylation (Schultheiß et al. 2015). The second hypothesis suggests that dosage compensation has not evolved in the stickleback VariationinDNA methylationhas also been suggestedtobe but that there is purifying selection to maintain dosage sensi- an important component of an organism’s response to envi- tive genes in stratum one of the Y chromosome (White et al. ronmental change (Bossdorf et al. 2008; Hofmann 2017; 2015). The differential methylation patterns between sexes Flores et al. 2013). Changes in environmental salinity are along the X chromosome that we observe are not entirely known to cause substantial changes in gene expression in consistent with either of these prevailing hypotheses. We ob- many species of ﬁsh, including stickleback (Gibbons et al. served apparent hypermethylation in stratum two in females. 2017; Zhang et al. 2017; Wang et al. 2014). In order to ex- This might be expected to result in reduced transcription or plore whether changes in DNA methylation may be involved partial silencing of genes in this region, which has not been in environmental regulation of gene expression, we identiﬁed observed in stickleback (Schultheiß et al. 2015; White et al. differentially methylated loci in stickleback reared at two sal- 2015). The less extensive and highly localized pattern of hyper- inities (2 and 21 ppt), and we compared genes near DMCs methylation in stratum one that we observe is suggestive of identiﬁed in ﬁsh reared at different salinities to genes that 782 Genome Biol. Evol. 10(3):775–785 doi:10.1093/gbe/evy034 Advance Access publication February 6, 2018 Downloaded from https://academic.oup.com/gbe/article-abstract/10/3/775/4840697 by Ed 'DeepDyve' Gillespie user on 16 March 2018 DNA Methylation Landscape of Stickleback GBE B 0.04 hypomethylated hypermethyated 0.03 0.02 −20 0.01 −40 −60 0.00 FIG.4.—Differentially methylated CpG loci between stickleback reared at low and high salinity. (A) Differentially methylated CpG (DMC) loci between stickleback reared at a salinity of 2 ppt compared with 21 ppt. Each point represents an individual DMC. The y-axis indicates the percent difference in methylation between salinity rearing treatments. A positive value on the y-axis indicates a DMC this is hypermethylated in stickleback reared at 21 ppt relative 2 ppt. A negative value on the y-axis indicates a DMC that is hypomethylated in stickleback reared at 21 ppt relative to 2 ppt. The x-axis indicates the position of the DMC in the stickleback genome. Chromosome boundaries are represented by vertical dashed lines. Only DMCs for which a change in methylation of >10% are presented. (B) The percentage of CpG loci on a given chromosome that were detected as being differentially methylated between stickleback reared at 21 ppt and 2 ppt. The light shading represents DMCs that are hypermethylated and dark shading represents DMCs that are hypomethylated in stickleback reared at 21 ppt compared with 2 ppt respectively. have been previously identiﬁed as salinity-responsive using for regulating cellular ion concentrations in hyper and hypo- RNA-seq (Gibbons et al. 2017; Wang et al. 2014). osmotic conditions such as the calcium pump atp2b4,the 1,259 CpG loci were differentially methylated between sodium/chloride cotransporter slc12a3, and the sodium/po- salinity treatments (supplementary data set S4, tassium/2 chloride cotransporter slc12a1. Taken together, Supplementary Material online), and these DMCs were dis- these data suggest that changes in DNA methylation could tributed across all chromosomes, with an average of 0.01% play a role in facilitating the transition between marine and of the CpG loci on each chromosome being differentially freshwater environments. methylated (ﬁg. 4). No DMCs were identiﬁed between indi- In stickleback, a variety of genomic regions have been iden- viduals from low and high salinities in the mitochondrial ge- tiﬁed as having been subject to positive selection following nome. The distribution of DMCs across genomic features colonization of freshwater habitats by ancestral marine ﬁsh (e.g., promoters, exons, CpG islands) did not differ from the (Jones et al. 2012b). Because epigenetic variation has been relative proportions of these features in the genome (supple- suggested to be a driver of adaptive evolution (Flores et al. mentary table S2, Supplementary Material online). The ma- 2013), we screened our data set of salinity responsive DMCs jority of DMCs (1,051) was hypomethylated in individuals in marine ﬁsh to identify those associated with genes found in from high salinity relative to low salinity. Analysis of the genes regions under positive selection in freshwater populations located close to these DMCs revealed several genes known to (Jones et al. 2012b). Very few of the DMCs identiﬁed in our be involved in the response to salinity in ﬁsh (supplementary study were near these genes (supplementary data set S7, data set S5, Supplementary Material online). However, GO Supplementary Material online), suggesting that salinity- enrichment analysis did not detect signiﬁcant enrichment responsive changes in DNA methylation are unlikely to have for any GO categories following FDR correction. The ten GO played a role in driving genetic divergence in these regions terms with the lowest P-values are listed in supplementary between marine and freshwater populations of stickleback. tables S3–S5, Supplementary Material online. Comparison of genes located near DMCs in ﬁsh reared at Conclusions different salinities to previously identiﬁed as salinity-responsive In this study we used whole-genome bisulﬁte sequencing usingRNA-seq(Gibbons et al. 2017; Wang et al. 2014) iden- to identify novel DNA methylation features in the stick- tiﬁed 126 candidate genes with changes in both expression leback epigenome. Apparent hypermethylation of stra- and methylation in response to salinity (supplementary data tum two on the female X chromosome compared with set S6, Supplementary Material online). Among the candidate levels in males suggests that DNA methylation could play genes that we identiﬁed are ion channels that are important Genome Biol. Evol. 10(3):775–785 doi:10.1093/gbe/evy034 Advance Access publication February 6, 2018 783 Downloaded from https://academic.oup.com/gbe/article-abstract/10/3/775/4840697 by Ed 'DeepDyve' Gillespie user on 16 March 2018 Change in % methylation chr 1 chr 2 chr 3 chr 4 chr 5 chr 6 chr 7 chr 8 chr 9 chr 10 chr 11 chr 12 chr 13 chr 14 chr 15 chr 16 chr 17 chr 18 chr 19 chr 20 chr 21 % CpG per chromosome chr 1 chr 2 chr 3 chr 4 chr 5 chr 6 chr 7 chr 8 chr 9 chr 10 chr 11 chr 12 chr 13 chr 14 chr 15 chr 16 chr 17 chr 18 chr 19 chr 20 chr 21 Metzger and Schulte GBE Coulondre C, Miller JH, Farabaugh PJ, Gilbert W. 1978. Molecular basis of an important role in the suppressing recombination be- base substitution hotspots in Escherichia coli.Nature tween the X and Y chromosome, and potentially in reg- 274(5673):775–780. ulating sex-biased gene expression patterns. We also Devlin RH, Nagahama Y. 2002. Sex determination and sex differentiation detected signiﬁcant changes in DNA methylation in re- in ﬁsh: an overview of genetic, physiological, and environmental inﬂu- sponse to rearing salinity, some of which were associated ences. Aquaculture 208(3–4):191–364. Feng S, et al. 2010. Conservation and divergence of methylation pattern- with genes known to be differentially regulated in re- ing in plants and animals. Proc Natl Acad Sci U S A. sponse to changes in environmental salinity. This 107(19):8689–8694. epigenetic change reﬂects a response to environmental Flores KB, Wolschin F, Amdam GV. 2013. The role of methylation of DNA salinity that could facilitate the accumulation of epige- in environmental adaptation. Integr Comp Biol. 53(2):359–372. netic variation between natural populations, and thus be Franks SJ, Hoffmann AA. 2012. Genetics of climate change adaptation. Annu Rev Genet. 46:185–208. implicated in long-term responses to environmental Gertz J, et al. 2011. Analysis of DNA methylation in a three-generation change. family reveals widespread genetic inﬂuence on epigenetic regulation. PLoS Genet. 7(8):e1002228. Gibbons TC, Metzger DCH, Healy TM, Schulte PM. 2017. Gene expression Data Accessibility plasticity in response to salinity acclimation in threespine stickleback The sequencing FASTQ ﬁles from the whole genome bisulﬁte ecotypes from different salinity habitats. Mol Ecol. 26(10):2711–2725. Glazer AM, Killingbeck EE, Mitros T, Rokhsar DS, Miller CT. 2015. Genome sequencing can be downloaded from the NCBI sequence read assembly improvement and mapping convergently evolved skeletal archive (SRA accession: SRP127356). traits in sticklebacks with genotyping-by-sequencing. G3 5(7):1463–1472. Gorelick R. 2003. Evolution of dioecy and sex chromosomes via methyla- Supplementary Material tion driving Muller’s ratchet. Biol J Linn Soc. 80(2):353–368. Supplementary data areavailableat Genome Biology and Graves JAM. 2016. Evolution of vertebrate sex chromosomes and dosage compensation. Nat Rev Genet. 17(1):33–46. Evolution online. Heard E, Martienssen RA. 2014. Transgenerational epigenetic inheritance: myths and mechanisms. Cell 157(1):95–109. Hofmann GE. 2017. Ecological epigenetics in marine metazoans. Front Acknowledgments Mar Sci. 4:doi:10.3389/fmars.2017.00004. Ishikawa A, et al. 2017. Different contributions of local- and distant- We thank Timothy M. Healy for insightful discussions and regulatory changes to transcriptome divergence between stickleback comments during the preparation of this manuscript and ecotypes. Evolution (N Y) 71(3):565–581. for assistance with the analysis in R. Support for this research Jablonka E, Raz G. 2009. Transgenerational epigenetic inheritance: prev- alence, mechanisms, and implications for the study of heredity and was provided by a Natural Sciences and Engineering Research evolution. Q Rev Biol. 84(2):131–176. Council of Canada (NSERC) Discovery Grant to P.M.S., a UBC Jaenisch R, Bird A. 2003. Epigenetic regulation of gene expression: how Faculty of Graduate Studies Four Year Fellowship to D.C.H.M, the genome integrates intrinsic and environmental signals. Nat Genet. and a Zoology Graduate Fellowship to D.C.H.M. 33(3s):245–254. Jeong M, et al. 2014. Large conserved domains of low DNA methylation maintained by Dnmt3a. Nat Genet. 46(1):17–23. Literature Cited Jones FC, et al. 2012a. The genomic basis of adaptive evolution in threes- Akalin A, et al. 2012. methylKit: a comprehensive R package for the anal- pine sticklebacks. Nature 484(7392):55–61. Jones FC, et al. 2012b. A genome-wide SNP genotyping array reveals ysis of genome-wide DNA methylation proﬁles. Genome Biol. 13(10):R87. patterns of global and repeated species-pair divergence in stickle- Artemov AV, et al. 2017. Genome-wide DNA methylation proﬁling reveals backs. Curr Biol. 22(1):83–90. Kawahara R, Miya M, Mabuchi K, Near TJ, Nishida M. 2009. Stickleback epigenetic adaptation of stickleback to marine and freshwater con- ditions. Mol Biol Evol. 34(9):2203–2213. phylogenies resolved: evidence from mitochondrial genomes and 11 nuclear genes. Mol Phylogenet Evol. 50(2):401–404. Bell MA, Foster SA. 1994. The evolutionary biology of the threespine stick- leback. New York: Oxford University Press. Klose RJ, Bird AP. 2006. Genomic DNA methylation: the mark and its Bell MA, Stewart JD, Park PJ. 2009. The world’s oldest fossil threespine mediators. Trends Biochem Sci. 31(2):89–97. Kucharski R, Maleszka J, Foret S, Maleszka R. 2008. Nutritional control of stickleback ﬁsh. Copeia 2009(2):256–265. Bossdorf O, Richards CL, Pigliucci M. 2008. Epigenetics for ecologists. Ecol reproductive status in honeybees via DNA methylation. Science 319(5871):1827–1830. Lett. 11(2):106–115. Leder EH, et al. 2010. Female-biased expression on the X chromosome as a Boyko A, et al. 2010. Transgenerational adaptation of Arabidopsis to stress requires DNA methylation and the function of dicer-like proteins. PLoS key step in sex chromosome evolution in threespine sticklebacks. Mol Biol Evol. 27(7):1495–1503. One 5(3):e9514. Chatterjee A, et al. 2016. Sex differences in DNA methylation and expres- McCairns RJS, Bernatchez L. 2009. Adaptive divergence between fresh- sion in zebraﬁsh brain: a test of an extended ‘male sex drive’ hypoth- water and marine sticklebacks: insights into the role of phenotypic plasticity from an integrated analysis of candidate gene expression. esis. Gene. 590:307–316. Cooney CA, Dave AA, Wolff GL. 2002. Maternal methyl supplements in Evolution 64(4):1029–1047. McRae AF, et al. 2014. Contribution of genetic variation to transgenera- mice affect epigenetic variation and DNA methylation of offspring. J Nutr. 132(8 Suppl):2393S–2400S. tional inheritance of DNA methylation. Genome Biol. 15(5):R73. 784 Genome Biol. Evol. 10(3):775–785 doi:10.1093/gbe/evy034 Advance Access publication February 6, 2018 Downloaded from https://academic.oup.com/gbe/article-abstract/10/3/775/4840697 by Ed 'DeepDyve' Gillespie user on 16 March 2018 DNA Methylation Landscape of Stickleback GBE Melamed-Bessudo C, Levy AA. 2012. Deﬁciency in DNA methylation Takai D, Jones PA. 2002. Comprehensive analysis of CpG islands in human increases meiotic crossover rates in euchromatic but not in heterochro- chromosomes 21 and 22. Proc Natl Acad Sci U S A. 99(6):3740–3745. matic regions in Arabidopsis.Proc Natl Acad Sci U S A. Toli E-A, Calboli FCF, Shikano T, Merila J. 2016. A universal and reliable 109(16):E981–E988. assay for molecular sex identiﬁcation of three-spined sticklebacks Metzger DCH, Schulte PM. 2016. Maternal stress has divergent (Gasterosteus aculeatus). Mol Ecol Resour. 16(6):1389–1400. effects on gene expression patterns in the brains of male and Trucchi E, et al. 2016. BsRADseq: screening DNA methylation in natural female threespine stickleback. Proc R Soc B Biol Sci. populations of non-model species. Mol Ecol. 25(8):1697–1713. 283(1839):20161734. Varriale A. 2014. DNA methylation, epigenetics, and evolution in verte- Mirouze M, et al. 2012. Loss of DNA methylation affects the recombina- brates: facts and challenges. Int J Evol Biol. 2014:1–7. tion landscape in Arabidopsis. Proc Natl Acad Sci. 109(15):5880–5885. Vidalis A, et al. 2016. Methylome evolution in plants. Genome Biol. Morris MRJ, et al. 2014. Gene expression plasticity evolves in response to 17(1):264. colonization of freshwater lakes in threespine stickleback. Mol Ecol. Wan ZY, et al. 2016. Genome-wide methylation analysis identiﬁed sexually 23(13):3226–3240. dimorphic methylated regions in hybrid tilapia. Sci Rep. 6:35903. Nakamura R, et al. 2014. Large hypomethylated domains serve as strong Wang G, et al. 2014. Gene expression responses of threespine stickleback repressive machinery for key developmental genes in vertebrates. to salinity: implications for salt-sensitive hypertension. Front Genet. Development 141(13):2568–2580. 5:312. Peichel CL, et al. 2004. The master sex-determination locus in threespine Warton DI, Hui FKC. 2011. The arcsine is asinine: the analysis of propor- sticklebacks is on a nascent Y chromosome. Curr Biol. tions in ecology. Ecology 92(1):3–10. 14(16):1416–1424. White MA, Kitano J, Peichel CL. 2015. Purifying selection maintains Qu W, et al. 2012. Genome-wide genetic variations are highly correlated dosage-sensitive genes during degeneration of the threespine stickle- with proximal DNA methylation patterns. Genome Res. back Y chromosome. Mol Biol Evol. 32(8):1981–1995. 22(8):1419–1425. Wreczycka K, et al. 2017. Strategies for analyzing bisulﬁte sequencing Richards EJ. 2006. Inherited epigenetic variation—revisiting soft inheri- data. bioRxiv doi:10.1101/109512. tance. Nat Rev Genet. 7(5):395–401. Xie W, et al. 2013. Epigenomic analysis of multilineage differentiation of Ross JA, Peichel CL. 2008. Molecular cytogenetic evidence of rearrange- human embryonic stem cells. Cell 153(5):1134–1148. ments on the Y chromosome of the threespine stickleback ﬁsh. Yelina NE, et al. 2015. DNA methylation epigenetically silences crossover Genetics 179(4):2173–2182. hot spots and controls chromosomal domains of meiotic recombina- Ross JA, Urton JR, Boland J, Shapiro MD, Peichel CL. 2009. Turnover of sex tion in Arabidopsis. Genes Dev. 29(20):2183–2202. chromosomes in the stickleback ﬁshes (Gasterosteidae). PLoS Genet. Zemach A, McDaniel IE, Silva P, Zilberman D. 2010. Genome-wide evolu- 5(2):e1000391. tionary analysis of eukaryotic DNA methylation. Science Schultheiß R, Viitaniemi HM, Leder EH. 2015. Spatial dynamics of evolving 328(5980):916–919. dosage compensation in a young sex chromosome system. Genome Zhang W, Wang X, Yu Q, Ming R, Jiang J. 2008. DNA methylation and Biol Evol. 7(2):581–590. heterochromatinization in the male-speciﬁc region of the primitive Y Shao C, et al. 2014. Epigenetic modiﬁcation and inheritance in sexual chromosome of papaya. Genome Res. 18(12):1938–1943. reversal of ﬁsh. Genome Res. 24(4):604–615. Zhang X, et al. 2017. RNA-Seq analysis of salinity stress–responsive tran- Shen JC, Rideout WM, Jones PA. 1994. The rate of hydrolytic deamination scriptome in the liver of spotted sea bass (Lateolabrax maculatus). PLoS of 5-methylcytosine in double-stranded DNA. Nucleic Acids Res. One 12(3):e0173238. 22(6):972–976. Ziller MJ, Hansen KD, Meissner A, Aryee MJ. 2015. Coverage recommen- Smith G, Smith C, Kenny JG, Chaudhuri RR, Ritchie MG. 2015. Genome- dations for methylation analysis by whole-genome bisulﬁte sequenc- wide DNA methylation patterns in wild samples of two morphotypes ing. Nat Methods. 12(3):230–232. of threespine stickleback (Gasterosteus aculeatus). MolBiolEvol. 32(4):888–895. Associate editor: Kateryna Makova Genome Biol. Evol. 10(3):775–785 doi:10.1093/gbe/evy034 Advance Access publication February 6, 2018 785 Downloaded from https://academic.oup.com/gbe/article-abstract/10/3/775/4840697 by Ed 'DeepDyve' Gillespie user on 16 March 2018
Genome Biology and Evolution – Oxford University Press
Published: Mar 1, 2018
It’s your single place to instantly
discover and read the research
that matters to you.
Enjoy affordable access to
over 18 million articles from more than
15,000 peer-reviewed journals.
All for just $49/month
Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly
Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place.
All the latest content is available, no embargo periods.
“Whoa! It’s like Spotify but for academic articles.”@Phil_Robichaud