Grandparental immune priming in the pipefish Syngnathus typhle

Background Phenotypic changes in response to environmental influences can persist from one generation into the next. In many systems parental parasite experience influences offspring immune responses, known as transgenerational immune priming (TGIP). TGIP in vertebrates is mainly maternal and short-term, supporting the adaptive immune system of the offspring during its maturation. However, if fathers and offspring have a close physical connection, evolution of additional paternal immune priming can be adaptive. Biparental TGIP may result in maximized immunological protection. Here, we investigate multigenerational biparental TGIP in the sex-role reversed pipefish Syngnathus typhle by exposing grandparents to an immune challenge with heat-killed bacteria and assessing gene expression (44 target genes) of the F2-generation. Results Grandparental immune challenge induced gene expression of immune genes in one-week-old grandoffspring. Similarly, genes mediating epigenetic regulation including DNA-methylation and histone modifications were involved in grandparental immune priming. While grand-maternal impact was strong on genes of the complement component system, grand-paternal exposure changed expression patterns of genes mediating innate immune defense. Conclusion In a system with male pregnancy, grandparents influenced the immune system of their grandoffspring in a sex-specific manner, demonstrating multigenerational biparental TGIP. The involvement of epigenetic effects suggests that TGIP via the paternal line may not be limited to the pipefish system that displays male pregnancy. While the benefits and costs of grandparental TGIP depend on the temporal heterogeneity of environmental conditions, multigenerational TGIP may affect host-parasite coevolution by dampening the amplitude of Red Queen Dynamics. Electronic supplementary material The online version of this article (doi:10.1186/s12862-017-0885-3) contains supplementary material, which is available to authorized users.


