This is the accepted version of the article: Moraleda-Prados, J.; Caballero-Huertas, M.; Valdivieso, A.; [et al.]. «Epigenetic differences in the innate response after immune stimulation during zebrafish sex differentiation». Developmental and Comparative Immunology, Vol. 114 (January 2021), art. 103848. DOI 10.1016/j.dci.2020.103848 This version is available at https://ddd.uab.cat/record/233725 under the terms of the license Manuscript File Click here to view linked References 1 Epigenetic differences in the innate response after immune stimulation during 2 zebrafish sex differentiation 3 4 Moraleda-Prados J1Ψ, Caballero-Huertas M1Ψ, Valdivieso A1, Joly S1, Ji J2,3, Roher N2, 5 Ribas L1* 6 1 7 Institut de Ciències del Mar, Spanish National Research Council (CSIC) 8 Passeig Marítim de la Barceloneta, 37–49. 08003 Barcelona,Spain 2 9 Institut de Biotecnologia i Biomedicina (IBB) and Dep. de Fisiologia, Facultat de Biologia, 10 Universitat Autònoma de Barcelona, 08003 Barcelona, Spain 3 11 Present address: Department of Marine Science and Engineering, Nanjing Normal University, 12 Nanjing, China 13 14 Ψ 15 authors participate equally to this manuscript 16 17 *correspondence: 18 Dr Laia Ribas. Tel. +34–932309547; e-mail lribas@icm.csic.es 19 20 21 Keywords: gonads, immune, LPS, methylation, reproduction, sex ratio, zebrafish 22 23 24 25 26 27 28 29 30 31 32 33 1 34 Abstract (150 words) 35 Infections are able to trigger epigenetic modifications; however, epigenetic-mediating 36 infections in the immune system in fish is currently unavailable. Within this purpose, 37 zebrafish were immune-stimulated with three lipopolysaccharides (LPS) during sex 38 differentiation. Methylation patterns of three immune genes were studied by a candidate 39 gene approach together with gene expression analysis, and in adulthood, sex ratios were 40 determined. It was shown that the entrance of LPS was through the gills and 41 accumulated in the pronephros. Significant hypomethylation levels of CASP9 and a 42 significant CpG site for IL1β after Pseudomonas aeruginosa LPS exposure were found. 43 No methylation difference was observed for TNFα. Gene expression and correlation 44 data differed among studied genes. Sex ratios showed a feminization in dose and LPS 45 strain-dependent manner. Here, it is provided epigenetic regulatory mechanisms derived 46 by innate response and the first evidence of possible epigenetic interactions between the 47 immune and reproductive systems. 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 2 65 Introduction 66 Epigenetic modifications occur in the host genome after infections altering the 67 transcriptome and the corresponding signaling pathways (Marr et al., 2014). The 68 understanding of epigenetic host-pathogen interactions is a flourishing research field to 69 untangle the sophisticated cellular strategies (Gomez-Diaz et al., 2012; Zhang and Cao, 70 2019). In humans, plenty of studies give evidence of rapid change in DNA methylation 71 of the innate immune cells to face infections (Sinclair et al., 2015; Wiencke et al., 2016; 72 Pacis et al., 2019) and the importance of the epigenetic mechanisms in the immunity 73 reprograming (Netea et al., 2016; Binnie et al., 2020). In contrast, less data is found in 74 fish. Findings in the guppy (Poecilia reticulata) showed changes in DNA methylation 75 dynamics after host-parasite interactions (Hu et al., 2018). In Atlantic salmon (Salmon 76 salar) it is shown that stress response together with immune challenges altered 77 transcriptome and methylome of the gills through immune responses (Uren Webster et 78 al., 2018). In zebrafish (Danio rerio), viral infections showed that histone modification 79 were able to increase methylation levels of the gene promoters associated with an innate 80 immune response in head kidney, liver, spleen and heart (Medina-Gali et al., 2018). 81 Surprisingly, to our knowledge, there is a lack of data regarding the epigenetic- 82 mediating infections in the immune system in fish, and more in particular, in the fish 83 gonads. 84 85 The interaction between reproduction and the immune system is present in fish although 86 it is far to be understood. To prevent germ cells against pathogens, gonads are 87 considered as immune-privileged organs (Maddocks and Setchell, 1990). The presence 88 of immune cells has been described along the gonadal development in gilthead 89 seabream (Sparus aurata L.) (Chaves-Pozo et al., 2008; Chaves-Pozo et al., 2009), in 90 turbot (Scophthalmus maxima) (Ribas et al., 2016) and in European sea bass 91 (Dicentrarchus labrax) (Ribas et al., 2019). Further, the pro-inflammatory cytokines, 92 tumor necrosis factor alpha (TNFα) and interleukin 1 beta (IL1β) are important to 93 regulate goldfish (Carassius auratus) testicular steroid biosynthesis (Lister and van der 94 Kraak, 2002). In zebrafish, it is known that the undifferentiated gonad, first developed 95 as an ovarian-like organ to latter, half of the population develop into testicular 96 phenotype throughout the activation of apoptotic pathways (Liew and Orban, 2014). 97 The nuclear factor kappa beta (NF-κB) has a crucial role by blocking apoptotic 98 pathways promoting cell survival and not allowing to activate the testicular pathway 3 99 and, consequently, the ovarian-like organ became ovaries (Orban et al., 2009). By 100 treating zebrafish with heat-killed Escherichia coli during the sex differentiation (~15– 101 40 days post fertilization, dpf), a skewed sex ratio towards females was observed, thus 102 demonstrating the existence that the activation of the immune system has consequences 103 in the final sexual phenotype (Pradhan et al., 2012). Strikingly, this is the only available 104 study showing that the activation of the immune system during gonadal development is 105 able to skew sex ratios in fish. 106 107 In contrast to mammals in which one master gene leads the sex of an individual (the sex 108 determination region of the Y chromosome, Sry), in fish, the genetic but also 109 environmental factors are responsible to determine the final phenotypic sex (Penman 110 and Piferrer, 2008). In the last decade, it has been explored the role of epigenetics in the 111 fish gonads, most of them focus on the temperature effects, demonstrating that the final 112 sexual phenotype depends on the epigenetic-environment interactions (Piferrer, 2019). 113 114 In fish, the innate immune response is the most important mechanism to cope with 115 immune challenges (Press and Evensen, 1999), although the adaptive immunity cannot 116 be underestimated by gaining shreds of evidence during the last years of research 117 (Secombes and Belmonte, 2016). The innate immune response is the first response 118 when an infection occurs and requires the pathogen recognition through germline 119 encoded by the pathogen recognition receptors (PRR) (Rauta et al., 2014). This 120 recognition is fast and characterized by activating inflammatory pathways (Bayne and 121 Gerwick, 2001) in which several cytokines, such as IL1β together with interleukin 6 122 (IL6) and TNFα, are ones of the main players (Smith et al., 2000; Brocker et al., 2010). 123 Caspases, including CASP9, are also important actors in the innate responses in fish by 124 activating apoptotic pathways and the signaling pathways controlled by the 125 inflammasome (McIlwain et al., 2013; Shalini et al., 2015). 126 127 Bacterial lipopolysaccharide (LPS) is the major outer surface membrane component 128 present in almost all the Gram-negative bacteria and LPS is considered as a good model 129 for studying acute inflammation in fish as most of the regulatory genes are conserved 130 between zebrafish and mammals (Forn-Cuni et al., 2017). Nevertheless, the receptors 131 triggering cellular activation are not the same. Once it is recognized by the cells, the 132 response to LPS in fish is immediate as gene expression is altered in already 30 minutes 4 133 (min) (Novoa et al., 2009) and can last for several hours (h) (MacKenzie et al., 2008) or 134 even for days (Abdel-Mageid et al., 2020). 135 136 The main objective of the present study was to evaluate epigenetic modifications 137 undertaken/initiated by the activation of the innate immune system during the gonadal 138 development. Consequently, influences of early immune activation during sex 139 differentiation on the final sexual phenotypes were evaluated. To address it, zebrafish, a 140 well-known experimental animal, was used to develop an in vivo model by using 141 different LPS strains from three of the most common bacteria (i.e., E. coli, Aeromonas 142 hydrophila, and Pseudomonas aeruginosa). First, we ensured LPS entrance in the 143 zebrafish larva by using LPS conjugated with fluorescence; secondly, we performed 144 dose-response experiments to determine the sensitivity to LPS treatments and the final 145 consequences on sex ratios. Following, the DNA methylation patterns of three major 146 players of the fish innate immune system (i.e., IL1β, TNFα and CASP9) were studied in 147 larvae after LPS immune stimulation. For this purpose, a candidate gene approach 148 method (Multiplex Bisulfite Sequencing, MBS) that allows determining epigenetic 149 changes at the level of single CpG sites (Anastasiadi et al., 2018b), were performed. In 150 parallel, gene expression of the three immune genes together with their correlations with 151 methylation levels were determined. 152 153 Materials and methods 154 1. Ethics statement 155 The experimental protocol was approved by the Spanish National Research Council 156 (CSIC) Ethics Committee within the project AGL2015-73864-JIN and licensed by the 157 Bioethical Committee of the Generalitat de Catalunya under reference code 9977. 158 European regulations of animal welfare (ETS N8 123, 01/01/91) were respected 159 regarding fish maintenance. Likewise, ICM facilities were validated for animal 160 experimentation by the Ministry of Agriculture and Fisheries (certificate number 161 08039–46–A) in accordance with the Spanish law (R.D. 223 of March 1988). 162 163 2. Animal rearing conditions 164 AB (ZFIN ID: ZDB-GENO-960809-7) laboratory strain zebrafish were housed in the 165 animal facilities of the experimental aquariums at the Institute of Marine Sciences 166 (ICM-CSIC, Barcelona, Spain). Zebrafish were reared in an ad hoc closed recirculating 5 167 system built for this purpose with a water pump of 3000 l h-1 and a UV light system to 168 eliminate any possible bacteria in the water. This system was placed in a chamber with a 169 photoperiod of 12 h of light and 12 h of darkness, an air temperature of 26 ± 1°C and a 170 humidity of 50 ± 3%. Physicochemical parameters were monitored daily, staying at 171 appropriate conditions (Ribas and Piferrer, 2014), which included temperature (28 ± 172 0.2°C), pH (7.2 ± 0.5), conductivity (750–900 μS), and dissolved oxygen (6.5–7.0 mg l- 1 173 ). Sulfite, sulfate, nitrate and ammonia quality parameters were checked weekly using 174 commercial kits and periodically by the water analysis service of the ICM-CSIC. 175 176 At the time of mating for natural egg fertilization, one female per male coming from 177 different families were placed in crossing tanks in order to maintain the interfamily 178 variation observed in this fish species (Ribas et al., 2017a). The total number of 179 fertilized eggs was counted to guarantee fertility according to the reference values for 180 this species and fish post-hatch survival agreed with OECD’s guidelines for the Fish 181 Sexual Development Test (OECD, 2011). About 50 eggs were arranged per Petri dish 182 with E3 embryonic medium (pH 7.2 ± 0.5) supplemented with 0.1% methylene blue 183 (Sigma-Aldrich, Madrid, Spain) at 26 ± 1°C until 6 dpf. Then, 25–35 larvae were 184 transferred into 2.8 liter plastic tanks (Aquaneering, mod. ZT280) to avoid 185 masculinization due to rearing density effects (Ribas et al., 2017b). 186 187 Fishes were fed three times a day according to their stages of development: 6–15 dpf 188 received specific larvae dried food (17.FO.CE.0611, Sparos I&D Nutrition in 189 Aquaculture) and from 15 dpf onwards, they were fed with pellets of increasing size 190 (AquaSchwarz, Göttingen, Germany). In all stages up to 10 dpf, diets were 191 complemented with Artemia nauplii (AF48, INVE Aquaculture, Dendermonde, 192 Belgium). 193 194 3. Purification of LPS 195 Commercial LPS from E. coli and P. aeruginosa (L2630 and L9143, respectively, 196 Sigma-Aldrich, Madrid, Spain) were dissolved in Milli-Q water at 1 µg ml-1 and, 197 subsequently, kept at −20ºC. In order to obtain A. hydrophila LPS, bacteria were grown 198 on LB agar plates containing 50 µg ml-1 ampicillin. A single isolated colony was 199 inoculated first in 50 ml LB overnight (O/N) culture to later grown O/N in 2.5 l of LB 200 for the large-scale production. The A. hydrophila culture was centrifuged at 6,000 rpm, 6 201 4°C, 30 min in a Beckman JLA 8.1000 rotor. Then resuspended sequencially in PBS, 202 ethanol, acetone and diethyl ether to finally prepare the dried biomass. Then a phenol- 203 water extraction method was used to separate the LPS from the biomass into two phases 204 after performing four centrifugations at 4ºC. The LPS was further dialyzed at 4ºC for 205 seven days and its concentration was checked using NanoDrop One spectrophotometer 206 (Thermo Fisher Scientific Inc., Wilmington, DE, USA). The solution was diluted to 207 obtain a final working 1 μg μl−1 concentration. 208 209 4. LPS entrance to larvae body 210 The LPS conjugated with fluorescein isothiocyanate (E. coli O111: B4 LPS-FITC, 211 F3665, Sigma-Aldrich, Madrid, Spain) was dissolved in Milli-Q water at a 212 concentration of 1 µg ml-1 and, subsequently, kept in a refrigerator at 4ºC for further 213 uses. All procedures were done in the dark to avoid FITC degradation. To visualize the 214 entrance of LPS-FITC, 20 dpf larvae were bath exposed to E. coli LPS-FITC in 1.5 ml 215 glass vials for 3 h (1 larva per tube). The four groups were: control group (CT) (without 216 LPS) (N = 4), 150 µg ml-1of E. coli LPS (N = 4), and 150 (N = 4) and 750 µg ml-1 (N = 217 4) of E. coli LPS-FITC (Fig. S1). After exposures, larvae were fixed with 500 μl of 4% 218 paraformaldehyde in 0.2 M phosphate buffer saline (PBS) in glass vials. Then, the 219 samples were washed with 0.2 M PBS three times for ten min each. Once the fixation 220 was completed, the 4′, 6-diamidino-2-phenylindole dihydrochloride (DAPI) 221 fluorescence staining (D9542, Sigma-Aldrich, Madrid, Spain) was added. For this 222 purpose, 1 μl of DAPI was diluted in 3 ml of 0.2 M PBS remaining in darkness. Then, 223 500 μl of this solution was placed in borosilicate vials for 5 min. Finally, larvae were 224 washed twice for one minute each with 0.2 M PBS and placed on a glass slide for 225 microscopical analysis under the Zeiss LSM880 microscope using the software ZEN 2 226 pro. 227 228 5. Immune stimulation experiments 229 5.1. Dose-response experiments 230 These dose-response experiments were performed to obtain a non-lethal dose capable of 231 trigger an immune response. LPS from three different bacteria species, E. coli, P. 232 aeruginosa and A. hydrophila, were used for dose-response experiments. A total of 16 233 larvae at 15 dpf per each of the two-three technical replicates were bath-immersed in 234 different LPS concentrations of the different strains (Fig. S1). Two, four, and three 7 235 biological family pairs were used for E. coli, P. aeruginosa and A. hydrophila, 236 respectively. The total number of larvae used for dose-response experiments were 192, 237 768, 1008 for E. coli, P. aeruginosa and A. hydrophila, respectively. For testing LPS 238 from E. coli, larvae were bathed in 0, 150 and 300 µg ml-1 LPS dilutions. For LPS from 239 P. aeruginosa, larvae were subjected to concentrations of 0, 150, 175 and 200 µg ml-1. 240 Lastly, the tested concentrations of A. hydrophila LPS were 0, 25, 50, 75, 150, 200 and 241 300 µg ml-1. 242 243 Larvae were placed in 3 ml glass containers for 24 h with the corresponding LPS 244 concentration for the three LPS strains, except for the control group. Glass vessels 245 containing larvae were incubated at 28 ± 1ºC in an oven chamber allowing circadian 246 rhythms (12 h light:12 h dark). After 3 h of P. aeruginosa LPS treatment, larvae were 247 washed with PBS and flash-frozen in liquid nitrogen and kept at -80ºC for further DNA 248 methylation (N = 10) and gene expression (N = 10) analyses. The same procedure was 249 followed for A. hydrophila LPS treated larvae after 3 h for gene expression analyses (N 250 = 10). During the following hours, survival for the threes LPS strains was recorded. 251 252 5.2. Immune stimulations during sex differentiation 253 Based on the observed survival results from dose-response experiments, to study sex 254 ratio resulting from immune stimulations during sex differentiation we repeated the E. 255 coli and P. aeruginosa treatments as previously reported with the 150 µg ml-1 256 concentration (Fig. S1). A set of larvae of 15 dpf were treated with E. coli LPS and 257 other two batches of larvae of 15 and 25 dpf were treated with P. aeruginosa LPS for 24 258 h. Between two to three technical replicates in one and two different biological family 259 pairs for E. coli and P. aeruginosa LPS, respectively, were used. The total initial larvae 260 used for this experiment was 96 and 320 for E. coli and P. aeruginosa LPS, 261 respectively. Once LPS treatments finished, larvae were reared in tanks and placed in 262 the rack system until adulthood to study sex ratios (120 dpf). 263 264 Due to the lack of sex ratio differences in adulthood observed from larvae exposed for 265 24 h of immune stimulation from both E. coli and P. aeruginosa strains, we extended 266 the treatment along gonadal development based on the methodology reported in Ribas 267 et al. (2017d) to obtain sex ratio differences. Thus, 18 dpf fish were immersed 7 times 268 during sex differentiation period (18–32 dpf) for 3 h using A. hydrophila LPS using 8 269 concentrations of 25, 75 and 150 µg ml-1 (Fig. S1). A total of seven family pairs (N = 2, 270 N = 2 and N = 3) for concentrations of 25, 75 and 150 µg ml-1 of A. hydrophila LPS, 271 respectively) were used. For each family pair, three technical replicates were performed 272 and a total of 816 initial number of fish was used. After 3h of immune stimulation 273 baths, larvae were washed using the same procedure as described above. At the end of 274 the treatments, fish of 32 dpf were reared in tanks and placed in the rack system until 275 120 dpf for sex ratio analyses. To evaluate possible stress originated from the 276 experimental procedures, a negative control group (mock) was added to the experiment 277 consisting of offspring of the same batch that was never removed from the tanks. 278 279 6. Methylation analysis 280 Methylation levels of larvae at 15 dpf treated 3 h with P. aeruginosa LPS at 150 281 μg ml−1 were studied for three immune genes i.e., IL1β, CASP9 and TNFα by MBS 282 technique following the procedures described elsewhere in (Anastasiadi et al., 2018b; 283 Caballero-Huertas et al., 2020). Briefly, genomic DNA from five larvae for each group 284 was extracted by 1 μg of proteinase K (Sigma-Aldrich, St. Louis, Missouri) overnight at 285 65 °C and the following day, a standard phenol-chloroform-isoamyl alcohol protocol 286 with ribonuclease A (PureLink RNase A, Life Technologies) was performed and 287 bisulphite converted. The targeted portion of the promoter was amplified with two PCR 288 rounds by using costumed primers for the three immune genes (Table S1). Primer 289 specificities were previously validated by Sanger sequencing of amplicons from a pool 290 of samples. Amplicon regions included the promoter, the first exon and the first intron 291 were included to the extent possible. Adaptor sequences for 16S metagenomic library 292 preparation (Illumina) were added to the 5ʹ ends of the primers designed: 293 Forward-TCGTCGGCAGCGTCAGATGTGTATAAGAGACAG and Reverse- 294 GTCTCGTGGGCTCGGAGATGTGTATAAGAGACAG. The target 295 regions included a total of nine for IL1β, 19 for CASP9 and seven for TNFα CpGs, 296 respectively. 297 298 Resulting PCR products were indexed by Nextera XT index Kit Set A (Illumina; FC- 299 131–2001) according to Illumina's protocol for 16 S metagenomic library preparation 300 and were pooled in an equimolar manner to obtain a single multiplexed library which 301 was sequenced in a MiSeq (Illumina, San Diego, California) using the paired-end (PE) 302 reads 250 bp protocol at the National Center of Genomic Analysis (CNAG, Barcelona). 9 303 Raw sequencing data were submitted in Gene Expression Omnibus (GEO) from NCBI 304 (https://www.ncbi.nlm.nih.gov) with the accession number: GSE134400. 305 306 7. Gene expression analysis 307 RNA was extracted from larvae treated for 3 h with P. aeruginosa LPS at 150, 175 and 308 200 µg ml-1 concentrations and with A. hydrophila LPS at 25, 75 and 150 µg ml-1. RNA 309 was individually extracted from 5 larvae each group with TRIzol (T9424, Sigma- 310 Aldrich, St. Louis, Missouri) according to manufactured procedures. RNA pellets were 311 suspended in 25 µl DEPC–water and kept at −80°C. RNA concentration was 312 determined by ND-1000 spectrophotometer (NanoDrop Technologies) and RNA quality 313 was checked on a 1% agarose/formaldehyde gel. By following supplier protocols, 200 314 ng of RNA were treated with DNAse I, Amplification Grade (Thermo Fisher Scientific 315 Inc., Wilmington, DE, USA) and retrotranscribed to cDNA with SuperScript III RNase 316 Transcriptase (Invitrogen, Spain) with Random hexamer (Invitrogen, Spain). 317 Quantitative PCR (qPCR) was carried out in technical triplicates for each sample with 318 the SYBR Green chemistry (Power SYBR Green PCR Master Mix; Applied 319 Biosystems). The conditions in the thermocycler: 50°C for 2 min, 95°C for 10 min, 320 followed by 40 cycles of 95°C for 15 sec and 60°C for 1 min in a 384-well plate (CFX- 321 386, Touch BioRad). The qPCRs were run in optically clear 384-well plates. Finally, a 322 temperature-determining dissociation step was performed at 95°C for 15 seconds (sec), 323 60°C for 15 sec and 95°C for 15 sec at the end of the amplification phase. Dissociation 324 step, primers efficiency curves and PCR product sequencing confirmed the specificity 325 for each primer pair. The qPCR primers used in this analysis are shown in Table S2. 326 327 8. Bioinformatics analysis 328 Raw data from the sequencer were demultiplexed based on the index codes by the 329 Illumina software and adapters were removed using the Trim Galore! software (v. 0.4.5) 330 (Babraham Bioinformatics). Quality controls of the samples were carried out during pre- 331 and post-trimming using the FastQC software (v. 0.11.8) (Andrews, 2010) to ensure 332 adapters were cut off correctly (Ewels et al., 2016). Low-quality bases were filtered 333 (Phred score <20) and only PE reads were used for this analysis. 334 10 335 We used in silico bisulfite-converted zebrafish genome (danRer11, GCA_000002035.4) 336 as a reference to align the PE using both steps procedures Bismark software (v.20.0) 337 (Krueger and Andrews, 2011). Bisulfite conversion efficiency was calculated for each 338 sample with a minimum threshold of 99.0%. All samples passed the minimum 339 threshold. The “BSgenome.Drerio.UCSC.danRer11” package was used to obtain the 340 coordinate positions of all CpG sites (TBD, 2019). Data were cleaned up, labelled and 341 tabulated using the Python 3 web-based environment Jupyter Notebook (v. 5.7.4). Every 342 single CpG site coordinate was checked after deleting those whose coverage was below 343 10 reads. Then, we calculated the methylation levels by averaging the CpG site values 344 in each gene for each sample and then were averaged by treatment. For individual CpGs 345 sites, the methylation values were averaged by the coordinate CpG position by the 346 treatment. 347 348 9. Statistical analysis 349 All statistical analyses were performed using RStudio (v. 1.1.456). Data were 350 expressed as mean ± S.E.M and the differences were considered significant when P 351 <0.05. Graphs were generated either by Sigma Plot software (v. S13.0) and by using the 352 “ggplot2” package (v. 3.1.0) (Wickham, 2009). 353 354 DNA methylation data 355 To work with a methylation database, the “dyplir” package were used (Wickham et al., 356 2020). DNA methylation differences between control and LPS groups were determined 357 by Student t-tests. Previously, homoscedasticity was checked by Levene’s test for every 358 single group, as well as normality was tested by the Shapiro-Wilk test for each group. 359 When normality was not followed, a Kruskal-Wallis test was performed. 360 361 Gene expression and correlation analyses 362 Data obtained from qPCR were collected by SDS 2.3 and RQ Manager 1.2 software. 363 For each sample, the relative quantity (RQ) values of the genes of interest were used to 364 normalize against the geometric mean value of two references genes (EF and RPL3A) 365 validated for zebrafish (Tang et al., 2007) and the fold change was calculated using the 366 2∆∆Ct method (Schmittgen and Livak, 2008). One-way ANOVA was used to detect 367 differences in gene expression between treatments. Previously, homoscedasticity was 368 checked by Levene’s test for every single group, as well as normality was tested by the 11 369 Shapiro-Wilk test for each group. When normality was not followed, a Kruskal-Wallis 370 test was performed. Tukey’s test was used to perform post hoc multiple comparisons. 371 Correlation analyses between methylation and gene expression were performed by a 372 Spearman’s rank correlation coefficient (𝑟𝑠 ) test using the corrplot package (v. 0.84). 373 374 Survival data from dose-response experiments 375 For dose-response experimental data, the Kolmogorov–Smirnov’s and Levene’s tests 376 were used to check data normality and the homoscedasticity of variances, respectively. 377 Then, a one-way analysis of variance (ANOVA) was used to detect possible differences 378 among groups. Tukey’s test was used to perform post hoc multiple comparisons. 379 380 Sex ratio statistics from immune stimulation experiments 381 Statistical significance of the resulting sex ratio was calculated by using the Chi-square 382 test (χ²) with the application of the Yates correction (Yates, 1934). 383 384 Results 385 1. LPS penetrates through the gills and recruited in the pronephros 386 Positive LPS-FITC signal was observed inside the larvae, in particular in the pronephros 387 area, after 3 h of bath incubation at 150 μg ml−1 (Fig. 1C). No signal was observed in 388 the negative control group neither in none-LPS conjugated (Fig. 1B and D). Likewise, 389 after bathing with a higher concentration (750 μg ml−1), fluorescence was observed in 390 the gill arches (Fig. 1E). To facilitate the observation of LPS localization inside the 391 larvae, figure 1A is a diagram of the dorsal view of illustrating the location of 392 pronephros, branchial arches and developing gonad. 393 394 2. Fish survival depends on LPS strain 395 E. coli LPS treatments did not affect larva survival as both concentrations used showed 396 similar results as the control group after 24, 48 and 72 h post-treatment (Fig. S2A). No 397 differences in larva survival were observed during the 24 h LPS treatment (data not 398 shown). 399 400 In P. aeruginosa LPS immune stimulation experiments, the larva survival rate behaved 401 in a dose-dependent manner (Fig. S2B). After 3 h onwards at 200 µg ml-1 a significant 402 decrease in survival was observed when compared to control group. At 175 µg ml-1 a 12 403 significant decrease in larva survival was observed after 24 h whereas no differences 404 were found at 150 µg ml-1. 405 406 Survival observed from A. hydrophila LPS experiments showed a clear dose-dependent 407 pattern although not enough biological replicates were obtained for all the 408 concentrations (Fig. S2C). Nevertheless, with the dose of 150 µg ml-1, a significant 409 decrease (P < 0.05) in fish survival was observed after 24 h. 410 411 3. LPS altered the DNA methylation of immune genes 412 After treating 15 dpf larvae for 3 h with P. aeruginosa LPS at 150 µg ml-1, no 413 significant differences in IL1β mean methylation were found between control and LPS 414 groups (Fig. 2A). However, at individual CpG sites, a significant (P = 0.00902) 415 methylation difference were found in the CpG9 site (position +189 after first ATG) with 416 hypomethylation in the LPS group (Fig. 2B). 417 418 Significant differences (P = 0.001305) in mean DNA methylation between control and 419 LPS groups were detected in CASP9 (Fig. 2C), being methylation levels in LPS group 420 lower than the control group. The same pattern was found in all CpG sites were 421 individually studied and concretely, CpG11 site (position -95 to first ATG) was 422 significantly hypomethylated in LPS group when compared to control (P = 0.01337) 423 (Fig. 2D). TNFα did not show significant differences in DNA methylation levels 424 between control and treatment groups both by the mean or by CpG sites analysis as 425 DNA methylation patterns were not altered (Fig. 2E and F). 426 427 4. Gene expression varies for each immune gene and depends on LPS strain 428 No significant up or downregulation in gene expression of the three studied genes were 429 found at 15 dpf larvae treated with P. aeruginosa LPS. IL1β gene was repressed more 430 than four and seven times after 3 h of incubation at 150 µg ml-1 and 175 µg ml-1 being 431 the former concentration almost significant (P = 0.0512) (Fig. 3A). 432 433 Gene expression of IL1β gene in larvae treated 3h with A. hydrophila LPS upregulated 434 significantly (P <0.01) for 75 and 150 µg ml-1 concentrations when compared to control 435 group. Gene expression was significantly upregulated in a dose-dependent manner with 436 more than twenty-fold change when compared to the lowest concentration (i.e., 25µg 13 437 ml-1) (Fig. 3B). No significance in gene expression differences was observed in CASP9 438 (Fig. 3D) while TNF expression was significantly downregulated in the highest dose 439 (Fig. 3F). 440 441 5. Negative correlation observed are mostly not significant 442 Correlation analyses were studied between DNA methylation and gene expression 443 levels in P. aeruginosa LPS treated larvae during 3h at 150 µg ml-1. Results for five out 444 of the six studied groups showed a negative tendency (Fig. 4). Nevertheless, only one 445 correlation was significant in TNFα LPS group (P = 0.01667; 𝑟𝑠 = -1). No significant 446 correlations between methylation and gene expression levels were found in IL1 (P = 447 0.35; 𝑟𝑠 = 0.6 and P = 0.6833; 𝑟𝑠 = -0.3, in control and LPS, respectively) or in CASP9 448 (P = 0.35; 𝑟𝑠 = -0.1 and P = 0.6833; 𝑟𝑠 = -0.3, in control and LPS, respectively). 449 450 6. Feminization depends on the duration of the exposure and the LPS strain 451 Percent of males of control groups of the eleven families used for immune stimulation 452 ranged among standard levels in AB strain (53–78%). Larvae of zebrafish treated during 453 24 h with E. coli or P. aeruginosa LPS at 150 µg ml-1 was not able to skew sex ratios in 454 populations (Fig. 5A and B). However, in the P. aeruginosa LPS treatments a 455 feminization tendency was observed. Likewise, 25 µg ml-1 repeated exposures for 3 h 456 during seven times along the gonadal development with A. hydrophila LPS did not 457 interfere sex ratio (Fig. 6A). Interestingly, those larvae treated with a concentration of 458 75 µg ml-1 of A. hydrophila LPS showed a significant (P <0.01) increase in the number 459 of females (Fig. 6B). In contrast, this feminization tendency was not observed after 460 increasing the concentration up to 150 µg ml-1 of A. hydrophila LPS but without 461 statistical significance (Fig. 6C). Nevertheless, a significant difference was found in the 462 mock group when compared to control and LPS treated (P <0.05). 463 464 Discussion 465 In the present work, we used the LPS, a molecule described as a good inflammatory 466 model in zebrafish (Yang et al., 2014), to activate the immune system in the developing 467 gonads to study DNA methylation impact and the possible consequences on the final 468 sexual phenotype. To ensure that the LPS penetrated the larvae body, as well as to 469 visualize which tissues it interacts with, larvae were submerged during 3 h in water with 14 470 E. coli LPS conjugated with FITC. It was observed that LPS penetrated through the 471 gills, with clear evidence in the highest LPS concentration. Penetration of labelled 472 antigens by bath immersion through the gills was observed in other fish species such as 473 rainbow trout (Oncorhynchus mykiss) or common carp (Cyprinus carpio L.) (Kiryu et 474 al., 2000; Huising et al., 2003). It is known that gills, together with the nose, the skin, 475 the gut and the urogenital system constitute the first line of fish defense (Rombout et al., 476 2014). Gills are mucosal tissues with rapid renewing allowing to recover from minor 477 epithelial damages (Lyndon and Houlihan, 1998), and after infections, locally secreted 478 factors allowing the migration of circulating leukocytes such as neutrophils (Huising et 479 al., 2003) or lymphocytes (Castro et al., 2014). In addition, LPS fluorescence was 480 detected near the head of the larvae where pronephros is located (Outtandy et al., 2019). 481 Pronephros consists of glomerulus connected by two bilateral ducts and persist in the 482 same position during the first month after fertilization (Drummond et al., 1998). The 483 role of pronephros is essential in the osmoregulation in fish but also in functions as an 484 immune organ with the presence of B lymphocytes and immunoglobulins two weeks 485 after fertilization (Zapata et al., 2006). In common carp, pronephros increased up to 486 53% of the platelet-forming cells when infected (Rijkers et al., 1980). After ontogenic 487 development, pronephros becomes the head kidney which is considered the principal 488 immune organ in fish and it is homologous to mammalian bone marrow (Meseguer et 489 al., 1995; Rauta et al., 2012). In adult fish, LPS insults in the head kidney increased cell 490 proliferation to cope with the infection (Ribas et al., 2008). Therefore, our observation 491 in the pronephros belongs to the recruitment of labeled LPS-FITC by the activation of 492 the immune cells responsible to cope with the insult. No fluorescence was detected in 493 the abdominal area where the differentiating gonads are located at any of the 494 concentrations of LPS tested, maybe due to the short exposure time or the experimental 495 procedure. Liposomes containing labeled LPS showed accumulation in zebrafish gills 496 after bath immersion but an intraperitoneal injection was required to be visible in the 497 internal tissues (i.e., spleen) (Ruyra et al., 2013). With the help of flow-cytometry, 498 observation of the liposome positive cells were also possible in the head kidney (Ruyra 499 et al., 2014). Thus, to identify labeled LPS in the developing gonads further 500 experimental procedures should be accomplished. 501 502 Once we knew that LPS was inside the larvae, we studied the DNA methylation levels 503 of three important innate immune-related genes at 15 dpf zebrafish when sex 15 504 differentiation takes place in this fish species. In all of the genes studied, LPS was able 505 to reduce DNA methylation levels on the promoter regions after only 3 h of treatment, 506 although only a significant alteration of mean DNA methylation was observed for 507 CASP9 gene. These results showed the rapid activation of the immune system by 508 epigenetic regulatory mechanisms. When looking in CASP9 at individual CpG sites, 509 only the CpG11 was differentially hypomethylated when compared to the control group 510 although hypomethylation tendency was observed in all the 19 CpGs studied. In 511 Atlantic salmon fry, in silico associations between transcriptomic data obtained after A. 512 aeruginosa LPS treatment with methylation data obtained after chronic stress treatments 513 showed that some immune genes, including CASP3a, suffered hypermethylation and 514 reduced expression, although other were hypomethylated (Uren Webster et al., 2018). 515 Thus, denoting that in fish there is a crosstalk between the environment and DNA 516 methylation of genes related to the immune system. 517 518 Here, despite not finding differences in IL1β global methylation levels between control 519 and LPS groups, there were significant DNA hypomethylation differences in the CpG9 520 site after LPS treatment. The rest of the CpG sites showed similar methylation levels 521 between control and LPS groups. In adult zebrafish, dimorphic differences in the 522 methylation patterns of IL1β between ovaries and testes were observed (Caballero- 523 Huertas et al., 2020). Strikingly, all CpGs studied in the mature zebrafish gonads 524 showed significant dimorphic differences except for the CpG9 site, in which similar 525 methylation levels were found in both sexes. Recent findings in mice gonads showed 526 dimorphic epigenetic marks in the Wnt/beta-catenin signalling pathway that allows 527 activating the ovary formation and maintenance (Koth et al., 2020). Overall, and 528 although more research is needed specifically in infected adult gonads, CpG9 site might 529 function as an important player of the immune response in fish and acts independently 530 to the sex of an individual. 531 532 LPS mimics bacterial infections and consequently, a cascade of inflammatory pathways 533 are triggered, including the expression of interleukins and cytokines (Bayne and 534 Gerwick, 2001; Novoa et al., 2009). Gene expression levels of IL1 gene was 535 significantly increased after A. hydrophila LPS immersion treatments but not 536 significance was found after P. aeruginosa LPS thus showing that A. hydrophila likely 537 stimulated the innate response to a greater extent than P. aeruginosa. Nevertheless, 16 538 IL1 gene expression at 175 µg ml-1 in P. aeruginosa LPS was almost significant. Thus, 539 in both LPS, IL1β expression was upregulated after only 3h of LPS treatments in a dose 540 dependent manner while no difference was found in CASP9. In contrast, TNFα was 541 inhibited at the highest dose after 3h of A. hydrophila LPS treatment. Another published 542 study in zebrafish larvae treated with E. coli LPS observed a decrease of IL1β 543 expression in a low dose but an increase in a high dose at 3 h (Novoa et al., 2009). 544 These dose-related differences could be linked with the tolerance response to LPS in 545 which it may protect the host from developing a shock syndrome caused by 546 hyperactivation of the immune system (Henricson et al., 1991; Medvedev et al., 2000). 547 Thus, when the immune system is hyperactivated, there is a suppression of the 548 production of many cytokines which can reprogramme the immune cells setting up an 549 adaptation to future bacterial infections (Shnyra et al., 1998). Thus, low doses of LPS 550 may exert beneficial effects by enhancing the protection against different fish diseases 551 and, in fact, LPS is used as adjuvant in some fish vaccines (LaFrentz et al., 2004; 552 Selvaraj et al., 2009; Vallejos-Vidal et al., 2016). 553 554 The MBS methodology here performed for DNA methylation analyses was based on a 555 high resolution at locus-specific but amplifies limited regions of DNA (~500 bp), and 556 although, targeted regions intended to include promotors, first intron and first exon, the 557 entire gene body and possible enhancer were not fully covered. Consequently, our 558 correlation analyses were not significant in any of the genes tested except for TNF in 559 the LPS group. Nevertheless, data showed slight inverse correlations between DNA 560 methylation levels and gene expression for CASP9 and TNFgenes following the 561 pattern dictated by the classical dogma; low DNA methylation of CpG-rich promoters is 562 associated with the activation of the gene transcription machinery (Jones and Takai, 563 2001; Deaton and Bird, 2011). Nevertheless, current research shows that DNA 564 methylation patterns are more dynamic and complex as hypomethylation of the 565 promoters is associated with gene silencing through the blockage of the transcription 566 initiation machinery (Ambrosi et al., 2017; Smith et al., 2020). Thus, other genomic 567 elements rather than the promoters contribute to transcriptional regulation like the exons 568 (Brenet et al., 2011), gene body (Blattler et al., 2014), the introns (Anastasiadi et al., 569 2018a), as well as post transcription modifications (Shilatifard, 2006). That could be the 570 case of the pattern observed in the IL1 correlation analysis. 17 571 572 In fish, abiotic factors such as temperature (Ribas et al., 2017a) or hypoxia (Shang et al., 573 2006) and biotic factors such as density (Ribas et al., 2017b, c) are able to skew sex 574 ratios towards males in the final population. In contrast, zebrafish populations treated 575 with chemical compounds i.e. demethylation agent 5-aza-2′-deoxycytidine (5-aza-dC) 576 (Ribas et al., 2017e), or with biotic factors such as E. coli heat-killed bacteria (Pradhan 577 et al., 2012) during sex differentiation, showed sex ratio bias towards females. In this 578 study, we have tested different bacterial strains at different concentrations to determine 579 whether activation of the immune system was able to alter the final phenotypic sex in 580 the developing zebrafish larvae. Larvae treated with E. coli and P. aeruginosa LPS at 581 150 µg ml-1 for a short period (i.e., 24 h) during sex differentiation did not have any 582 effect in the final sex ratios. However, repeated exposures along the sex differentiation 583 with A. hydrophila LPS at 75 µg ml-1 demonstrated a significant feminization of the 584 populations (P <0.01) in a similar manner previously described by Pradhan et al. 585 (2012). Although we did not find significant feminization at lower (25 µg ml-1) and 586 higher doses (150 µg ml-1), a tendency towards females was observed in the former. 587 Thus, higher number of biological family replications should be further tested. 588 Nevertheless, overall sex ratio results showed the existence of a crosstalk between the 589 immune system and the reproduction system in zebrafish. 590 591 Needs to be discussed that despite an observed trend towards the increase in the number 592 of females under 150 μg ml−1 A. hydrophila LPS, significant results were only found 593 when compared LPS treated and control groups to mock control group. Results 594 indicated that experimental procedure performed to treat larvae of zebrafish with LPS 595 along sex differentiation caused stress. It is known that cortisol released by stress 596 response is able to masculinize fish populations (Fernandino et al., 2012) and in the case 597 of zebrafish, synthetic cortisol was able to fully masculinized populations (Ribas et al., 598 2017b). Thus, based on the present data, LPS treatment might attenuate the 599 masculinization by counteracting the stress, thus reinforcing that the activation of the 600 immune system during sex differentiation is able to feminize zebrafish populations. Not 601 finding sex ratio differences in other LPS strains tested or other doses might be 602 explained by the fact that sensitivity to the LPS response depends on the LPS bacterial 603 origin species, serotype, doses, experimental concentrations and the time of duration of 604 the exposure (Novoa et al., 2009). It is known that fish, together with amphibians, are 18 605 resistant to the toxic effect of LPS as high doses are required to trigger an immune 606 response when compared to mammals (Berczi et al., 1966; Sepulcre et al., 2009) and 607 these differences were originated by diverse molecular strategies to recognize LPS (Chu 608 and Mazmanian, 2013; Yang et al., 2018). 609 610 Conclusions 611 Epigenetics integrates environmental cues and genomic information to determine the 612 final phenotype. Here, we have shown that the immune stimulation by LPS in the 613 developing zebrafish rapidly hypomethylates CASP9 and IL1 genes identifying, in the 614 former, a CpG which might play an important role in the immune system in zebrafish. 615 Gene expression data is altered only after 3h of immune challenges with higher 616 sensitivity of A. hydrophila when compared to P. aeruginosa. LPS entrance by 617 immersion is through the gills and recruited in the pronephros. Further, LPS is able to 618 skew sex ratios towards females but in a dose and strain-dependent manner. The present 619 work emphasizes the importance of epigenetic control in the response of zebrafish 620 immune system during gonadal development and reinforces the existence of the 621 crosstalk between reproductive and immune systems. 622 623 Acknowledgements 624 This study was mainly supported by the Spanish Ministry of Science and Innovation 625 (MICINN) projects AGL2015-73864-JIN “AMBISEX” together with RyC2018- 626 024017-I and 202030E004 “INTEROMICS” from Spanish National Research Council 627 (CSIC) grants to LR. Further, by Spanish Ministry of Science, European Commission 628 and AGAUR funds to NR (RTI2018-096957-B-C21, MINECO/FEDER and 2014SGR- 629 345 AGAUR). Special thanks to Dr. Francesc Piferrer for his mentoring in this study. 630 We thank Gemma Fuster, Dani Ortega, Esther Hoyo and Clara Barcelón for their help 631 with fish husbandry. 632 633 Conflicts of interest 634 The authors declare no conflict of interest. 635 636 Figure legends 637 Figure 1. Immune stimulation with E. coli LPS-FITC labelled. Images obtained from 638 fluorescence microscopy of 20 days post fertilization zebrafish larvae (blue by DAPI 19 639 labelling). (A) Dorsal view of a larva illustrating the location of pronephros, branchial 640 arches and developing gonad. (B) Control (without the addition of LPS). (C) Immune 641 stimulated larva with 150 µg ml-1 of LPS-FITC (green). (D) Immune stimulated larva 642 with 150 µg ml-1 of LPS. (E) Immune stimulated larva with 750 µg ml-1 of LPS-FITC 643 (green). 644 Figure 2. Percent of DNA methylation levels of immune genes of 15 dpf larvae 645 exposed to Pseudomonas aeruginosa LPS for 3 h. (A) IL1β percentage of mean DNA 646 methylation levels under P. aeruginosa immersion treatment. (B) DNA methylation 647 percentage for each of the nine CpG analyzed in IL1β. (C) CASP9 percentage of mean 648 DNA methylation under P. aeruginosa immersion treatment. (D) DNA methylation 649 percent for each of the 19 CpG analyzed in CASP9. (E) TNF percentage of mean DNA 650 methylation levels under P. aeruginosa immersion treatment. (F) DNA methylation 651 percentage for each of the seven CpG analyzed in TNFSignificant differences among 652 groups (P <0.05) were analyzed by Student t-tests and when normality was not 653 followed, a Kruskal-Wallis test was performed. 654 655 Figure 3. (A) IL1β, (C) CASP9 and (E) TNFα expression profiles after 3 h of treatment 656 with Pseudomonas aeruginosa LPS at different concentrations in 15 dpf larvae. (B) 657 IL1β, (D) CASP9 and (F) TNF expression profiles after 3 h of treatment with 658 Aeromonas hydrophila LPS at different concentrations in 15 dpf larvae. Data are shown 659 as mean ± SEM of fold change using control values set at 1. Sample size N = 5 larvae 660 per treatment. Significant differences (P <0.05) are symbolized by letters between 661 treated groups and by asterisks between treated and control groups: * (P <0.05), ** (P 662 <0.01). Statistics were analyzed by one-way ANOVA. 663 664 Figure 4. Correlations among DNA methylation and gene expression of (A) IL1β, (B) 665 CASP9 and (C) TNF genes. Correlations were analyzed by Spearman’s rank 666 correlation coefficient (𝑟𝑠 ). 667 668 Figure 5. Zebrafish sex ratio at 120 days post fertilization (dpf) after 24 h of 150 µg ml- 1 669 LPS challenges. (A) Sex ratio of fish treated with Escherichia coli LPS at 15 dpf. The 670 total number of fish are 35 and 64 in control and LPS groups, respectively. Each data 20 671 point is the mean ± SEM, corresponding to three technical replicates for one breeding 672 pair. (B) Sex ratio of fish treated Pseudomonas aeruginosa LPS at 15 and 25 dpf. Total 673 numbers of fish are 98 and 150 in 15 and 25 dpf, respectively. Each data point is the 674 mean ± SEM, corresponding to two to three technical replicates for each of the two 675 breeding pairs. Significant differences among groups (P <0.05) were analyzed by the 676 Chi-squared test with Yate’s correction. 677 678 Figure 6. Sex ratio of fish treated with Aeromonas hydrophila LPS during sex 679 differentiation (18-32 dpf). (A) Sex ratio results obtained at 25 µg ml-1 of A. hydrophila 680 LPS. Total numbers of fish are 121 and 116 in control and treated groups, respectively. 681 (B) Sex ratio results obtained at 75 µg ml-1 of A. hydrophila LPS. Total numbers of fish 682 are 84 and 87 in control and treated groups, respectively. (C) Sex ratio results obtained 683 at 150 µg ml-1 of A. hydrophila LPS. Total numbers of fish are 143, 146 and 72 in 684 control, treated and mock control group, respectively. Each data point is the mean ± 685 SEM, corresponding to breeding families (N = 2 families, N = 2 families, N = 3 families 686 for 25, 75 and 150 µg ml-1, respectively) for each of the two or three technical replicates 687 used for each concentration. Significant differences among groups (P <0.05) were 688 analyzed by the Chi-squared test with Yate’s correction. 689 690 691 692 693 694 695 696 697 698 699 700 701 702 703 704 21 705 References 706 Abdel-Mageid, A.D., Zaki, A.G., El Senosi, Y.A., Fahmy, H.A., El Asely, A.M., Abo- 707 Al-Ela, H.G., El-Kassas, S., 2020. Modulatory effect of lipopolysaccharide on immune- 708 related gene expression and serum protein fractionation in grey mullet, Mugil cephalus. 709 Aquacult. Res., 10. 710 711 Ambrosi, C., Manzo, M., Baubec, T., 2017. Dynamics and context-dependent roles of 712 DNA methylation. J. Mol. Biol. 429, 1459-1475. 713 714 Anastasiadi, D., Esteve-Codina, A., Piferrer, F., 2018a. Consistent inverse correlation 715 between DNA methylation of the first intron and gene expression across tissues and 716 species. Epigenet. Chromatin 11, 37. 717 718 Anastasiadi, D., Vandeputte, M., Sanchez-Baizan, N., Allal, F., Piferrer, F., 2018b. 719 Dynamic epimarks in sex-related genes predict gonad phenotype in the European sea 720 bass, a fish with mixed genetic and environmental sex determination. Epigenetics 13, 721 988-1011. 722 723 Andrews, S., 2010. FastQC: a quality control tool for high throughput sequence data. 724 Babraham Institute, Cambridge, United Kingdom. 725 726 Bayne, C.J., Gerwick, L., 2001. The acute phase response and innate immunity of fish. 727 Dev. Comp. Immunol. 25, 725-743. 728 729 Berczi, I., Bertok, L., Bereznai, T., 1966. Comparative studies on toxicity of 730 Escherichia coli lipopolysaccharide endotoxin in various animal species. Can. J. 731 Microbiol. 12, 1070-+. 732 733 Binnie, A., Tsang, J.L.Y., Hu, P., Carrasqueiro, G., Castelo-Branco, P., Dos Santos, 734 C.C., 2020. Epigenetics of sepsis. Crit. Care Med. 48, 745-756. 735 736 Blattler, A., Yao, L.J., Witt, H., Guo, Y., Nicolet, C.M., Berman, B.P., Farnham, P.J., 737 2014. Global loss of DNA methylation uncovers intronic enhancers in genes showing 738 expression changes. Genome Biol. 15. 22 739 740 Brenet, F., Moh, M., Funk, P., Feierstein, E., Viale, A.J., Socci, N.D., Scandura, J.M., 741 2011. DNA methylation of the first exon is tightly linked to transcriptional silencing. 742 Plos One 6. 743 744 Brocker, C., Thompson, D., Matsumoto, A., Nebert, D.W., Vasiliou, V., 2010. 745 Evolutionary divergence and functions of the human interleukin (IL) gene family. Hum. 746 Genomics 5, 30. 747 748 Caballero-Huertas, M., Moraleda-Prados, J., Joly, S., Ribas, L., 2020. Immune genes, 749 IL1 beta and Casp9, show sexual dimorphic methylation patterns in zebrafish gonads. 750 Fish Shellfish Immunol. 97, 648-655. 751 752 Castro, R., Bromage, E., Abos, B., Pignatelli, J., Granja, A.G., Luque, A., Tafalla, C., 753 2014. CCR7 is mainly expressed in teleost gills, where it defines an IgD(+)IgM(-) B 754 lymphocyte subset. J. Immunol. 192, 1257-1266. 755 756 Chaves-Pozo, E., Liarte, S., Fernandez-Alacid, L., Abellan, E., Meseguer, J., Mulero, 757 V., Garcia-Ayala, A., 2008. Pattern of expression of immune-relevant genes in the 758 gonad of a teleost, the gilthead seabream (Sparus aurata L.). Mol. Immunol. 45, 2998- 759 3011. 760 761 Chaves-Pozo, E., Liarte, S., Mulero, I., Abellan, E., Meseguer, J., Garcia-Ayala, A., 762 2009. Early presence of immune cells in the developing gonad of the gilthead seabream 763 (Sparus aurata Linnaeus, 1758). J. Reprod. Dev. 55, 440-445. 764 765 Chu, H., Mazmanian, S.K., 2013. Innate immune recognition of the microbiota 766 promotes host-microbial symbiosis. Nat. Immunol. 14, 668. 767 768 Deaton, A.M., Bird, A., 2011. CpG islands and the regulation of transcription. Genes 769 Dev. 25, 1010-1022. 770 771 Drummond, I.A., Majumdar, A., Hentschel, H., Elger, M., Solnica-Krezel, L., Schier, 772 A.F., Neuhauss, S.C.F., Stemple, D.L., Zwartkruis, F., Rangini, Z., Driever, W., 23 773 Fishman, M.C., 1998. Early development of the zebrafish pronephros and analysis of 774 mutations affecting pronephric function. Development 125, 4655-4667. 775 776 Ewels, P., Magnusson, M., Lundin, S., Kaller, M., 2016. MultiQC: summarize analysis 777 results for multiple tools and samples in a single report. Bioinformatics 32, 3047-3048. 778 779 Fernandino, J.I., Hattori, R.S., Kishii, A., Strussmann, C.A., Somoza, G.M., 2012. The 780 cortisol and androgen pathways cross talk in high temperature-induced masculinization: 781 the 11 beta-hydroxysteroid dehydrogenase as a key enzyme. Endocrinol. 153, 6003- 782 6011. 783 784 Forn-Cuni, G., Varela, M., Pereiro, P., Novoa, B., Figueras, A., 2017. Conserved gene 785 regulation during acute inflammation between zebrafish and mammals. Sci. Rep. 7. 786 787 Gomez-Diaz, E., Jorda, M., Peinado, M.A., Rivero, A., 2012. Epigenetics of host- 788 pathogen interactions: the road ahead and the road behind. PLoS Pathog. 8, e1003007. 789 790 Henricson, B.E., Neta, R., Vogel, S.N., 1991. An interleukin-1 receptor antagonist 791 blocks lipopolysaccharide-induced colony-stimulating factor production and early 792 endotoxin tolerance. Infect. Immun. 59, 1188-1191. 793 794 Hu, J.T., Perez-Jvostov, F., Blondel, L., Barrett, R.D.H., 2018. Genome-wide DNA 795 methylation signatures of infection status in Trinidadian guppies (Poecilia reticulata). 796 Mol. Ecol. 27, 3087-3102. 797 798 Huising, M.O., Guichelaar, T., Hoek, C., Verburg-van Kemenade, B.M.L., Flik, G., 799 Savelkoul, H.F.J., Rombout, J., 2003. Increased efficacy of immersion vaccination in 800 fish with hyperosmotic pretreatment. Vaccine 21, 4178-4193. 801 802 Jones, P.A., Takai, D., 2001. The role of DNA methylation in mammalian epigenetics. 803 Science 293, 1068-1070. 804 24 805 Kiryu, I., Ototake, M., Nakanishi, T., Wakabayashi, H., 2000. The uptake of fluorescent 806 microspheres into the skin, fins and gills of rainbow trout during immersion. Fish 807 Pathol. 35, 41-48. 808 809 Koth, M.L., Garcia-Moreno, S.A., Novak, A., Holthusen, K.A., Kothandapani, A., 810 Jiang, K., Taketo, M.M., Nicol, B., Yao, H.H.C., Futtner, C.R., Maatouk, D.M., 811 Jorgensen, J.S., 2020. Canonical Wnt/beta-catenin activity and differential epigenetic 812 marks direct sexually dimorphic regulation of Irx3 and Irx5 in developing mouse 813 gonads. Development 147. 814 815 Krueger, F., Andrews, S.R., 2011. Bismark: a flexible aligner and methylation caller for 816 Bisulfite-Seq applications. Bioinformatics 27, 1571-1572. 817 818 LaFrentz, B.R., LaPatra, S.E., Jones, G.R., Cain, K.D., 2004. Protective immunity in 819 rainbow trout Oncorhynchus mykiss following immunization with distinct molecular 820 mass fractions isolated from Flavobacterium psychrophilum. Dis. Aquat. Organ. 59, 17- 821 26. 822 823 Liew, W.C., Orban, L., 2014. Zebrafish sex: a complicated affair. Brief. Funct. 824 Genomics 13, 172-187. 825 826 Lister, A., van der Kraak, G., 2002. Modulation of goldfish testicular testosterone 827 production in vitro by tumor necrosis factor alpha, interleukin-1 beta, and macrophage 828 conditioned media. J. Exp. Zool. 292, 477-486. 829 830 Lyndon, A.R., Houlihan, D.F., 1998. Gill protein turnover: Costs of adaptation. 831 Comparative Biochemistry and Physiology a-Molecular and Integrative Physiology 119, 832 27-34. 833 834 MacKenzie, S., Balasch, J.C., Novoa, B., Ribas, L., Roher, N., Krasnov, A., Figueras, 835 A., 2008. Comparative analysis of the acute response of the trout, O-mykiss, head 836 kidney to in vivo challenge with virulent and attenuated infectious hematopoietic 837 necrosis virus and LPS-induced inflammation. BMC Genom. 9. 838 25 839 Maddocks, S., Setchell, B.P., 1990. Recent-evidence for immune privilege in the testis. 840 J. Reprod. Immunol. 18, 9-18. 841 842 Marr, A.K., MacIsaac, J.L., Jiang, R., Airo, A.M., Kobor, M.S., McMaster, W.R., 2014. 843 Leishmania donovani infection causes distinct epigenetic DNA methylation changes in 844 host macrophages. PLoS Pathog. 10. 845 846 McIlwain, D.R., Berger, T., Mak, T.W., 2013. Caspase functions in cell death and 847 disease. CSH Perspect Biol 5. 848 849 Medina-Gali, R., Bello-Perez, M., Martinez-Lopez, A., Falco, A., Ortega-Villaizan, 850 M.M., Encinar, J.A., Novo, B., Coll, J., Perez, L., 2018. Chromatin 851 immunoprecipitation and high throughput sequencing of SVCV-infected zebrafish 852 reveals novel epigenetic histone methylation patterns involved in antiviral immune 853 response. Fish Shellfish Immunol. 82, 514-521. 854 855 Medvedev, A.E., Kopydlowski, K.M., Vogel, S.N., 2000. Inhibition of 856 lipopolysaccharide-induced signal transduction in endotoxin-tolerized mouse 857 macrophages: Dysregulation of cytokine, chemokine, and Toll-like receptor 2 and 4 858 gene expression. J. Immunol. 164, 5564-5574. 859 860 Meseguer, J., López-Ruiz, A., García-Ayala, A., 1995. Reticulo-endothelial stroma of 861 the head-kidney from the seawater teleost gilthead seabream (Sparus aurata L): an 862 ultrastructural and cytochemical study. The Anatomical Record 241, 303-309. 863 864 Netea, M.G., Joosten, L.A.B., Latz, E., Mills, K.H.G., Natoli, G., Stunnenberg, H.G., 865 O'Neill, L.A.J., Xavier, R.J., 2016. Trained immunity: A program of innate immune 866 memory in health and disease. Science 352. 867 868 Novoa, B., Bowman, T.V., Zon, L., Figueras, A., 2009. LPS response and tolerance in 869 the zebrafish (Danio rerio). Fish Shellfish Immunol. 26, 326-331. 870 26 871 OECD, 2011. Test No. 234: Fish sexual development test, OECD guidelines for the 872 testing of chemicals, https://doi.org/10.1787/9789264122369-en. OECD Publishing, 873 Paris. 874 875 Orban, L., Sreenivasan, R., Olsson, P.E., 2009. Long and winding roads: Testis 876 differentiation in zebrafish. Mol. Cell. Endocrinol. 312, 35-41. 877 878 Outtandy, P., Russell, C., Kleta, R., Bockenhauer, D., 2019. Zebrafish as a model for 879 kidney function and disease. Pediatr. Nephrol. 34, 751-762. 880 881 Pacis, A., Mailhot-Léonard, F., Tailleux, L., Randolph, H.E., Yotova, V., Dumaine, A., 882 Grenier, J.-C., Barreiro, L.B., 2019. Gene activation precedes DNA demethylation in 883 response to infection in human dendritic cells. PNAS 116, 6938-6943. 884 885 Penman, D.J., Piferrer, F., 2008. Fish gonadogenesis. Part I: Genetic and environmental 886 mechanisms of sex determination. Rev. Fish. Sci. 16(S1), 16-34. 887 888 Piferrer, F., 2019. Epigenetics of sex determination and differentiation in fish. Sex 889 Control in Aquaculture, 65-83. 890 891 Pradhan, A., Khalaf, H., Ochsner, S.A., Sreenivasan, R., Koskinen, J., Karlsson, M., 892 Karlsson, J., McKenna, N.J., Orban, L., Olsson, P.E., 2012. Activation of NF-kappa ß 893 protein prevents the transition from juvenile ovary to testis and promotes ovarian 894 development in zebrafish. J. Biol. Chem. 287, 37926–37938. 895 896 Press, C.M., Evensen, O., 1999. The morphology of the immune system in teleost 897 fishes. Fish Shellfish Immunol. 9, 309-318. 898 899 Rauta, P.R., Nayak, B., Das, S., 2012. Immune system and immune responses in fish 900 and their role in comparative immunity study: A model for higher organisms. Immunol. 901 Lett. 148, 23-33. 902 27 903 Rauta, P.R., Samanta, M., Dash, H.R., Nayak, B., Das, S., 2014. Toll-like receptors 904 (TLRs) in aquatic animals: Signaling pathways, expressions and immune responses. 905 Immunol. Lett. 158, 14-24. 906 907 Ribas, L., Crespo, B., Sanchez-Baizan, N., Xavier, D., Kuh, H., Rodriguez, J.M., Diaz, 908 N., Boltana, S., MacKenzie, S., Moran, F., Zanuy, S., Gomez, A., Piferrer, F., 2019. 909 Characterization of the European Sea Bass (Dicentrarchus labrax) Gonadal 910 Transcriptome During Sexual Development. Mar. Biotechnol. 21, 359-373. 911 912 Ribas, L., Liew, W.C., Diaz, N., Sreenivasan, R., Orban, L., Piferrer, F., 2017a. Heat- 913 induced masculinization in domesticated zebrafish is family-specific and yields a set of 914 different gonadal transcriptomes. PNAS 114, E941-E950. 915 916 Ribas, L., Piferrer, F., 2014. The zebrafish (Danio rerio) as a model organism, with 917 emphasis on applications for finfish aquaculture research. Rev. Aquacult. 6, 209-240. 918 919 Ribas, L., Robledo, D., Gomez-Tato, A., Vinas, A., Martinez, P., Piferrer, F., 2016. 920 Comprehensive transcriptomic analysis of the process of gonadal sex differentiation in 921 the turbot (Scophthalmus maximus). Mol. Cell. Endocrinol. 422, 132-149. 922 923 Ribas, L., Roher, N., Martinez, M., Balasch, J.C., Donate, C., Goetz, F.W., Iliev, D., 924 Planas, J.V., Tort, L., MacKenzie, S., 2008. Characterization and expression of the 925 transcription factor PU.1 during LPS-induced inflammation in the rainbow trout 926 (Oncorhynchus mykiss). Fish Shellfish Immunol. 24, 35-45. 927 928 Ribas, L., Valdivieso, A., Diaz, N., Piferrer, F., 2017b. Appropriate rearing density in 929 domesticated zebrafish to avoid masculinization: links with the stress response. J. Exp. 930 Biol. 220, 1056-1064. 931 932 Ribas, L., Valdivieso, A., Diaz, N., Piferrer, F., 2017c. Response to "The importance of 933 controlling genetic variation remarks on 'Appropriate rearing density in domesticated 934 zebrafish to avoid masculinization: links with the stress response'". J. Exp. Biol. 220, 935 4079-4080. 936 28 937 Ribas, L., Vanezis, K., Antonio Imues, M., Piferrer, F., 2017d. Treatment with a DNA 938 methyltransferase inhibitor feminizes zebrafish and induces long-term expression 939 changes in the gonads. Epigenetics & Chromatin 10. 940 941 Ribas, L., Vanezis, K., Imues, M.A., Piferrer, F., 2017e. Treatment with a DNA 942 methyltransferase inhibitor feminizes zebrafish and induces long-term expression 943 changes in the gonads. Epigenet. Chromatin 10, 59. 944 945 Rombout, J.H.W.M., Yang, G., Kiron, V., 2014. Adaptive immune responses at 946 mucosal surfaces of teleost fish. Fish Shellfish Immunol. 40, 634-643. 947 948 Ruyra, A., Cano-Sarabia, M., Garcia-Valtanen, P., Yero, D., Gibert, I., Mackenzie, 949 S.A., Estepa, A., Maspoch, D., Roher, N., 2014. Targeting and stimulation of the 950 zebrafish (Danio rerio) innate immune system with LPS/dsRNA-loaded nanoliposomes. 951 Vaccine 32, 3955-3962. 952 953 Ruyra, A., Cano-Sarabia, M., MacKenzie, S.A., Maspoch, D., Roher, N., 2013. A novel 954 liposome-based nanocarrier loaded with an LPS-dsRNA cocktail for fish innate immune 955 system stimulation. Plos One 8. 956 957 Schmittgen, T.D., Livak, K.J., 2008. Analyzing real-time PCR data by the comparative 958 CT method. Nat. Protoc 3, 1101-1108. 959 960 Secombes, C.J., Belmonte, R., 2016. Overview of the Fish Adaptive Immune System, 961 in: Adams, A. (Ed.), Fish Vaccines, pp. 35-52. 962 963 Selvaraj, V., Sampath, K., Sekar, V., 2009. Administration of lipopolysaccharide 964 increases specific and non-specific immune parameters and survival in carp (Cyprinus 965 carpio) infected with Aeromonas hydrophila. Aquaculture 286, 176-183. 966 967 Sepulcre, M.P., Alcaraz-Perez, F., Lopez-Munoz, A., Roca, F.J., Meseguer, J., Cayuela, 968 M.L., Mulero, V., 2009. Evolution of lipopolysaccharide (LPS) recognition and 969 signaling: fish TLR4 does not recognize LPS and negatively regulates NF-kappa B 970 activation. J. Immunol. 182, 1836-1845. 29 971 972 Shalini, S., Dorstyn, L., Dawar, S., Kumar, S., 2015. Old, new and emerging functions 973 of caspases. Cell Death Differ. 22, 526-539. 974 975 Shang, E.H.H., Yu, R.M.K., Wu, R.S.S., 2006. Hypoxia affects sex differentiation and 976 development, leading to a male-dominated population in zebrafish (Danio rerio). 977 Environ. Sci. Technol. 40, 3118-3122. 978 979 Shilatifard, A., 2006. Chromatin modifications by methylation and ubiquitination: 980 Implications in the regulation of gene expression. Annu. Rev. Biochem. 75, 243-269. 981 982 Shnyra, A., Brewington, R., Alipio, A., Amura, C., Morrison, D.C., 1998. 983 Reprogramming of lipopolysaccharide-primed macrophages is controlled by a 984 counterbalanced production of IL-10 and IL-12. J. Immunol. 160, 3729-3736. 985 986 Sinclair, S.H.G., Yegnasubramanian, S., Dumler, J.S., 2015. Global DNA methylation 987 changes and differential gene expression in Anaplasma phagocytophilum-infected 988 human neutrophils. Clin. Epigenetics 7. 989 990 Smith, D.E., Renshaw, B.R., Ketchem, R.R., Kubin, M., Garka, K.E., Sims, J.E., 2000. 991 Four new members expand the interleukin-1 superfamily. J. Biol. Chem. 275, 1169- 992 1175. 993 994 Smith, J., Sen, S., Weeks, R.J., Eccles, M.R., Chatterjee, A., 2020. Promoter DNA 995 hypermethylation and paradoxical gene activation. Trends in Cancer. 996 997 Tang, R.Y., Dodd, A., Lai, D., McNabb, W.C., Love, D.R., 2007. Validation of 998 zebrafish (Danio rerio) reference genes for quantitative real-time RT-PCR 999 normalization. Acta Biochimica Et Biophysica Sinica 39, 384-390. 1000 1001 TBD, T., 2019. BSgenome.Drerio.UCSC.danRer11: Full genome sequences for Danio 1002 rerio (UCSC version danRer11). R package version 1.4.2 1003 30 1004 Uren Webster, T.M., Rodriguez-Barreto, D., Martin, S.A.M., Van Oosterhout, C., 1005 Orozco-terWengel, P., Cable, J., Hamilton, A., Garcia De Leaniz, C., Consuegra, S., 1006 2018. Contrasting effects of acute and chronic stress on the transcriptome, epigenome, 1007 and immune response of Atlantic salmon. Epigenetics 13, 1191-1207. 1008 1009 Vallejos-Vidal, E., Reyes-Lopez, F., Teles, M., MacKenzie, S., 2016. The response of 1010 fish to immunostimulant diets. Fish Shellfish Immunol. 56, 34-69. 1011 1012 Wickham, H., 2009. ggplot2: elegant graphics for data analysis. Springer-Verlag, New 1013 York. 1014 1015 Wickham, H., François, R., Henry, L., Müller, K., 2020. dplyr: A grammar of data 1016 manipulation. R package version 0.8.5. R package version 0.8.5. https://CRAN.R- 1017 project.org/package=dplyr. 1018 1019 Wiencke, J.K., Butler, R., Hsuang, G., Eliot, M., Kim, S., Sepulveda, M.A., Siegel, D., 1020 Houseman, E.A., Kelsey, K.T., 2016. The DNA methylation profile of activated human 1021 natural killer cells. Epigenetics 11, 363-380. 1022 1023 Yang, D., Zheng, X., Chen, S., Wang, Z., Xu, W., Tan, J., Hu, T., Hou, M., Wang, W., 1024 Gu, Z., 2018. Sensing of cytosolic LPS through caspy2 pyrin domain mediates 1025 noncanonical inflammasome activation in zebrafish. Nat. Commun 9, 1-12. 1026 1027 Yang, L.-L., Wang, G.-Q., Yang, L.-M., Huang, Z.-B., Zhang, W.-Q., Yu, L.-Z., 2014. 1028 Endotoxin molecule lipopolysaccharide-induced zebrafish inflammation model: a novel 1029 screening method for anti-inflammatory drugs. Molecules 19, 2390-2409. 1030 1031 Yates, F., 1934. Contingency tables involving small numbers and the χ 2 test. 1032 Supplement to the Journal of the Royal Statistical Society 1, 217–235. 1033 1034 Zapata, A., Diez, B., Cejalvo, T., Frias, C.G., Cortes, A., 2006. Ontogeny of the 1035 immune system of fish. Fish Shellfish Immunol. 20, 126-136. 1036 31 1037 Zhang, Q., Cao, X.T., 2019. Epigenetic regulation of the innate immune response to 1038 infection. Nature Reviews Immunology 19, 417-432. 1039 32 Figure 1 A Control B C 200 mm 200 mm 150 μg/ml LPS-FITC D E 200 mm 200 mm 150 μg/ml LPS 750 μg/ml LPS-FITC Figure 2 A B Control (CT) Lipopolysaccharide (LPS) a a IL1β a b C D a a b b CASP9 E F a a TNFα Figure 3 A B 0 30 ** b ** b (IL1β) 25 (IL1) *** ) change(IL1β) -2 ** change (IL1 20 -4 Fold change Foldchange a 15 -6 a ** 10 Fold Fold -8 a a* 5 -10 0 150 175 200 25 75 150 Concentrations (μg/ml) Concentrations (μg/ml) C D 3.0 3.0 a (CASP9) (CASP9) a change (CASP9) change (CASP9) a 2.5 ** a 2.5 2.0 ** 2.0 a a a Fold change Fold change 1.5 1.5 1.0 1.0 Fold 0.5 0.5 0.0 0.0 150 175 200 25 75 150 Concentrations (μg/ml) Concentrations (μg/ml) E F 3.0 0 (TNFα) (TNFα) ) aa ) 2.5 a aa change (TNF change(TNF a -2 2.0 aa Fold change Fold change 1.5 -4 1.0 a a -6 Fold 0.5 aa 0.0 -8 * 150 175 200 25 75 150 Concentrations (μg/ml) Concentrations (μg/ml) Figure 4 A B C Figure 5 A 100 a a 80 males Percentmales 60 Percent 40 20 0 Control LPS B a A A a Percent Percent males males Days post fertilization Figure 6 A 100 80 a males aa Percentmales a 60 Percent 40 20 0 Control 25 μg/ml LPS B 60 aa 50 Percent males Percent males 40 30 bb 20 10 0 Control 75 μg/ml LPS C 100 80 a aa males a Percent males bb 60 Percent 40 20 0 S e T iv C LP at Supplementary Material Click here to access/download Supplementary Material Moraleda et al 27082020 Supplementary.docx