Environmental enrichment induces intergenerational behavioural and epigenetic effects on fish

Parental effects influence offspring phenotypes through pre‐ and post‐natal routes but little is known about their molecular basis, and therefore their adaptive significance. Epigenetic modifications, which control gene expression without changes in the DNA sequence and are influenced by the environment, may contribute to parental effects. We investigated the effects of environmental enrichment on the behaviour, metabolic rate and brain DNA methylation patterns of parents and offspring of the highly inbreed mangrove killifish (Kryptolebias marmoratus). Parental fish reared in enriched environments had lower cortisol levels, lower metabolic rates and were more active and neophobic than those reared in barren environments. They also differed in 1,854 methylated cytosines (DMCs). Offspring activity and neophobia were determined by the parental environment. Among the DMCs of the parents, 98 followed the same methylation patterns in the offspring, three of which were significantly influenced by parental environments irrespective of their own rearing environment. Our results suggest that parental environment influences the behaviour and, to some extent, the brain DNA methylation patterns of the offspring.

The parental rearing environment can induce phenotypic modifications during early development, which can be long-lasting and potentially intergenerational (Burton & Metcalfe, 2014). A well known example is the effect of structural environmental complexity on behaviour (Braithwaite & Salvanes, 2005;Roberts, Taylor, & Garcia de Leaniz, 2011), physiology (Näslund et al., 2013), cognitive capacity (Salvanes et al., 2013) and brain structure (Kihslinger, Lema, & Nevitt, 2006) in fish. Fish make use of physical structures at different points of their life cycle (e.g., for spawning, sheltering, foraging), suggesting that structural complexity is an important ecological factor of their natural environment (Näslund & Johnsson, 2016). Captive fish reared in enriched environments have shown increased survival in the wild compared to those reared in impoverished environments (D'Anna et al., 2012;Roberts, Taylor, Gough, Forman, & Garcia de Leaniz, 2014), as well as enhanced cognitive capacity and behavioural flexibility (Salvanes et al., 2013;Spence, Magurran, & Smith, 2011;Strand et al., 2010). Environmental enrichment affects different aspects of fish biology from behaviour to growth performance or resistance to disease (Näslund & Johnsson, 2016) some of which, like brain growth, can be plastic (Näslund, Aarestrup, Thomassen, & Johnsson, 2012). In captive salmonids, environmental enrichment alters not only behaviour but also gene expression (Evans, Hori, Rise, & Fleming, 2015) and liver methylation profiles (Gavery et al., 2019), while in zebrafish it changes avoidance behaviour and gene expression related to neurodevelopmental plasticity, some of which depends on age (Manuel et al., 2015). The use of enrichment structures such as shelters lowers basal levels of cortisol (Näslund et al., 2013), a corticosteroid commonly used as stress indicator (Barton, 2002), together with resting metabolic rate (Vaz-Serrano et al., 2011). While it seems clear that enrichment has various fitness benefits, whether these are inherited is unclear. What is needed are controlled experiments in common and contrasting environments, but these are rare, particularly in fish (Burton & Metcalfe, 2014).
Kyrptolebias marmoratus (Poey, 1880) is a predominantly self-fertilising fish living in mangrove forests in North and Central America (Tatarenkov et al., 2017), occupying a varied range of mangrove fossorial microhabitats under tidal control .
We reared two generations of genetically-identical K. marmoratus in matched or mismatched parent-offspring environments with different levels of structural complexity to examine the intergenerational influence of environmental enrichment on metabolic rate and behaviour, and the potential role of epigenetic mechanisms (brain DNA methylation) in mediating environmentally-induced parental effects.

| Experimental design
We used a highly inbred (at least 20 generations of inbreeding) strain of K. marmoratus (R [Ellison, Allainguillaume, et al., 2012]) kept under standard laboratory conditions (25-27°C, 16%-18%, 12 hr light: 12 hr dark photoperiod). For the parental generation, eggs from five fish of similar size and age were reared individually until hatching, when larvae were randomly transferred to tanks enriched with shelter and artificial plants (enriched environment, n = 14) or to barren tanks (poor environment, n = 13), where they were kept for 10 months (Berbel-Filho, Rodríguez-Barreto, Berry, Garcia De Leaniz, & Consuegra, 2019). All individuals analysed were hermaphrodites and no males were produced by any of the parents for the duration of the experiment.
To standardise potential age-related parental effects, eggs were only collected from parents of similar age (7-10 months). The offspring of five genetically inbred parents (three from enriched and two from poor environments) were set up following a factorial design with matched or mismatched parent-offspring environments ( Figure S1; Table S1). The offspring consisted of 15 mismatched individuals (seven poor to enriched, and eight enriched to poor), and 13 matched individuals (eight enriched to enriched and seven poor to poor).

| Measurement of metabolic rate and cortisol
We measured basal metabolic rate (oxygen consumption) of 55 adults (14 parental enriched; 13 parental poor, eight offspring from enriched parents reared in enriched environment; 10 offspring from enriched parents reared in poor environment; five offspring from poor parents reared in poor environment; five offspring from poor parents reared in enriched environment), 30 of which were also analysed for epigenetic variation (Table S1). Basal metabolic rates were measured at 8 months of age using four identical close respirometers consisting of 120 ml sealed dark chambers filled with oxygen saturated autoclaved water. Two blank trials were carried out to confirm no leakage of oxygen. Oxygen levels were calibrated in trials using saturated oxygenated water (100% dissolved oxygen) and anoxic water (2% dissolved oxygen). Fish were fasted for 48 hr prior to acclimation for 20 hr. Oxygen consumption was measured once for each fish for 40 min after acclimation, with oxygen levels always above 60%. Chambers were drained and cleaned between runs. Basal metabolic rate was calculated taking into account the rate of oxygen decrease in the chamber, mass of the individual, volume of water and time of measurement (mg O 2 g −1 min −1 ). Averaged background respiration levels across runs was 12.34% (SD = ±9.71).
We used ultrasensitive graphene immunosensors (Barton et al., 2018) for measuring waterborne cortisol noninvasively from parents reared in both enriched and poor environments. For this, 120 ml of water were taken from the respirometer after each individual measurement of metabolic rate and kept at -80°C until the analysis. A total of 10 ml were centrifuged at 109 g rpm for 5 min, and 10 µl of the supernatant were pipetted onto the modified sensor surface. Electrochemical measurements were conducted with a potentiostat/galvanostat (Autolab), controlled with NOVA software as in Barton et al. (2018).

| Behavioural analyses
Neophobia (i.e., fear of novelty, measured as number of contacts and inspections of a novel object) and exploratory behaviour were assessed using a plastic test arena (7 cm depth × 7 cm width × 30 cm length) filled with 0.7 L of water. The arena was divided into six equally spaced zones: a covered acclimatisation section (zone 0) with a sliding opaque door, and five open test zones (5 cm each) without cover (zones 1, 2, 3, 4 and 5) delineated by marks at the top margins of the arena walls. A coloured toy block (0.5 cm depth × 0.5 cm width × 3.5 cm height) was glued at the middle of zone 3 to serve as a novel object ( Figure S2). Nine-month old fish were placed individually into zone 0 for 15 min acclimatisation, after which the removable gate was slowly lifted, and the fish behaviour recorded for 20 min with an overhead camera fixed 0.5 m above the arena. After the test period, tanks were drained, rinsed with ethanol and distilled water. Videos were analysed by the same person using BORIS v. 7.
1. 4 (Friard & Gamba, 2016) to ensure consistency. The following four behaviours were quantified for both parents and offspring: (a) latency (s) to exit the acclimatisation zone; (b) number of inspections within 3 cm of the novel object (i.e., individual facing towards the novel object for more than 3 s); (c) number of contacts with the novel object; and (d) number of movements between zones (activity).

| Statistical analysis
All statistical analyses were run in r v. 3. 4. 3. Cortisol levels and basal metabolic rate were analysed for the parents using a linear model with environment (poor versus enriched), and bodyweight as predictors. We used a GLM with a quassipoisson link to account for overdispersion for the parental behavioural count data (No. contacts, No. inspections and activity) and a gaussian link for latency as a function of environment and bodyweight.
To test for parental effects on the offspring phenotype, we only analysed those phenotypes significantly different between parental environments, using the same model structure as described above but including also the parental values and environment as predictors.
We used the multimodel approach implemented in the r package glmulti v 1.0.7 (Calcagno & de Mazancourt, 2010) for model selection, which tests all possible models and all interactions, and considered models within 2 AIC units as being equivalent. To take into account potential parentage effects, we first selected the best-fit model (highest Akaike weight) using glmulti and then ran generalized mixed-models including parent of origin as a random factor using mlmrev v.1.0-7. Models were tested for overdispersion and individual observations (fish ID) were also taken into account when models displayed overdispersion. Outliers were identified using the function aout.pois in the package alphaOutlier.

| Genome-wide DNA methylation data
All individuals were analysed at the same age (10 months old A second library was created using 14 individuals from the offspring (five from enriched to enriched environments, three from enriched to poor, three from poor to poor, and three from poor to enriched). The library followed the same procedures and sequencing conditions as the first library.

| Sequence quality and alignment
We assessed the quality of the sequences using fastqc (Andrews, 2010), trimming of adaptors and low-quality reads was done using the RRBS default parameters (function: --rrbs) in TrimGalore! (Krueger, 2016).
Reads were aligned to the Kryptolebias marmoratus reference genome (NCBI ASM164957v1) (Rhee et al., 2017) prior to in silico bisulphite conversion using bismark v0.17.0 (Krueger & Andrews, 2011), which was also used to perform cytosine methylation calls. We only considered methylation within CpG context for the downstream analysis (Feng et al., 2010) and included CpGs with a minimum coverage of ≥10 reads in each sample across the 30 individuals sequenced. Individuals were grouped into generations (parents/offspring) and environments (own/parental) (Table S2). Mapped reads from parents and offspring were simultaneously processed and compared using the r package methylkit v. 1. 10 (Akalin et al., 2012). All analyses were conducted on a local server running NEBC Bio-Linux 8.

| Differentially methylated cytosines and methylation patterns
We first assessed differentially methylated cytosines (DMCs) between parental environments (enriched versus poor), using logistic regression on quantitated normalised data with q-value <0.01 after multiple testing correction and >20% minimal CpG methylation difference (|ΔM|), using methylKit. To test whether the number of DMCs between environments were different from random expectations, we generated 4,000 random combinations of 16 parental individuals and tested for the number of DMCs for each combination following the same parameters as the ones described for the original grouping. All code used for methylation analysis and the permutation pipeline is avaiable at: https://github.com/waldirmbf/DifferentialMethylation_MethylKit.
We then analysed whether the DNA methylation patterns (hy- gene orthologs in panther v. 11 (Mi et al., 2016). panther was used to assess the probability of overrepresentation of pathways, biological processes and molecular functions within the list of zebrafish orthologs using, using Fisher's exact test false discovery ratio (p < .05).
We used genemania (Assenov, Ramirez, Schelhorn, Lengauer, & Albrecht, 2008) to construct interaction networks of the annotated genes affected by DMCs between parental environments. We used NetworkAnalyser from cytoscape v. 3.7.1 (Assenov et al., 2008) to identify central genes within the molecular network and PANTHER GO terms were used to identify biological process and pathways for the most connected genes (>10 connections) within the network.

| Parental physiology and behaviour
Parents reared in the enriched environment were larger than those from the poor environment at the time of analysis (t = 2.84, df = 1, p = .008). For the parental individuals, the multi-model selection ap-  Figure 1d; Table S3). The number of inspections of the novel object was significantly explained by bodyweight (estimate: -9.50; z-value: −2.21, df = 1, p = .02; Figure 1e; Table S3). The number of contacts with the novel object was significantly affected by environment, with individuals from poor environments having higher number of contacts (estimate: 0.92; z-value: 2.98, df = 1, p = .001) (Figure 1f; Table S3e). Two individuals, one form the enriched and one from the poor environment, were identified as outliers (p < .01) for the number of contacts (7 and 14 respectively) but rerunning the analyses without these two individuals did not change the significance or the direction of the difference in the number of contacts between groups. The same individual from the poor environment was identified as an outlier for inspections (10) and after its removal from the analysis, neither bodyweight (estimate: -2.82; z-value: -0.51, df = 1; p = .89) nor the environment (estimate: 0.05; z-value: 0.12; df = 1; p = .60) significantly influenced the number of inspections.
Variation in parental latency to leave the acclimatisation zone was not affected by body weight or rearing environment (estimate: 0.12; z-value: 1.08, df = 1, p = .10).

| Parental effects on offspring behaviour
The only plausible model of offspring activity after correcting for overdispersion included parental activity, which had a positive effect on the activity of the offspring (estimate = 5.28, z-value: 2.61, df = 1, p = .009) (Figure 2a,b; Table S4a).

| Differential methylation between parental environments
After quality filtering, approximately 378 million reads were retained (range: 6-25 million), averaging 12.5 million reads per sample. Of those, approximately 61.1% were uniquely mapped reads to the reference genome. Overall bisulphite conversion efficiency was 99.5% (Table S1). In total, we identified 5.5 million cytosine sites of (SD ± 158.5), indicating that the DMCs identified between environments was significantly higher than expected by chance (p < .001).

| Methylation patterns for parents and offspring
Of the 1,854 DMCs identified between parental environments, 724  (Table 1). Four of the DMCs overlapped annotated genes in the mangrove killifish reference genome (insbr, transmembrane protein 230-like, col25a1, col18a1) (Table 1). In addition, the methylation pattern of three of the five DMCs which maintained the parental methylation was significantly influenced by parental environment but not by the offspring environment or their interaction (Table 2; Figure 3).
When analysed separately by environment-specific context, 30 DMCs in the offspring originated from enriched environment, and 19 in the offspring originated from poor environment maintained their methylation percentage relatively to its parents regardless of the offspring environment within less than 10% change. Within those, the terms found to be overrepresented were related to multicellular organism development (GO:0007275) and cell differentiation (GO:0030154) ( Table S6). The main pathways identified were represented by 3-7 genes each, but none of them were found to be overrepresented (Table S5b).

| Gene ontology and network analysis
Network analyses revealed a highly integrated molecular network, with 303 genes matching with ortholog genes in the zebrafish database and 235 genes showing at least one interaction, mainly based on co-expression (87.29%) and shared protein domains (9.32%) ( Figure S4). The average number of connections among the 235 genes was 4.18 (SD ± 3.40). A total of 22 input genes had >10 connections. Gene ontology analysis indicated that some of these genes are involved in mechanisms of plastic responses, such as regulation of gene expression (hopx, smad7), cell differentiation (fl4) and perception of extracellular signals and activation of intracellular cascades (map2k1, dld). Genes related to nervous system development (nav3), involved in its processing and perception of external information (penkb) and cellular reaction (axon guidance, efnb2a) were amongst the most connected genes in the network affected by the DMCs (Table S7).

| D ISCUSS I ON
The
Our results indicate that parental fish reared in enriched environments have lower basal metabolic rates and waterborne cortisol levels. While metabolic rate did not appear to be related to the rearing environment of the parents, the tight correlation between waterborne cortisol levels and basal metabolic rates in parental fish, suggests that parents reared in enriched environments were less stressed and spent less energy to maintain basal metabolic rate than individuals reared in barren environments.
In fish, structural environmental enrichment tends to decrease activity, mainly due to increased sheltering (von Krogh, Sorensen, Nilsson, & Overli, 2010;Moberg, Braithwaite, Jensen, & Salvanes, 2011;Roberts et al., 2011), and exploratory activity and boldness tend to be positively correlated (Champneys, Castaldo, Consuegra, & Garcia de Leaniz, 2018;Mazué, Dechaume-Moncharmont, & Godin, 2015). Here, parents reared in enriched environments were more active, but also more neophobic, than individuals reared in poor environments, suggesting no clear boldness-exploratory relationship in response to environmental enrichment in K. mamoratus. In this sense, plastic behavioural responses during the ontogeny of this species have been previously suggested, based on its variable responses to conspecific presence and simulated predation risk (Edenbrow & Croft, 2013).
Although behavioural effects of environmental enrichment on fish are well known (Jonsson & Jonsson, 2014;Näslund & Johnsson, 2016), the understanding of their potential inter-or transgenerational effects is limited, as most studies have focused on one generation (Näslund & Johnsson, 2016). Our results indicate that offspring activity and neophobia were influenced by the parental rearing environment. In general, offspring from parents reared in enriched environments had higher activity levels, regardless of their own environment, suggesting a sustained parental effect on activity levels. Previous studies in this species suggested that life-history traits (offspring size), but not behaviour (exploration, boldness and aggression), were affected by the parental environment (Edenbrow & Croft, 2013). However, in mammals there is ample evidence of parental effects caused by environmental enrichment, where the offspring from enriched environments tend to be more exploratory (Mychasiuk et al., 2012), have increased learning capacity and memory formation (Bygren, 2013;Dell & Rose, 1987), than those reared in non-enriched environments. While in fish increased cognitive capacity due to environmental enrichment was known (Roberts et al., 2011;Salvanes et al., 2013), this is the first evidence of behavioural intergenerational (parental) effects.

| Environmental enrichment effect on DNA methylation
Our results revealed a strong effect of environmental enrichment on brain DNA methylation patterns, with 1,854 differentially meth- poor environments. Several studies have reported effects of environmental enrichment on brain growth (Näslund et al., 2012), cell proliferation (von Krogh et al., 2010, cognitive capacity (Salvanes et al., 2013;Spence et al., 2011), and gene expression levels (Evans et al., 2015) in fish. The functional analysis of the genes affected by DMCs in the parents indicated that biological processes involved with cellular differentiation/communication, brain and organism structural development were significantly overrepresented. Among the 42 genes overrepresented for cell differentiation, eight were involved on signalling transduction. Among metabolic processes, 12 genes affected by DMCs were involved on the transcription by RNA polymerase II, suggesting that could be regulating differential gene expression levels in the brain of individuals reared in enriched and poor environments.
Network analysis identified dld (deltaD) among the most connected genes. dld is involved in the Notch signalling system which in the zebrafish brain is required to maintain the normal cyclic patterns of gene expression (Artavanis-Tsakonas, Rand, & Lake, 1999).
Disruption of dld leads to defective development during embryogenesis (Oates & Ho, 2002). Notch signalling is also involved in blood-vessel formation, together with flt4 (fms-related tyrosine kinase 4) (Siekmann & Lawson, 2007), also identified as part of the network of genes affected by the DMCs. The transcription factor Smad7 and the hopx gene (involved in regulation of transcription by RNA polymerase II) were also two of the most connected genes of the network, potentially related to changes in expression patterns in the brain, together with penkb (proenkephalin b), an opioid receptor binding expressed in the zebrafish brain, associated to chemical synaptic transmission. The differential methylation patterns observed between rearing environments could therefore be related to differential transcription levels (or transcriptional noise reduction, depending of the genetic context) (Evans et al., 2015) in pathways related to brain development that might explain the physiological and behavioural differences observed between individuals reared in enriched and poor environments. Alternatively, differences in brain structure related environmental enrichment similar to those observed in other fish (Kihslinger et al., 2006) could have resulted in different proportions of cell-types in the brain of fish reared in the different environments, which could explain the differences in methylation patterns. Further research analysing gene expression and targeting specific epigenetic variants, for example using CRISPR/ Cas9 to examine the effects of selective methylation of CpG sites in specific promoters (Vojta et al., 2016), would expand the information on how the specific epigenetic variants found here may be affecting downstream phenotypes.
Due to epigenetic reprogramming during embryogenesis, only a small subset of epigenetic variants on the parents are likely to be transmitted to the offspring (Burggren, 2016;Illum, Bak, Lund, & Nielsen, 2018). DNA methylation changes during embryogenesis in K. marmoratus has a longer and later DNA methylation reprogramming period when compared to other fish and mammals (Fellous et al., 2018), which might represent an epigenetic window of environmental sensitivity. In the offspring, most of our results indicated a stronger effect of their own rearing environment than that of their parents on DNA methylation patterns. Thus, although there were clear effects of environmental enrichment on the brain DNA methylation patterns of the parents, these changes may not have influenced the germline to the same extent, as for steelhead trout reared in captive or semi-natural environments (Gavery et al., 2019), suggesting limited potential for epigenetically-mediated parental effects being transmitted transgenerationally and/ or scape epigenetic reprogramming. Yet, three DMCs maintained the same methylation patterns in both parents and offspring while additional DMCs maintained the parental methylation patterns in the offspring in a more environment-specific manner. To our knowledge, this is the first evidence of parental effects on the offspring epigenetic patterns caused by environmental enrichment in fish, extending previous results in mice, in which parental enrichment has shown to affect offspring brain weight, global methylation levels (Mychasiuk et al., 2012) and learning capacity (Arai & Feig, 2011).
In summary, our results reveal behavioural and, limited but significant, epigenetic parental effects in the offspring caused by environmental enrichment which, if maintained, could have long-term evolutionary implications.

ACK N OWLED G EM ENTS
We thank Chis Cunningham and Ed Pope for advice with methylation and respirometer analyses. We are also thankful to Ben Jarrett for help with behavioural and metabolic measurements and CSAR staff for help with the rearing of killifish and to Larissa Rodrigues for her support with bioinformatic pipelines. This work was sup-

All experiments were approved by Swansea University Ethics
Committee (permit STU_BIOL_30484_110717192024_3).