Small ocean temperature increases elicit stage-dependent changes in DNA methylation and gene expression in a fish, the European sea bass

In natural fish populations, temperature increases can result in shifts in important phenotypic traits. DNA methylation is an epigenetic mechanism mediating phenotypic changes. However, whether temperature increases of the magnitude predicted by the latest global warming models can affect DNA methylation is unknown. Here, we exposed European sea bass to moderate temperature increases in different periods within the first two months of age. We show that increases of even 2 °C in larvae significantly changed global DNA methylation and the expression of ecologically-relevant genes related to DNA methylation, stress response, muscle and organ formation, while 4 °C had no effect on juveniles. Furthermore, DNA methylation changes were more marked in larvae previously acclimated to a different temperature. The expression of most genes was also affected by temperature in the larvae but not in juveniles. In conclusion, this work constitutes the first study of DNA methylation in fish showing that temperature increases of the magnitude predicted by the latest global warming models result in stage-dependent alterations in global DNA methylation and gene expression levels. This study, therefore, provides insights on the possible consequences of climate change in fish mediated by genome-wide epigenetic modifications.


Further details on temperature treatments of experiment 1
To achieve the desired temperature, and in order to avoid inducing a thermal shock, the nylon mesh bottom of the 19-liter PVC cylindrical container was blocked to prevent water loss and the container was gently moved to the next +1°C or -1°C mesocosm, depending on the polarity of change, until the desired temperature was achieved. These floating containers had an individual water supply to maintain water renewal at 30% vol . h -1 and two individual stone air diffusers to maintain dissolved oxygen saturation levels ~85-90%. The 650-liter tank had in turn its own air and water supplies (saturation levels ~85-100%). Thus, there were a total of twelve 650-liter mesocosms, the two at 14±0.5°C, which always received the fertilized eggs, and 10 additional mesocosms for the temperature experiments: two at 15±0.5°C, two at 16±0.5°C, two at 17±0.5°C, two at 18±0.5°C and two at 19±0.5°C (Fig. S2). Temperature changes were performed at steps of 1°C/2 h. In this way, groups passed through different 650-liter mesocosms. This procedure was repeated until the final designed temperature was reached (Fig. S2). Also, the tank effect was avoided since the 19-liter PVC cylindrical container passed through all the 650-liter mesocosms with the different temperatures (except for the 17°C group, which just passed through half of the 650-liter mesocosms).

DNA extractions
DNA was extracted from pools of larvae (experiment 1) and individually from lower body trunk of juveniles (experiment 2). The samples were immersed in digestion buffer (0.1 M NaCl, 10 mM Tris-HCl, 1 mM EDTA at pH=8, 0.5% SDS) with 1 μg of proteinase K (Sigma-Aldrich) and incubated overnight at 55ºC. DNA was extracted using phenol/chloroform/isoamyl alcohol (PCI) after treatment with 0.5 μg of ribonuclease A (PureLink RNase A; Life Technologies). DNA was precipitated by 95% ethanol, eluted in Milli-Q® water (Merck, Millipore) and measured by ND-1000 spectrophotometer (NanoDrop Technologies).