Background
In sexually produced offspring, genotypes are determined by both maternal and paternal genetic contributions. An offspring phenotype is also influenced by a plethora of environmental factors experienced during its ontogeny and by its parents [1][2][3]. Such transgenerational plasticity of phenotypes is often adaptive, can promote efficient and rapid acclimatization to environmental changes, and even has the potential to modify evolutionary dynamics [4][5][6]. Anti-predator defenses [7], tolerance of abiotic environmental change [8][9][10], and induced disease resistance in offspring [11,12] are amongst the most studied transgenerational effects that are not inherited via DNA, but through a diversity of alternative mechanisms [13].
The transmission of parental parasite experience that subsequently leads to an enhanced offspring immune defense is known as transgenerational immune priming (TGIP) [14][15][16][17][18][19]. TGIP enables a faster or stronger offspring immune reaction that matches the current parasite environment [20]. Environmental variation can result in heterogeneous parasite distributions across environments [21] persisting through host generations. Under such matching environmental conditions where host-dispersal is limited and hosts have a long lifespan, selection for TGIP is predicted to be strong and evolutionarily adaptive [22,23].
In vertebrates, studies of TGIP have mainly focused on the transfer of maternal antibodies [14,17,24], while also substances of the innate immune system are involved [25][26][27]. Usually mothers deposit immune defense components into the eggs, transfer them during development (e.g. via the placenta), or, in mammals, after birth via lactation [14,17,28]. The classical view is that male sperm only contributes to heredity, i.e. via DNA, to the offspring. Considering recently discovered hereditary mechanisms that are not based on the pure DNA sequence itself, like DNA-methylation, histone acetylation pattern or tRNA, this view is challenged into question [29][30][31][32]. Recently, a growing number of examples indicate influences beyond pure transfer of DNA via sperm on offspring and, thus, underline the fathers' role beyond the determination of the offspring genotype [32][33][34][35][36][37][38]. In an invertebrate system with only an ejaculate-based connection between father and offspring, the paternal environment influences offspring immune phenotype [36,39]. Further, a strong paternal contribution to immune phenotype of progeny was found in vertebrates with intense paternal care or investment [19,35,40,41].
The advantages of TGIP in vertebrates were considered to be strongest during early life stages by strengthening the developing adaptive offspring immune system, with the effect fading upon maturation [42,43]. However, some ecological conditions may select for multigenerational TGIP modulating immune responses beyond the F1-generation [44]. Persistent TGIP should be favoured when the parasite environment is stable over time and, hence, across host generations. So far, our knowledge about TGIP in vertebrates past the early phase of an offspring's life is limited. Only scarce evidence supports the existence of TGIP beyond the maturation of the adaptive immune system of vertebrates [45,46], while in invertebrates TGIP can apparently cross the borders of more than just one generation [47].
To explore the potential for long lasting and multigenerational effects on immunity, we experimentally assessed grandparental TGIP in the sex-role reversed pipefish Syngnathus typhle. In this fish species, males have evolved a unique placenta-like structure [48]. Not mothers but fathers are thus the pregnant sex [49]. During male pregnancy, embryos are provided with nutrients and oxygen over this placenta-like structure [50][51][52][53], which may mechanistically enable a paternal investment into offspring immune defense. Usually in teleosts, females transfer immune components such as immunoglobulins, complement components, lectins, lysozymes and soluble antimicrobial peptides across follicle cells during the early stage of vitellogenesis into the oocyte [26,54]. However, apart from maternally derived immunity syngnathids profoundly rely on supplemental paternal immune priming [19,41]. As such, the pipefish system was chosen as here not only mothers but also fathers can induce offspring immune response over biparental TGIP [19,41]. So far, it was considered that TGIP in vertebrates is only of major importance during early development, to bridge the time of maturation of the acquired immune system, when selection pressure due to high mortality is greatest [55]. In contrast, recent studies indicate that in Syngnathus typhle the persistence of immune priming lasts past the maturation of the adaptive immune system in four-month-old juveniles [41], coupled with a high degree of bacteria specificity [40]. In the current study we aimed to address the impact of biparental parental immune priming beyond the generation border, affecting the immune dynamics of the grandoffspring generation.
Mature pipefish males and females (F0-generation) were exposed to two heat-killed bacteria (Vibrio spp. and Tenacibaculum maritimum) or a control prior to mating in a fully reciprocal mating design (Fig. 1). By leaving the F1-generation untreated, we were able to examine grandparental immune priming effects in the F2-offspring while challenging them with the same bacteria treatment as their grandparents (F0-generation). For the evaluation of grandparental sex-specific influences on grandoffspring immunity, either only grandmothers, only grandfathers or both grandparents (grandbiparental) were expsed to the bacteria treatment. In the F2-generation we assessed expression of 44 target genes functionally associated to different pathways of the immune system (innate and adaptive immune system, complement component system) and epigenetic regulation processes (DNA-methylation and histone modifications), to test whether complementing sex-specific contribution as previously found to exist for the F1-generation [40] may have been trans generationally maintained over two generations. Here, we found strong grandparental effects that influenced the immune gene expression of grandoffspring upon bacterial exposure. This grandparental TGIP is sex-specific (grandfather vs. grandmother) with regard to immune pathway activation and the involvement of epigenetic regulation genes.

Results
By using multivariate data analyses differential gene expression patterns were explored in one-week-old F2juveniles upon the applied F0-grandparental bacteria treatment (Vibrio: 'V+' and Tenacibaculum: 'T+' bacteria) in contrast to the naïve control group (Naïve: 'N'). We evaluated with a Permutational Multivariate Analysis of Variance (PERMANOVA) whether gene expression (44 target genes) of F2-juvenile pipefish revealed grandparental sex-specific influences ('F0-sex') and grandoffspring bacteria treatment effects ('F2-bacteria') including their interaction ('F0-sex x F2-bacteria'), while setting the family structure as random term. The multivariate PERMANOVA model was based on an Euclidean distance matrix and applied for 29 immune genes and 15 genes associated to epigenetic regulation, but also divided into following specific functional gene categories: (i) innate immune system (13 genes connected to the pro-inflammatory immune response), (ii) adaptive immune system (eight genes associated to the antibody-mediated immune defense), (iii) innate and adaptive immune genes (five genes connected to both immune pathways), (iv) complement system (three complement component genes that assist the antibody and phagocytic cell mediated immune response), (v) DNA methylation (five genes coding for DNA-methyltransferases), (vi) histone de/methylation (four histone de/methyltransferase genes), (vii) histone deacetylation (three histone deacetylation genes), and (vii) histone acetylation (two histone acetylation genes) [41]. Finally, we evaluated the contribution of variance explained by each target gene to identify central genes driving the grandparental bacteria treatment effect.

