DOI: 10.5061/dryad.5qfttdzgq
title: Data for: Parent-specific transgenerational immune priming enhances offspring defense – unless heat stress negates it all

description: Trans-generational immune priming (TGIP) adjusts offspring’s immune
 responses based on parental immunological experiences. It is predicted to
 be adaptive when parent-offspring environmental conditions match, while
 mismatches negate those advantages, rendering TGIP potentially costly. We
 tested these cost-benefit dynamics in the pipefish Syngnathus typhle
 (Syngnathidae). Because of their unique male pregnancy, egg production,
 and rearing occur in different sexes, providing both parents multiple
 avenues for TGIP. Parental bacteria exposure in our pipefish was simulated
 through vaccinations with heat-killed Vibrio aestuarianus before mating
 the fish to each other or to controls. The resulting offspring were
 exposed to V. aestuarianus in control or heat-stress environments, after
 which transcriptome and microbiome compositions were investigated.
 Transcriptomic TGIP effects were only observed in Vibrio-exposed offspring
 at control temperatures, arguing for low costs of TGIP in non-matching
 microbiota environments. Transcriptomic phenotypes elicited by maternal
 and paternal TGIP had limited overlap and were not additive.
 Parentally-induced transcriptomic responses were associated with immune
 functions, and specifically the paternal response to the innate immune
 branch, possibly hinting at trained immunity. TGIP of both parents reduced
 the relative abundance of the experimental Vibrio in exposed offspring,
 showcasing its ecological benefits. Despite TGIP’s significance in
 matching biotic environments, no TGIP-associated phenotypes were observed
 for heat-treated offspring, illustrating its limitations. Heat spikes
 caused by climate change thus threaten TGIP benefits, potentially
 increasing susceptibility to emerging marine diseases. We demonstrate the
 urgent need to understand how animals cope with climate-induced changes in
 microbial assemblages to assess their vulnerability in light of climate
 change.