MSAP: laboratory methods
The MSAP method was a modification of the protocol used by Morán and Pérez-Figueroa 1 , in turn based on Reyna-López et al. 2 and Xu et al. 3 . Briefly, for each sample, digestion with restriction enzymes and ligation of adapters were performed simultaneously. For digestion, one aliquot of 100 ng of DNA was digested with 1 U of MspI and 5 U of EcoRI in Multicore buffer (Promega) containing 1 µg of BSA, whereas a second aliquot of 100 ng of DNA was digested with 1 U of HpaII and 5 U of EcoRI (Promega) in buffer A (Promega) also containing 1 µg of BSA. For ligation, 0.2 U of T4 ligase (Promega), 5 pmol of EcoRI adapter and 50 pmol of HpaII/MspI adapter (Table  S4) were included in each digestion reaction. The mixtures were incubated at 37ºC for 2 h. Adapters were prepared by carefully mixing equal volumes of oligos, denaturing at 95ºC for 5 min and annealing at room temperature for 10 min. After 1:5 dilution, the produced fragments were amplified using pre-selective primers (Table S4) containing one additional nucleotide by 1 U of GoTaq® Flexi DNA Polymerase (Promega). PCR cycling parameters were as follows: 72ºC for 2 min, followed by 20 cycles of 94ºC for 20 s, 56ºC for 30 s and 72ºC for 2 min, and finally 30 min at 60ºC. Two combinations of selective primers were used to produce a sufficient number of bands for robust statistical scoring. HpaII-MspI-TC primers end-labeled using 6-FAM reporter dye were combined with EcoRI-ACT primers (Table S4). On the other hand, HpaII-MspI-TC primers end-labeled using HEX reporter dye were combined with EcoRI-AAG primers (Table S4). Selective amplifications were performed by 1 U of GoTaq® Flexi DNA Polymerase using touchdown cycling parameters: 94ºC for 5 min, followed by 12 cycles of 94ºC for 20 s, 0.7ºC decreasing in each cycle starting at 66ºC for 30s and 72ºC for 2 min, in turn followed by 23 cycles at 94ºC for 20 s, 56ºC for 30s and 72ºC for 2 min and finally 30 min at 60ºC. For each sample and pair of restriction enzymes, 1 μl of selective PCR was loaded with 15 μl of Hi-Di™ Formamide (Applied Biosystems) and 0.5 μl of GeneScan™ 500 ROX™ dye Size Standard (Applied Biosystems) on an ABI 3130xl or 3730xl Genetic Analyzer (Applied Biosystems).

MSAP: raw data processing
Peaks along the electropherograms were detected and sized by the Peak Scanner™ Software v1.0 (Applied Biosystems). All subsequent analyses were performed using R (v. 3.2.5) and Rstudio (v. 1.0.136) 4,5 . The definition of bins, i.e., the amplicon size categories, final filtering and further quality controls were performed with the R package RawGeno 6 . Bins with over 100 standard raw fluorescence units, inside the 50-500 bp size range and reproducible over the 80% of replicated samples, were considered valid. Genotyping error per primer pair 7 was calculated by running two treatment groups, the 15ºC and 19ºC larvae (total of 35 samples), in duplicate and eight samples in triplicate from experiment 1 and seven samples in duplicate from experiment 2.

MSAP: classification of fragments according to methylation status
HpaII and MspI are isoschizomers that recognize the same nucleotide sequence (5'-CCGG-3') and have different sensitivity to the DNA methylation of the cytosines of the recognition sequence. Both enzymes cut unmethylated sequences and neither enzyme cuts hypermethylated sequences. In addition, MspI cuts the 5'-C me CGG-3'/3'-GG me CC-5' and HpaII cuts the 5'-me CCGG-3'/3'-GGCC-5' sequence [8][9][10][11][12] . Therefore, bands were scored as follows: if they were present in both the HpaII and MspI digestions they were considered unmethylated (Type I); if they were present only in MspI digestion they were considered methylated at the internal cytosine (Type II); if they were present only in the HpaII digestion they were considered hemi-methylated in the outer cytosine (Type III); and if they were absent from both the HpaII and MspI digestions they were considered hypermethylated (Type IV). Regarding type IV loci, genetic variation was expected to be low and therefore the absence of restriction sites was more probable due to hypermethylation. Loci were classified as methylation-susceptible (MSL) or nonmethylated (NML) if they were superior or inferior, respectively, to the Error Ratebased threshold (ERT). The ERT is specific for each primer combination, estimated as the proportion of discordant HpaII-MspI scores across individuals and indicative of scoring errors 13,14 . ERT was calculated by RawGeno as the mismatch error rate 7 . The error rate per primer was higher than in AFLP markers due to cell and tissue heterogeneity 14 . NML were too few to be considered as indicators of genetic variation and consequently they were eliminated from any further analysis. The binary matrices indicating band presence/absence for each enzyme per sample were transformed into binary matrices indicating methylated/unmethylated MSL per sample as follows: unmethylated (Type I) fragments were considered unmethylated, while hemimethylated, methylated in the internal cytosine and hypermethylated fragments were considered methylated (Types II, III and IV). Due to the potential error in identifying hemi-methylated outer cytosines as HpaII cuts 10 , one additional metric of global methylation was calculated as the ratio: methylated loci (Types II and III)/scorable loci (Types I, II and III), as previously reported 15,16 .