Complement component system (three genes)
The expression of complement component genes (C3, C9 and C1Q-sco) that mediate between innate and adaptive immune system and also complement the antibody-mediated immune response, showed grandparental sex-specific influences (PERMANOVA-complement: F 3,174 = 5.31, p = 0.017, Table 1, Fig. 2d). The complement component factors were impacted predominantly by the grand-maternal treatment (ANOSIM-complement: F0-Mat vs. F0-N p = 0.002, Table 2) but not by the grandpaternal treatment (ANOSIM-complement: F0-Pat vs. F0-N p = 0.124, Additional file 1: Table S2). As the grandmaternal and grand-biparental groups are not significantly different from each other and have the same center of gravity in the PCA (ANOSIM-complement: F0-Bi vs. F0-Mat = 0.168, Table 2 & Additional file 1: Table S2; Fig. 2d) the gene expression of complement factors of F2-juveniles was more affected by the grandmothers.   For further details see Additional file 1: Table S2 Multivariate ANOSIM was performed following significant PERMANOVA effects to assess differences in the gene expression profiles per treatment groups applying pairwise comparison on relative gene expression data (−ΔCt-values) based on a Euclidean distance matrix and 10000 permutations. Pairwise comparison was conducted for following two fixed factors and their interactions:

Gene contribution
Immune gene expression (29 genes-total) In the factor map the arrows of 29 immune genes were generally all pointing into the direction of the first principle component (Fig. 3a), which explains 30% of the total variation and visualizes in the corresponding PCA plots the grandparental treatment effect (Fig. 2a). Important genes with an average contribution above 60% were the innate immune genes Calreticulin, Transferrin, and Natural resistance-associated macrophage protein (Nramp), the adaptive immune genes Integrin, HIVEP2, and HIVEP3 as well as Complement component 3 and 9 (Fig. 3a). Most of the analyzed innate immune genes showed a contribution between 40-60% such as Peptidoglycan recognition protein, Heat shock protein 60 (Hsp60), Coagulation factor II, Lectin protein II, Kinesin, Allograft inflammation factor, Tyroproteinkinase, Ik-cytokine, Interferon, and Translocator protein (TSPO) (Fig. 3a). Besides, also following adaptive immune genes revealed a contribution between 40-60%: CD45, Tapasin (TAP), Lymphocyte cytosolic protein 2, Lymphocyte antigen 75, and Immunoglobulin light chain (Fig. 3a). Finally, Chemokine 7, LPS induced TNFα factor, Complement component 1, Interleukin 10, and Interleukin 8 contributed below 40% of the average variance (Fig. 3a).

Epigenetic regulation genes (15 genes-total)
Epigenetic regulation genes with the highest average contribution of over 80% were DNA(cytosine-5)-methyltransferases 3a and 3b (Fig. 3B), both responsible for de novo transfer of methyl groups to specific CpG sites in the DNA, permitting the formation of new methylation marks on unmethylated DNA [56][57][58]. Further, Histone acetyltransferase KAT2A (BROMO) which promotes acetylation of core histone proteins and with it transcriptional activation revealed over 80% of inertia contribution (Fig. 3b). These important key genes were followed with 60-80% of contribution by the maintenance DNA methyltransferase DNMT1 (Fig. 3b), which copies complementary marks of newly-replicated DNA by recognizing the hemimethylated sequences inherited from daughter strands [57].