description: <em>Experimental design and fish
 sampling</em> <em>Syngnathus
 typhle</em> individuals were caught in the seagrass meadow in
 Falckenstein (54.394497, 10.188988) in April 2020
 (<em>prior</em> to mating season), brought to the Helmholtz
 Centre for Ocean Research Kiel (GEOMAR) and kept in sex-specific 200L
 barrels with Baltic Sea water flow through. The temperature was gradually
 increased from 10°C to 18 °C and a day:night rhythm of 16h:8h was achieved
 over three weeks to approximate summer conditions in the Baltic Sea,
 inducing fish’s mating readiness. Animals were transferred to 80l aquaria
 attached to a common water circulation system and kept in same-sex groups
 (Fig 1A). To induce immune priming in offspring, first, half of the males
 and females were vaccinated subcutaneously with 50 µl heat-killed
 bacterial cells (10<sup>8</sup> <em>Vibrio</em>
 <em>aestuarianus</em> bacteria/ml strain “I11Ma2”; originally
 isolated from a pipefish foregut caught in the Mediterranean Sea, Italy
 (Roth et al., 2011)) twice, with one week in between injections. One day
 after the second vaccination, eight mating pairs (“families”) were
 established for each of the intended parental mating groups in separate
 tanks attached to the same water flow-through system (32 pairs in total):
 neither female (maternal treatment; “Mat”) nor male vaccinated (paternal
 treatment; “Pat”; group
 Mat<sub>C</sub>Pat<sub>C</sub>), only female
 vaccinated (Mat<sub>V</sub>Pat<sub>C</sub>), only
 male vaccinated (Mat<sub>C</sub>Pat<sub>V</sub>),
 and both female and male vaccinated
 (Mat<sub>V</sub>Pat<sub>V</sub>). The pairs were
 given two days for mating, after which females were removed, resulting in
 23 families with sufficient offspring numbers after male pregnancy for the
 following experiments (final family numbers were
 Mat<sub>C</sub>Pat<sub>C</sub>: 5,
 Mat<sub>V</sub>Pat<sub>C</sub>: 7,
 Mat<sub>C</sub>Pat<sub>V</sub>: 5,
 Mat<sub>V</sub>Pat<sub>V</sub>: 6). Body length
 and weight were recorded for all parents. Immediately
 after parturition, the offspring were transferred to a 1l tank per
 parental family in a flow-through system at 18°C, with the same day and
 night cycle as before, and fed twice a day with live <em>Artemia
 salina</em> nauplii (Fig 1A). Nine to ten days after birth, 20
 offspring individuals per family were randomly selected for further
 treatments and, as offspring of different families varied in age (matings
 and parturitions were not perfectly synchronized across mating pairs),
 subsequent experiments were performed on three consecutive days to match
 offspring age. As family was treated as a random factor in our statistical
 analyses, a potential effect of the sampling day is also accounted for as
 all individuals per family were processed on the same day. The 20
 individuals per family were split into four groups of equal sizes and each
 group was moved into a beaker with 300 ml Baltic seawater, which was
 aerated by a weak air bubbler (Fig 1B). Two beakers per family were kept
 at 18°C (control temperature; “Temp<sub>18</sub>”), while the
 other two were kept at 25°C (heat stress situation;
 “Temp<sub>25</sub>”; temperature treatment = “Temp”). Finally,
 juveniles were subjected to a bacterial exposure (“BE”): 50µl of a
 solution containing live <em>Vibrio</em>
 <em>aestuarianus</em> (again strain I11M2;
 10<sup>8</sup> <em>Vibrio</em> bacteria/ml) was
 added to 50ml rinsed one day-old <em>Artemia</em> nauplii in
 suspension. After an incubation period of 30min, one milliliter of this
 mixture was added always to one of the two beakers per family and
 temperature (“BE<sub>V</sub>”), while the other received the
 same amount of non-enriched <em>Artemia</em> nauplii
 (“BE<sub>C</sub>”). After one day in the beakers, all
 juveniles were killed via a lethal dosage of MS222 (0.04% in PBS), an
 overview photo was taken of all 5 juveniles per beaker to estimate
 treatment batch average body lengths using the segmented line tool in
 ImageJ (v.1.52k). Subsequently, juveniles were moved to RNAlater, kept for
 1d at 4°C, and then stored at -20°C. Juveniles were not monitored for a
 longer period of time as the transcriptomic response of key immune genes
 to infection is most pronounced 6-48h after an infection/immune
 stimulation in these fish (unpublished data). Additionally, the “natural”
 mortality (i.e., under the most ideal lab conditions that could be
 provided) and heterogeneity in growth rates of <em>S.
 typhle</em> juveniles massively increase after approximately 3 weeks
 of age, which is when a food transition occurs in captivity. This high
 mortality and growth rate heterogeneity would have confounded any estimate
 of survivorship in juveniles and would have substantially increased the
 required number of experimental animals. Unfortunately, keeping and
 treating five juveniles within a single beaker induced a nesting factor we
 cannot correct for (affecting the microbiota dataset exclusively; see
 below), however, individual-specific beakers were not feasible considering
 the large sample size. We therefore continued the analysis as if this
 effect would not exist, but acknowledge that conclusions from the
 microbiota dataset have to be drawn cautiously. All experiments were
 conducted in accordance with local ethics regulations (University of Kiel
 Antrag §7 V242-35168/2018). <em>Gut
 dissections</em> To investigate the gut
 microbiome of juveniles, a transversal cut through the neck of each
 RNAlater-stored specimen was conducted, still leaving the ventral body
 half connected to the trunk. Using forceps, the head (and attached
 tissues) was then pulled from the trunk, which led the ventral body wall
 to rip off directly posterior to the head, while the intestines remained
 connected to the head and instead got dislodged from the trunk close to
 the vent. Subsequently, the majority of the intestine was dissected from
 the head part, leaving only the most anterior part of the gut and
 neighbouring organs, such as the liver, with the head. Intestines were
 stored in 100% EtOH for microbiome analyses, while all remaining tissues
 were immediately used for RNA extraction. <em>RNA
 extraction, library synthesis, and sequencing</em>
 Total RNA was extracted from one individual per family and
 treatment batch (i.e., “beaker”), leading to a total of 92 samples (23
 families, with one individual for each of the two Temp and BE groups). For
 this purpose, remaining juvenile tissues (i.e., everything but intestines,
 which still includes immunogenic tissues, such as the head kidney, which
 are expected to respond strongly to an immune challenge) were taken
 immediately after gut dissection to homogenization in lysis buffer and
 processed using the Qiagen DNA/RNA All-Prep Blood and Tissue Mini kit
 (Cat. No. 80284). Extracted RNA had an average concentration of 41ng/µl
 and an average RIN of 9.7. mRNA library synthesis and sequencing (DNBseq,
 150bp PE, stranded) was conducted by BGI (Hong Kong, China), and on
 average 58mio reads/sample were obtained, with all samples having
 &gt;30mio reads, despite some duplication levels being determined to
 be rather high, which might be reflected in some data noise.
 <em>Transcriptome assembly</em>
 The RNA reads were quality-trimmed with Fastp (v0.20.1 (Chen et
 al., 2018)). Alignments to the annotated genome were carried out with STAR
 (v2.7.9a (Dobin et al., 2013)) in a two-pass mode with the following
 options “--outFilterIntronMotifs RemoveNoncanonical --outSAMunmapped None
 --outFilterMultimapNmax 1”. The count information per sample was obtained
 using TPMCalculator with options “-c 150 -p -q 255 -e -a” and then merged
 into a single multi-sample file with a script provided by the developer
 (tpmcalculator2matrixes.py). The output files from TPMcalculator were
 processed using a custom Python script that for each <em>S.
 typhle</em> gene added an orthologous gene ID and gene description
 from <em>Danio rerio</em>, <em>Hippocampus
 erectus</em>, <em>Syngnathus acus,</em> and
 <em>Homo sapiens</em>. Missing orthology information was
 marked as NA. The orthology assignment was performed with Orthofinder
 (v2.4.0 (Emms & Kelly, 2019)) using protein sequences from the
 following species: <em>Latimeria chalumnae</em>,
 <em>Danio rerio</em>, <em>Xenopus tropicalis</em>,
 <em>Hippocampus comes</em>, <em>Hippocampus
 erectus</em>, <em>Syngnathus acus</em>,
 <em>Syngnathus typhle</em> (this annotation), <em>Homo
 sapiens</em>, <em>Mus musculus</em>. 26,072 genes were
 processed and orthology was identified for <em>Homo
 sapiens</em> in 14,435 genes, for <em>D. rerio</em> in
 16,383 genes, for <em>S. acus</em> in 18,197 genes, and for
 <em>Hippocampus erectus</em> in 19,055.
 <em>(Pre-)analysis of the gene expression data set and
 removal of dissection-related variation</em> A
 TPM (transcripts per million) expression table of 92 juvenile samples was
 imported into R (v. 4.1.3,(R Core Team, 2021)) for downstream analyses and
 genes that were expressed in fewer than three samples were discarded,
 leaving 19,055 genes. Density plots of gene expression profiles were
 generated for each individual (Fig 2), revealing rather heterogeneous
 expression patterns, which may be indicative of non-homologous tissue
 sampling. To explore the main axes of variation in the
 gene expression data set, a Principal Component Analysis (PCA) was used on
 the log(gene_expr+1) transformed, centered, and scaled gene expression
 data set (Fig 3). When investigating the scores and loadings of PC1, we
 noticed that none of the major treatments seemed to be reflected on PC1
 scores (Fig 3A-C), using Dunnett’s Test with the non-primed group as
 control as a reference level to test for differences among parental
 treatments, and t-tests to test for differences between Temp and BE levels
 (all p&gt;0.1; package “DescTools “, v. 0.99.47;(Signorell et al.,
 2023)). Additionally, genes heavily negatively loaded on this PC
 unexpectedly appeared to be associated with energy metabolism (and, to a
 lesser degree, muscle tissue; Fig 3D), and a histogram of the loadings
 suggested an atypical negative skew, indicating that PC1 might reflect
 unintended variation present in the data set (Fig 3E).
 A set of possible technical and biological covariates was
 investigated for correlations with PC1, which would suggest an influence
 on overall gene expression patterns (Fig 4). These included per sample (i)
 the total RNA concentration after extraction and (ii) after shipping (as
 determined by the sequencing company), (iii) the respective juvenile’s
 father, and (iv) mother's total body length, (v) the average body
 length of juveniles of the respective batch (i.e., “beaker”), and (vi) the
 median of gene expression across all genes (which reflects the variation
 in gene expression profiles mentioned above, Fig 2). These variables and
 PC1 were correlated to each other using Spearman Rank-sum correlations and
 obtained p-values were corrected using the false-discovery-rate (FDR)
 method. The results revealed that PC1 scores were highly significantly
 associated with the median of the log gene expression (ρ=0.99,
 p.adj&lt;0.001), indicating that the variation in gene expression
 profiles displayed in Fig 3 indeed is the conceptually unintended major
 source of variation in the data-set. To understand the cause of this
 variation, PC1 loadings were inspected. The large majority of genes
 (16,979) showed positive loading, while much fewer (2,076) showed negative
 loading (Fig 3E). The 1% of most positively loaded genes (n=191) were
 analyzed for Gene Ontology (GO) biological process enrichment. The human
 GO annotation was used as it is the most comprehensive, even though not
 for all genes human ortholog annotations were present (n=165 with a human
 gene annotation; n=157 of those with GO annotation) and their function in
 human (and thus GO association) may be different compared to fish
 (corresponding results for GO analysis using <em>Danio
 rerio</em> annotations are reported in Data files uploaded to Dryad,
 but remained less conclusive and thus is not discussed here). We tested
 for GO enrichment using the gost(…,ordered_query=F,
 organism="hsapiens", significant=T,
 correction_method="gSCS", sources=c("GO:BP"),
 custom_bg=all_expressed_genes, domain_scope="custom") function
 (“gprofiler2” v.0.2.1; (Kolberg et al., 2020); all_expressed_genes include
 14,435 human orthologs). Eight terms were found significantly enriched, of
 which four were cellular regulation/organization related, while the
 remaining were generic regulatory processes, such as chromatin
 organization. When this was done with the 1% most negatively loaded genes
 (n=132 of which had a human ortholog and n=121 of these a GO annotation),
 twenty-one processes were identified and almost all appeared to be related
 to energy metabolism. These GO terms, the overall skew in gene loading
 frequencies towards positive loadings (Fig 3E), and the common energy
 metabolism genes with negative loading (Fig 3D; also, several
 muscle-related GO terms narrowly missed significance) suggest that
 variation reflected in PC1 is likely linked to differences in tissue
 composition (possibly reflecting a varying amount of the functionally
 distinct foregut portion of the gut) of samples. This is likely also
 reflected in the heterogeneous expression profiles (Fig 2), possibly as a
 result of inconsistencies during gut dissections. To
 correct for this effect and to identify whether and how the intended
 treatment factors affected gene expression, an unconventional analysis
 approach was chosen over more orthodox pipelines. The aim was to remove
 unintended variation in gene expression reflected in this PC1 in our
 subsequent statistical analysis. The applied pipeline thus proceeded as
 follows: individual genes were investigated and a gene was removed from
 the pipeline if it was not expressed in at least all but one individual of
 at least two treatment groups (from among the 16 treatment combinations),
 and mean log(gene expression+1) levels among individuals expressing the
 gene had to be at least 0.1 to be deemed meaningful. In this manner,
 1,567, 7, and 832 genes were excluded for being too low expressed, being
 expressed in too few samples, or both, respectively, leaving 16,649 to be
 considered in the following steps. Then, for each gene, a linear model was
 computed using the aforementioned log-transformed gene expression as the
 dependent variable and PC1 scores as the independent variable, and its
 residuals were stored. <em>Pairwise comparisons
 of gene expression</em> Gene residuals were used
 as the dependent variable to compute gene-wise linear mixed models that
 also included treatment variables and their interactions as independent
 terms (the model in R: lme(log_gene_expression_residuals ~ Mat * Pat * BE
 * Temp, random = ~1|family, method="REML"); package “nlme”, v.
 3.1-161; (Pinheiro et al., 2013; Pinheiro & Bates, 2006)). Due to
 their large number, individual models were not inspected for violations of
 model assumptions. For subsequent pairwise comparisons the function glht()
 (package “multcomp”, v.1.4-19;(Hothorn et al., 2016)) was used to test
 seven treatment groups
 (Mat<sub>V</sub>Pat<sub>C</sub>BE<sub>C</sub>Temp<sub>18</sub>, Mat<sub>C</sub>Pat<sub>V</sub>BE<sub>C</sub>Temp<sub>18</sub>, Mat<sub>C</sub>Pat<sub>C</sub>BE<sub>V</sub>Temp<sub>18</sub>, Mat<sub>C</sub>Pat<sub>C</sub>BE<sub>C</sub>Temp<sub>25</sub>, Mat<sub>V</sub>Pat<sub>C</sub>BE<sub>V</sub>Temp<sub>18</sub>, Mat<sub>C</sub>Pat<sub>V</sub>BE<sub>V</sub>Temp<sub>18</sub> and Mat<sub>V</sub>Pat<sub>V</sub>BE<sub>V</sub>Temp<sub>18</sub>) against a reference group (Mat<sub>C</sub>Pat<sub>C</sub>BE<sub>C</sub>Temp<sub>18</sub>) based on the mixed model. Estimates and raw p-values were computed for each pairwise comparison. Also, to compute corresponding comparisons with the reference temperature of 25°C, the aforementioned procedure was repeated with the reference level of the Temp factor being set to 25°C. Raw p-values were corrected for multiple testing using the p.adjust() function and the “fdr” method per comparison across all genes. For each comparison that yielded genes with FDR-corrected p-values below α=0.05, Gene Ontology (GO) terms were assigned to all genes when available and tested for enrichment as described before. <em>Average immune cell expression rank analyses</em> To explore if sets of significantly expressed genes were associated with immune functions, we investigated their association with immune cells. First, a reference single-cell expression data set based on human data was downloaded from the human protein atlas website using the R package “HPAanalyze” (v. 1.16.0;(Tran et al., 2019)). The “RNA single cell type” data set was chosen, which contains normalized expression data for 79 different annotated cell types (incl. “undifferentiated”) collated from various studies’ data sets. Of these cells, 13 were identified as being immune cell-related, and while “dendritic cells” were not assigned to one of the two branches of the immune system, those specific to the innate branch were "granulocytes", "Hofbauer cells", "Kupffer cells", "Langerhans cells", "macrophages", “microglial cells” and "monocytes", while those specific to the adaptive branch were "B-cells", "NK-cells" (natural killer cells), "Plasma cells", "T-cells" and “thymic epithelial cells”. All 14,435 human orthologs from our study were matched by Ensembl ID to the genes within the downloaded database. For each of these genes, we aimed to estimate its association in expression (in humans) to immune cell types, i.e., a high association reflects that a gene is predominantly expressed across immune cell types and rarely or only at low levels in other cell types. For this, first, cell types within the downloaded data set with no gene expression associated with them were excluded and the remaining cells were sorted according to their average normalized expression values (in humans) and ranks were assigned (the cell type with the highest expression value obtained the highest rank). Obtained rank numbers were then divided by the number of non-zero cells and the rank sum of the designated immune cells (IM value), those assigned to the innate (IN value), and adaptive immune system (AD value) were calculated, resulting in three values that reflect how strong a given gene’s expression is associated with cells of the respective type relative to all others (higher values indicate stronger association). Additionally, to illustrate if a gene is associated rather with cell types of one branch of the immune system over the other, the rank sums of the adaptive cell types were subtracted from the rank sum of the innate cell types, resulting in a value (IA) in which more positive numbers indicate higher gene expression association with cell types of the innate immune system, while more negative values indicate this for cell types of the adaptive immune system. While these values by themselves are not very informative, their distribution across all genes of the background gene set served as a reference to investigate the distribution of these values in sets of target genes, e.g., sets of differentially expressed genes. To determine if target gene sets contained genes that on average had different IM, IN, AD, or IA values compared to the reference background, Wilcoxon signed rank tests were performed and the distribution of ranks was plotted for both backgrounds and target gene sets for IM and IA values. Still, for these analyses it should be kept in mind that gene sets are not independent of each other, many genes overlap among them and p-values derived from Wilcoxon Sign-Rank tests were not corrected. Also, not all genes had a human annotation, and thus actual sample sizes for comparisons were somewhat smaller than DEG lists suggest. <em>Microbiome sequencing, assembly and analysis</em> 317 gut samples were processed using the Qiagen DNA Blood and Tissue kit (Cat. No.: 69504). Obtained DNA was sent to the IKMB sequencing center, Kiel, Germany, which performed library preparation and Illumina MiSeq amplicon paired-end sequencing of the V3-V4 bacterial 16S rRNA region randomly distributed to three runs (run 54: n=30, 55: n=8 and 57: n=269). Sequences were demultiplexed, conjoined with their respective metadata, and processed, quality-filtered, and analyzed (package “QIIME2”, v.2021.8 (Bolyen et al., 2019)). Correction of the fastq files was achieved using the DADA2 software within the QIIME2 pipeline, which included denoising, merging, chimera slayer fitting, and trimming of the sequences (“ASVs”). Due to declining quality scores, sequences were truncated 50 bases from forward and 70 bases (run 54) or 100 bases (run 55, run 57), respectively, from reverse reads. Taxonomy was assigned using a Naïve Bayes classifier trained on SILVA138 99% OTUs full-length sequence database. Microbiota count and metadata of all 317 samples were imported into R using the function qza_to_phyloseq() (package “qiime2R” v.0.99.6 (Bisanz, 2018); 89 samples of these were also used for gene expression analysis), with reads being assigned to 7,100 microbial taxa (ASVs). Samples with log(counts+1)&lt;=8 (~3000 reads) were excluded after histogram inspection (leaving 307 samples; minimum sample size per subgroup was n=13 and maximum number n=23). The remaining samples’ counts were normalized by dividing them by the respective sample’s mean count across taxa and then multiplying them by the mean count across all samples and taxa. Thus, when the text refers to “counts” for brevity, “normalized counts” or “relative abundance” is meant subsequently. The mean log(count+1) of all <em>Vibrio</em> ASVs was visualized using a heatmap (Fig 5). One strain (formerly two ASVs of identical barcode sequences were merged into one) showed not only the highest overall abundance among ASVs (16.5% of all microbiome counts) but also featured a marked difference in abundance between bacterial exposure groups. indicating that it represents the experimental strain. The experimental strains abundance across samples that also had RNA-seq data (n=87) did not correlate with the bias-indicating original RNA-seq PCA PC1, suggesting that no correction was required (S = 127812, p-value = 0.1273). Treatment effects on the experimental <em>Vibrio</em> counts were analyzed by first log-transforming all 307 samples (log(counts+1)), followed by visual inspection across treatment groups nine outliers were excluded. The remaining samples’ transformed <em>Vibrio</em> counts were used as the dependent variable in a linear mixed model (lme() function in R, with maximum likelihood method). A QQ-plot for residuals was used to investigate the normality of residuals and some deviation from normally distributed residuals was noted, mostly due to eight outliers. Attempts to remove these outliers led to model results very similar to those including outliers (see below), except for stronger and more significant effects, illustrating their limited effect on the conclusions drawn from the model. Outliers were therefore retained. As in the gene expression analyses, independent variables were the factors Mat, Pat, BE, and Temp, and all their interactions. A random intercept for the family factor was included (“random=list(~1|family)”), but also a variance structure for the BE-treatment (“weights=varIdent(form=~1|BE)”), compensating for the heterogeneous variance due to this treatment. As the focus of the analysis was to determine if parental priming induced by parental vaccinations affects microbe counts, the <em>Vibrio</em> exposure treatment level (and not the control) was chosen as a reference level for the BE factor. The Anova() function was used again to obtain an ANOVA table, and type III ANOVA was chosen as interactions rather than main effects were of central interest. Adding large amounts of the experimental <em>Vibrio</em> may affect the obtained microbiome sequencing signal of the respective treatment group (BE<sub>V</sub>; i.e., if a large proportion of read depth was claimed by the experimental <em>Vibrio</em>, relatively fewer reads might have represented remaining <em>Vibrio</em> and <em>non-Vibrio </em>strains). Using Wilcoxon tests, we tested whether BE<sub>V</sub> samples showed a higher read count than BE<sub>C</sub> samples when all <em>Vibrio</em> strains were summed up per sample. Then, the same was tested after subtracting the experimental <em>Vibrio</em> counts in both groups. Subsequently, the sum of counts from samples that were exposed to the experimental bacteria (BE<sub>V</sub>) was divided by the sum of all samples (BE<sub>C</sub>+ BE<sub>V</sub>) to estimate <em>Vibrio</em> treatment count bias in a proportion (the result indicated which proportion of reads came from BE<sub>V</sub> samples). To understand if the obtained proportion was similar across microbiota genera (suggesting sequencing depth was a limiting factor across samples), for each sample, the sum of all counts of strains belonging to each genus identified in this study was calculated (853 genera). Then, for each genus, the proportion of counts found in samples exposed to the experimental <em>Vibrio</em> was determined. We excluded genera that had counts in less than 5% of the 307 individuals considered, and we excluded the genus <em>Vibrio</em>, leading to 226 genera considered. Finally, a one-sample Wilcoxon test was used to test if the genus <em>Vibrio</em> was showing a significantly higher treatment bias than other genera. For further analyses of the remaining microbiome, first, counts of the experimental <em>Vibrio</em> strain were removed and samples re-standardized. Then, all counts from within each microbial family taxon were summed up per individual, and families with a total of fewer than ten counts were removed (332 microbial families remaining). To test if the immune-challenge phenotypes elicited by parental immune priming also affected this remaining microbiome, two ANOSIM analyses were conducted using Bray-Curtis dissimilarity matrices: one on microbiome abundances from samples without bacterial challenge at 18°C (Mat<sub>C/V</sub>Pat<sub>C/V</sub>BE<sub>C</sub>Temp<sub>18</sub>), and one on samples with this exposure (Mat<sub>C/V</sub>Pat<sub>C/V</sub>BE<sub>V</sub>Temp<sub>18</sub>; anosim(); “vegan” package v.2.6-2), for both of which the factor level combinations of Mat and Pat were treated as one factor with four levels. The Shannon Diversity index was calculated per individual using log(counts+1) and its correlation with experimental <em>Vibrio</em> counts was assessed using a Spearman rank correlation. Additionally, to estimate how treatment groups affected alpha diversity, two linear mixed models were computed using the alpha diversity as the response variable and Mat, Pat, and Temp as predictor variables, plus either BE or the experimental <em>Vibrio</em> log(counts+1), as well as all interactions (“family” was added as a random intercept). For both models, model selection using the stepAIC() function was conducted, and while both models yielded similar BIC values (200.97 and 199.49, respectively), <em>Vibrio</em> counts were suggested as a slightly better predictor than BE and it was therefore evaluated using the ANOVA (…,type=”III”) function.