RNA extractions
Total RNA was extracted from 4 pools of ~10 whole larvae per group (experiment 1) and individually from upper body trunk of 4 juveniles per group (experiment 2) using the TRIzol reagent (Life Technologies) and precipitated by isopropyl alcohol. RNA quantity was measured by a ND-1000 spectrophotometer (NanoDrop Technologies) and diluted to 200 ng/μl. Equilibrated RNA concentrations were measured by a ND-1000 spectrophotometer and normalized a second time.

Quantitative real-time PCR (qRT-PCR)
One microgram of RNA was treated with 0.5 U of DNase I (Life Technologies) and reverse transcribed to cDNA by the SuperScript III Reverse Transcriptase (Life Technologies) and 100 μM of random hexamers according to the manufacturer's instructions. New primers were designed using Primer3Plus 17 . Primers used for nr3c1 were previously designed and validated in our group 18 , while primers targeting reference genes, the elongation factor-1 alpha (ef-1a) and the 40S ribosomal protein S30 (fau), had been previously validated in sea bass 19 . These two reference genes were used for normalization after evaluating also the ribosomal protein L13a (L13a) 19 and the ribosomal protein S18 (S18) and were chosen because of their stable expression among groups and an amplification efficiency close to 2. Primers sequences are shown in Table  S5.
All primers were validated separately in larvae and juveniles by performing serial dilutions (1, 1:5, 1:10, 1:50, 1:100, 1:500) of one pool of 5 cDNA samples from the 17ºC group and 5 cDNA samples from 15-19 (240) group and a second pool of 5 cDNA samples from the 17ºC group and 5 cDNA samples from the 21ºC group from experiment 2. These dilutions were used to perform calibration curves in order to calculate the slopes of log-linear regressions, correlation coefficients (r 2 ) and amplification efficiencies as E=10 (-1/slope) for each primer pair (Table S5). The specificity of the qPCR was confirmed by melting curve analysis (95ºC, 60ºC and 95ºC each for 15 s). SYBR Green (Life Technologies) chemistry was used in qRT-PCR. Reactions were performed in a total volume of 10 μl using 0.5 μl of each primer at 10 μM concentration, 5 μl of SYBR Green and 2 μl of cDNA, where the maximum dilution that produced Cq values <30 was used for each primer (Cq is used instead of Ct following MIQE guidelines). Reactions were run in triplicate and negative controls lacking cDNA were included in duplicate. qRT-PCR reactions were carried out on an ABI 7900HT (Applied Biosystems) unit with the following cycling parameters: UDG decontamination at 50ºC for 2 min, initial activation at 95ºC for 10 min and 40 cycles of denaturation at 95ºC for 15 s and annealing/extension at 60ºC for 1 min. A dissociation step (95ºC, 60ºC and 95ºC each for 15 s) was added at the end of each run.