Grandparental immune priming effect
The parental impact on offspring immune system persisted, once the adaptive immune system reached maturation [41] and even continued into the second generation, affecting immune dynamics of grandoffspring in the pipefish Syngnathus typhle. Grandparental immunological treatment with heat-killed bacteria epitopes altered the gene expression patterns in the grandoffspring, affecting the global expression of 29 immune genes. The patterns are crosscorrelated with the targeted functional gene groups, including genes of both innate and adaptive immune defense and complement component system. Upon immune challenge in the F2-generation, grandoffspring only induced immune gene expression if their grandparents already had experienced the bacterial epitopes. This result highlights the existence and importance of multigenerational TGIP, as only animals whose ancestors were exposed to a parasitic assemblage are able to quickly react towards an encounter with matching pathogens.
Genes of the innate immune system and complement component system were pre-dominantly affected, as they are essential in the pro-inflammatory response initiated 20 h after the bacterial injection. Central innate immune genes (>60% contribution) involved in driving this remarkable grandparental bacteria treatment effect were Calreticulin, Transferrin, Natural resistance-associated macrophage protein (Nramp), Complement component 3 and Complement component 9. The latter two genes are key players in the alternative pathway of the complement component system, necessary for the immediate immune defense against invaders [59]. Complement component 3 proteins recognize and tag bacteria and thereby activate the lytic pathway during which the membrane attack complex (MAC) is formed under the support of Complement component 9 [59]. Activation proteins such as Natural resistance-associated macrophage protein trigger macrophages activity that perform phagocytosis and remove pathogens with the production of reactive oxygen species (respiratory burst) and a release of hydrolytic lysozymes [60], while Calreticulin chaperones assist in promoting the phagocytosis process and clearance of apoptotic cells. A primitive but effective antimicrobial mechanism of the innate immune system during which actively nutriens (iron) are removed from bacterial pathogens, a process referred to as iron-withholding, is maintained by Transferrin and also intracellularly within the phagosome over Natural resistance-associated macrophage proteins [61,62]. Hence, the first line of innate immune defense is activated upon the grandparental treatment in the juvenile pipefish [63]. Further innate immune genes with a lower average contribution (between 40-60%) are also involved in immediately available and inducible pathways like bacteria recognition (C-type lectin II, Peptidoglycan recognition protein, Lectin protein II), antiviral response (Interferon), oxidative burst by macrophages (Translocator protein), stress response over molecular chaperone (Heat shock protein 60) as well as general inflammatory response (Allograft inflammation factor, Coagulation factor II, Kinesin, Tyroproteinkinase, Ik-cytokine) [64,65]. The assessed adaptive immune genes were not affected by the F0-bacteria treatment, potentially due to the ongoing maturation of the adaptive immune system in one-week-old pipefish [66] and/or the rather early time point of sampling (20 h after injection). Yet, our former data implys that the expression of the same immune genes used here positively correlates with an activation of the innate (monocytes) and adaptive (lymphocytes) cellular immune system [19,40,67], suggesting that induced gene expression is linked to a physiological impact of TGIP.
Grandparental sex-specific effect (grand-maternal vs. grand-paternal effects) Grandfathers and grandmothers might have evolved different strategies to achieve an optimal immunological protection of their grandoffspring. Here, we investigated complementing grandparental sex-specific contribution to different pathways of the immune system in oneweek-old F2-juveniles. Grand-paternal immunological experience boosted the activity of the innate immune system in their grandoffspring. This male-specific effect on innate immunity is consistent with the result of TGIP over one generation [40]. Grandfathers transfer immediate protection via innate immunity during male pregnancy against prevalent pathogens of their surrounding environment.
In contrast, grandmothers largely influenced genes of the complement component system. In teleosts mothers activate the complement system of their offspring via the deposition of a variety of diverse complement component proteins such as C1, C3, and C4 into the eggs [54,68]. Grandoffspring might even profit from mRNA transcripts that can likewise be transferred into the egg yolk [26]. A grand-maternal priming of the complement system can result in an earlier usage of acquired immune responses as it supports the antibody-mediated adaptive immune response. This implies that a stronger response against pathogens and parasites that linger in the environment over several generations is initiated. Since grandparents differentially influence the distinct immune pathways of the F2-offspring, grand-maternal and grand-paternal immune priming can complement each other leading to a balanced effect on total expression of 29 immune genes. Parental sex -specific influences on different immune pathway observed in former results of the F1-generation [40] lasted into the F2-generation and reflect an efficient strategy to maintain optimal protection against parasites by both parents correspondingly both grandparents.
Our data now suggest that information on prevalent bacterial epitopes is conserved and sex-specifically transferred, leading to complementing biparental TGIP over two generations. With this strategy parents not only deliver specific protection to their genetic offspring and grandoffspring, but they also transfer the opportunity to plastically adapt to the prevailing pathogen environment. In contrast to most species with conventional sex roles, in a sex-role reversed pipefish offspring are born into the paternal environment and, thus, share the paternal parasitic experience. This makes the transfer of immunity via the paternal line likely to be adaptive. Hence, both fathers and grandfathers will increase their fitness by altering their phenotype to optimally acclimatize offspring to the local parasitic environment [38]. Yet, female specific immune priming effects still remained, to a lower extent. Potentially this ancient evolution of beneficial maternal transfer of immune components into the eggs was still selected for, as a certain likelihood of matching parasite environment in their seasonal habitat (seagrass meadows) remained. When low dispersal might have resulted in habitat matching between grand-maternal and grandoffspring environment, maternal transfer of immunity should have been selected [22]. While most species boost their offspring immune response exclusively via maternal TGIP, pipefish rely on both, on maternal and paternal TGIP [41] that last for at least two generations. This dual developmental plasticity with sexspecific effector pathways gives the next generations an evolutionary advantage in reacting towards potentially virulent parasites. The grand-biparental TGIP is adaptive, provided that the maternal and paternal parasitic environment is experienced by their offspring correspondingly grandoffspring [22]. TGIP, accordingly, gives individuals an advantage whose ancestors successfully defended parasites and transfer this experience to the following generations.
Such heightened reaction to a pathogen should only be expressed under certain ecological circumstances as strong expression of immunity and the maintenance of inducible defense is costly due to high energy demand [69]. Bi-parental immune priming is traded off with delayed maturation time of adult F1-males, reduced fecundity and reproduction of the adult F1-generation, in case of parental bacteria exposure, indicating a compensational effect of reduced energy investment into reproduction [40,41]. These costs might constrain the overall beneficial net output of biparental immune priming [40,41]. Nevertheless, selection for grandparental immune priming designates that adaptive net influence and total benefits outweighed the associated costs.