description: # Data for: Parent-specific transgenerational immune priming enhances
 offspring defense – unless heat stress negates it all
 [https://doi.org/10.5061/dryad.5qfttdzgq](https://doi.org/10.5061/dryad.5qfttdzgq) ## Description of the data and file structure This supplemental data is associated with the article "Parent-specific transgenerational immune priming enhances offspring defense – unless heat stress negates it all" published in Ecology and Evolution. Raw transcriptome and microbiome read data was uploaded on NCBI, and the associated SRA number can be found in the manuscript's Data Accessibility statement. Adult *Syngnathus typhle* individuals were vaccinated twice with heat-killed *Vibrio* of an experimental strain. Vaccinated adults were then crossed with control adults or other vaccinated adults. Also, control crosses were produced and a fully reciprocal design was achieved, with either mother, father, both parents, or none of them having been vaccinated. Offspring from crosses were raised for ~10 days and then, from each family, five offspring individuals were treated for one day either in a controlled environment without further *Vibrio* exposure, in a controlled environment with experimental *Vibrio* exposure, in a heat stress environment without further *Vibrio* exposure, or in a heat stress environment with experimental *Vibrio* exposure. Subsequently, offspring were killed, guts were used for DNA extraction and subsequent microbiome analyses while the remaining bodies (including immunogenic organs) were used for RNA extraction and RNA-seq. ### Files and variables The supplemental data provided here contain the results for the preprocessing of transcriptomic data ("Data 1_transcriptomics_preprocessing.xlsx"), the results of statistical analyses (pairwise comparisons based on the mixed model) using this preprocessed data ("Data 2_transcriptomics_final_analysis.xlsx"), results of GO analyses (“Data 3_GO_analyses.xlsx”), results of the average immune cell rank association analysis (“Data 4_immune cell score analysis.xlsx”), and the count, meta and taxonomy data of the microbiota data set (“Data 5_microbiome_data.xlsx”). Treatment groups are abbreviated as "P" or "Pat" for paternal (usually vaccination), "M" or "Mat" for maternal (usually vaccination), "BE" for bacterial exposure, "T18" for the 18°C group, and "T25" for the 25°C group. For annotation, "DR" indicates Danio rerio and "HS" Homo sapiens. ## Code/software ".xlsx" files can be opened with LibreOffice or other free software able to handle the file format.