Data availability
All data generated and analysed during this study are included as Supplementary Information files. These include the binary files used as inputs for the msap program, those with the methylation-sensitive polymorphic loci (outputs of the msap) on which the DNA methylation statistics were based and the dCq values of RT-qPCR. Figure S1. Schematic representation of the temperature treatments during development. The thermosensitive period in the European sea bass expands from 0 to ~60 days post fertilization (dpf). In experiment 1, fish were either reared at constant 15°C, 17°C or 19°C, or they were exposed to temperature changes at 15, 120 or 240 hours post fertilization (hpf) after an acclimation either at 15°C or 19°C. The treatments lasted from 0 to 15 dpf when all fish were sampled as larvae. The whole experiment was repeated 5 times. In experiment 2, natural fluctuations of temperature (around 17°C) occurred until 20 dpf. Then one group of larvae was subjected to constant 17°C or constant 21°C until 60 dpf when all fish were sampled as juveniles. Each group was carried out in duplicate. Figure S2. Experimental set-up and temperature treatments of experiment 1. Eggs from 5 different spawns coming from different broodstock kept in 5000-liter tanks were used. The eggs and larvae of each spawn were kept separate. Each batch of eggs was split in between three and five 19-liter cylindrical containers, depending on batch size, and first placed inside 650-liter fiberglass tanks at 14°C. There were 10 additional 650-liter tanks for the temperature experiments. Constant temperature experiments involved three groups reared at constant 15°C, 17°C, or 19°C. Six additional groups were exposed to temperature changes after acclimation at either 15°C or 19°C. Temperature changes took place at 15, 120 or 240 hours post fertilization (hpf) and all fish were sampled at 15 days post fertilization. Within a tank, temperature fluctuations were always <1°C. Pseudoreplication was avoided by using different families and replicate 19-liter cylinders for each family. Tank effects were avoided by having each temperature in duplicate and because the 19-liter tanks were moved across several 650-liter tanks to achieve the desired temperature. The subset of samples includes genomic DNA from pools of larvae reared at constant 17ºC, switched from 15ºC to 19ºC and from 19ºC to 15ºC at 120 hours post fertilization (hpf). For each sample, two separate double enzyme digestions were produced: with MspI/EcoRI and with HpaII/EcoRI. After a pre-selective PCR step, for each digestion two selective PCRs were performed with primers that included three selective nucleotides (ACT or AAG). At the same time, the second primer of the PCR was labeled with HEX or 6-FAM and included two selective nucleotides (TC). For samples 1, 2 … 19, HpaII/EcoRI digestions are shown after selective amplification with the AAG selective primer and HEX labeled (a), and with the ACT selective primer and 6-FAM labeled (b). In each case, in the first well the TrackIt 100 bp DNA Ladder (Life Technologies) was loaded, and in the last well, a negative control (NC) ran through the entire MSAP protocol without genomic DNA was loaded to the 2% agarose gels. Five microliters of the 20 μl total selective PCR reaction volumes were loaded for each sample. (c-d) Electropherograms after fluorescence-based capillary electrophoresis were visualized by Peak Scanner™ Software v1.0 (Applied Biosystems). For one sample (sample 2.2) of the 17ºC larvae group, one MspI/EcoRI (c) and one HpaII/EcoRI (d) digestions are shown after selective amplification with the ACT selective primer and 6-FAM-labeled (blue), and with the AAG selective primer and HEX-labeled (green). In all cases, the MSAP protocol was initiated with 100 ng of genomic DNA. Figure S4. Frequencies of methylation-susceptible polymorphic loci (MSL) band types. Shown are the percent of bands corresponding to unmethylated (Type I), methylated in the inner cytosine (Type II), hemi-methylated in the outer cytosine (Type III) and the ratio of Type II+Type III over the scorable loci (Type I+Type II+Type III) in the 15 days post-fertilization (dpf) larvae reared at 15ºC (15ºC-15 dpf) or at 19ºC (19ºC-15 dpf) or in the juveniles at 60 dpf reared at 17ºC (17ºC-60 dpf) or at 21ºC (21ºC-60 dpf). The sample sizes were as follows: Larvae: 15ºC, n=21, 19ºC, n=12; Juveniles: 17ºC, n=18, 21ºC, n=18. Fully methylated MSL are not included.