Mechanism of immune priming (epigenetic regulation)
Our results cannot be explained by parasite-induced selection, as we used virulent heat-killed bacteria for the immune challenges. This presumes that the inheritance mechanism is non-DNA sequence based. The mechanisms permitting immunological information to be preserved via the paternal line over two generations most likely rely on a combination of small soluble immune components and epigenetic factors that are transferred via the sperm, the placenta-like structure or the fluid of the paternal brood-pouch tissue. As innate immune genes were predominantly influenced by the grandpaternal treatment, these genes might play a crucial role in the paternal transmission process.
DNA methylation and histone modifications are responsible for regulating packing and de-packing of the chromatin structure around histone molecules [70] and, consequently, the activation or deactivation of transcription processes for our targeted immune genes. That such epigenetic modifications of the genome can be responsible for paternal effects was recently demonstrated in zebrafish displaying paternal methylome transmission [34,37]. In our study, the expression of total 15 genes connected to DNA-methylation and histone de/methylation and de/acetylation in one-week-old F2 juveniles showed a significant change of expression in case of grandparental bacteria exposure. Epigenetic regulation genes that displayed a high contribution with over 80% were DNA-methyltransferase 3a, DNA-methyltransferase 3b. As de novo methylation via DNMT3a/b causes new chemical modifications of the DNA [56,58] and is essential for maternal and paternal imprinting [71], DNMT3a/b are potentially crucial mediators for epigenetic changes based on environmental stressors. In accordance to previous findings epigenetic regulation genes might not only be central regulators of parental immune priming [40,41], but also of grandparental immune priming, revealing a persisting effect into the second generation. In addition, we found strong evidence that histone acetylation genes (Histone acetyltransferase KAT2A (BROMO) and Histone acetyltransferase HAT1 (MYST)) regulating positively the accessibility of the DNA sequence for transcription processes by addition of acetyl groups to histone tails [72] were strongly influenced by the grandfathers solely. Histone modifications are supposed to be heritable across generations [73][74][75] and might carry epigenetic information [76]. The reaction to repeated pathogen exposure in macrophages involves positive histone marks and chromatin remodeling at specific promotors [77]. Moreover, it was suggested that histone modifications are associated with immune memory following a viral infection in CD8 T-cells [78]. The recent findings of parental [40] and grandpaternal influences on histone modification genes upon bacterial immune challenge suggests that the regulation of immune priming might be mediated with heritable marks stored on histones.
Our data, thus, propose that environmental stressors like pathogens leave an epigenetic mark on the genome affecting gene expression of genes associated with the immune system and transcriptional regulation that can be inherited over multiple generations. The fact that grandparental TGIP involves epigenetic mechanisms may result in a novel selection scenario for the evolution of TGIP along the paternal line, as the argument that male sperm is too small to transfer any more than just the DNA does not apply any longer. The sperm is thus potentially not only a sole messenger of "the other half of the offspring DNA", but also an important mediator for developmental plasticity and fast phenotypic acclimation to environmental changes [2,32,38,70].

Conclusions
Transgenerational effects on immunity in vertebrates are not short-term but can be sustained across two consecutive generations by the involvement of epigenetic regulation mechanism. These grandparental immune priming effects in the pipefish revealed complementing sex-specific contribution to different pathways of their grandoffspring immune system. Although TGIP might be beneficial on the individual level, it also involves ecological and evolutionary consequences on population level and has the potential to change disease dynamics and the spread of epidemics in a population [79,80]. Under negative-frequency dependent selection, rare parasite alleles may spread quickly in a population, while it takes time for the hosts as a population to counter-adapt under a Red-Queen dynamic [81]. In contrast, as an individual response, TGIP plays out within one generation, because the exposure to a novel parasite will already be met with an amplified immune reaction in the next generation. With this, the advantage of the novel parasite genotypes vanishes. This dampens the amplitude of predicted frequency dependent selection and may slow down Red Queen dynamics, giving the host an advantage in fast clearance of novel pathogen genotypes [79,80].
According to neutral genetic markers, Syngnathus typhle from different sites in the Baltic Sea all belong to the same population [82]. This is in line with the broad-nosed pipefish migratory behaviour: from open waters in winter to shallow seagrass meadows along the coastline in summer to exert mating and reproduction [83]. While in this sexrole reversed species, females display secondary sexual signals and are subject to multiple mating, males as the choosing sex are bound to their offspring during pregnancy. With respect to the seasonal migration pattern and the larval exposure to the pathogenic environment that their fathers already experienced, the investment into grandparental immune priming can be adaptive as individuals will be pre-adapted for the pathogen fauna in which the subsequent generations mate and release their offspring. On-going climate change with higher temperatures and lower salinity levels induces the abundance and virulence of pathogenic strains e.g. Vibrio particularly during the summer season [84]. The efficient transfer of immunological information about prevalent pathogenic threads is key for efficient short-term acclimation to changing virulence patterns [84] with benefits occurring particularly during the summer mating season in highly exposed shallow seagrass meadows. Future work should focus on the evaluation of resistance effects in a survival experiment, bacterial specificity processes, and an in depth analysis of the physiological mechanisms mediating grandparental TGIP.

Grandparental generation (F0-treatment)
Broad-nosed pipefish Syngnathus typhle were caught in the south-western Baltic Sea (54°44'N; 9°53'E, Germany) in spring 2013 and acclimatized to local summer conditions (15psu, 18°C, 14:10 h light:dark) within three weeks. We hosted the pipefish in local water out of Kiel Fjord, which was initially cleaned by a sand filter followed by 5, 20 and 50 μm mesh filter, UV-light, surface skimmer and biological filter to reduce the amount of microbes in the aquaria system. The parental generation (F0) was vaccinated as described previously [40], using a combination of heat-killed immunological novel Vibrio spp. and Tenacibaculum maritimum bacteria (F0-bacteria) or were left naïve (F0-N) as control. Immune-challenged mature pipefish were used in following mating design: 1.  Fig. 1). For the immune challenge, we used a combination of two distinct marine bacteria species to cover a potential wide range of immunological pathways, which could be differentially influenced by TGIP. The Vibrio spp. bacteria used in this experiment were an isolate of an Italian pipefish, allopatric and novel for the Baltic pipefish species [85]. The Tenacibaculum maritinum bacteria were an isolate of a pacific seabream species of Japan [86], and have, to our knowledge, not been in contact with the Baltic pipefish before. Both, Vibrio (s-shaped and flagellated) and Tenacibaculum (rod-shaped but non-flagellated), are common gram-negative marine bacteria causing the following diseases in teleost.
Tenacibaculum maritimum induces 'Flexibacteriosis' also known as 'black patch necrosis' in marine fish [87,88]. This disease is mainly characterised by haemorrhagic skin lesions, an ulcerative condition leading to important mortalities among marine fish species [87,88]. Vibrio bacteria can trigger 'Vibriosis' , a systemic disease of marine fishes [89], and e.g. Vibrio harveyi species are known to cause mass mortalities in captive bred seahorses [90]. The combination of Vibrio and Tenacibaculum permitted to cover an extended range of bacteria specific TGIP [41].

Filial generation 1 (F1-treatment)
F1-offspring were reared in 36×80 L aquaria and stayed separated in their tanks according to their parental treatment. Depending on their developmental stage, fish were fed with Artemia salina naupliae, copepods (Acarcia spec) and mysids (Mysis spec). F1-individuals were crossed within former parental treatment groups but left immunologically naïve (from each of the 4 parental treatment groups five families were chosen to do F1-crosses, resulting in 20 F1-families).

Filial generation 2 (F2-treatment)
In spring 2014, one-week old (post birth) F2-juveniles (F2) were exposed to the same heat-killed Vibrio (F2-V+) and Tenacibaculum (F2-T+) bacteria used for the F0generation or left naïve (F2-N) (per F1-crossing four families produced F2-offspring resulting in 16 F1-families). Out of each family, 12 individuals were chosen for the direct immune challenge. Per F2-offspring treatment (F2-V+, F2-T+, F2-N) four individual replicates were used, resulting in a total of 192 samples. Upon immune challenge, F2-juveniles were kept for 20 h in 10×10 cm tanks at 18°C and 15 psu, using one tank per F2offspring treatment and family. After the incubation time, juvenile body standard length [cm] was measured and animals were killed with MS 222. The body was transferred into 1 ml RNA-later, kept at 4°C for 24 h, and then frozen at −20°C.

Gene expression and data processing
We quantified the mRNA-level of 44 target genes and 4 housekeeping genes in 192 samples using quantitative real time PCR (qPCR) over a 96.96 dynamic array Fluidigm-BioMark™ system as described previously [40]. Thereby, the RNA extraction of 192 tissue samples and reverse transcription into cDNA was performed with a fixed amount of 800 [ng/μl] per sample as described previously [40].
For the following gene expression data analysisthe mean cycle time (Ct), standard deviation (SD), and the coefficient of variation (CV) were calculated. Samples with a CV larger than 4% were removed [91]. As the combination of the housekeeping genes ubiquitin (Ubi) and ribosome protein (Ribop) showed the highest stability (geNorm M > 0.85) [92], their geomean was used to quantify relative gene expression of each target gene by calculating − ΔCt-values [41]. Multivariate statistics were used to infer differences in the entire expression pattern of 29 immune genes and 15 epigenetic regulation genes, for more detailed evaluation the genes were also divided into following functional gene categories: (i) innate immune system, (ii) adaptive immune system. (iii) innate and adaptive immune genes, (iv) complement system, (v) DNA methylation, (vi) histone de/methylation, (vii) histone deacetylation, and (vii) histone acetylation [40,41].

Multivariate statistics
Statistical multivariate tests and plots were performed in R v 3.2.2 [93] and PRIMERv6 [94]. Grandparental sex specific influences ('F0-sex') (defined by four levels I: F0biparental bacteria treatment; II: F0-maternal bacteria treatment; III: F0-paternal bacteria treatment; IV: F0naïve no bacteria treatment) were evaluated by using F0sex as main factor and assessing its interaction with F2-bacteria treatment ('F0-sex x F2-bacteria'). Consequently, we fitted a PERMANOVA model ('vegan' package -'adonis' function in R) for each functional gene category (see last paragraph) based on an Euclidean distance matrix, by defining 'F0-sex' and 'F2-bacteria' treatments as fixed factors and stratifying permutations within each family replicate 10000 times (family was included as random factor). Standard length of F2juveniles was included as covariate in the model to correct for the dependence between gene expression and body size. Significant PERMANOVA tests were followed by an ANOSIM (Analysis of Similarity) with the software PRIMERv6 [94] which allowed pairwise comparisons between the different levels of F0-sex and F2-bacteria treatment groups as well as their interaction in a multivariate approach [95]. The ANOSIM was conducted likewise with an Euclidean distance matrix and 10000 permutations.
Principle component analysis (PCA) for graphical visualization was carried out based on an Euclidean distance matrix with the 'ade4' package in R [96]. For drawing the PCAs, we implemented the first three axes to obtain a projection of the whole data set onto a conveniently small dimension and to assess the clustering according to the F0-bacteria treatment due to differential gene expression. PCAs were solely performed for functional gene categories that revealed a significant F0-sex effect (Fig. 2). In addition, we evaluated the percentage of contribution of response variables (genes) in explaining the variations retained by the principle components (PCs) by applying the 'factoextra' package implemented in R [97]. The total contribution of a variable (gene) which explains the variations elicited by the principle components (PCs) was calculated within the function 'fvizcontrib' [97]. Gene contribution (%) was visualized using a factor map in which a implemented colour gradient highlights most important genes with the highest contribution of variance (Fig. 3) [97].

Additional file
Additional file 1: Table S1. Additional values of 2-way PERMANOVA output. Multivariate PERMANOVA analysis to assess the effect and interaction of two fixed factors F0-sex and F2-bacteria while including size as covariate and family as strata term on relative gene expression data (−ΔCt-values). Each analysis was based on an Euclidean distance matrix with p-values obtained by 10000 permutations. Significant p-values are marked in bold letters and asterix symbol (significance code: <0.001***, 0.001**, 0.01*, 0.1 > p-value ≥ 0.05 trend •). R 2 value indicate the percentage of variance explained by the model. Table S2. Results from PERMANOVA and ANOSIM analysis of one-week-old F2-juveniles per functional gene categories. Multivariate ANOSIM was performed following significant PERMANOVA effects to assess differences in the gene expression profiles per treatment groups applying pairwise comparison on relative gene expression data (−ΔCt-values) based on a Euclidean distance matrix and 10000 permutations. Pairwise comparison was conducted for following fixed factors and their interactions: F0-sex (grandparental (F0-Bi), grand-maternal (F0-Mat), grand-paternal (F0-Pat), grandparental control (F0-N)) and F2-bacteria (F2-bacteria control (F2-N), F2-bacteria Vibrio (F2-V+) and F2-bacteria Tenacibaculum (F2-T+)).(DOCX 56 kb)