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.

Experimental design and fish sampling

Syngnathus typhle individuals were caught in the seagrass meadow in Falckenstein (54.394497, 10.188988) in April 2020 (prior 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⁸ Vibrio aestuarianus 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_CPat_C), only female vaccinated (Mat_VPat_C), only male vaccinated (Mat_CPat_V), and both female and male vaccinated (Mat_VPat_V). 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_CPat_C: 5, Mat_VPat_C: 7, Mat_CPat_V: 5, Mat_VPat_V: 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 Artemia salina 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₁₈”), while the other two were kept at 25°C (heat stress situation; “Temp₂₅”; temperature treatment = “Temp”). Finally, juveniles were subjected to a bacterial exposure (“BE”): 50µl of a solution containing live Vibrio aestuarianus (again strain I11M2; 10⁸ Vibrio bacteria/ml) was added to 50ml rinsed one day-old Artemia 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_V”), while the other received the same amount of non-enriched Artemia nauplii (“BE_C”). 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 S. typhle 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).

Gut dissections

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.

RNA extraction, library synthesis, and sequencing

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 >30mio reads, despite some duplication levels being determined to be rather high, which might be reflected in some data noise.

Transcriptome assembly

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 S. typhle gene added an orthologous gene ID and gene description from Danio rerio, Hippocampus erectus, Syngnathus acus, and Homo sapiens. 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: Latimeria chalumnae, Danio rerio, Xenopus tropicalis, Hippocampus comes, Hippocampus erectus, Syngnathus acus, Syngnathus typhle (this annotation), Homo sapiens, Mus musculus. 26,072 genes were processed and orthology was identified for Homo sapiens in 14,435 genes, for D. rerio in 16,383 genes, for S. acus in 18,197 genes, and for Hippocampus erectus in 19,055.

(Pre-)analysis of the gene expression data set and removal of dissection-related variation

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>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<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 Danio rerio 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.

Pairwise comparisons of gene expression

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_VPat_CBE_CTemp₁₈, Mat_CPat_VBE_CTemp₁₈, Mat_CPat_CBE_VTemp₁₈, Mat_CPat_CBE_CTemp₂₅, Mat_VPat_CBE_VTemp₁₈, Mat_CPat_VBE_VTemp₁₈ and Mat_VPat_VBE_VTemp₁₈) against a reference group (Mat_CPat_CBE_CTemp₁₈) 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.

Average immune cell expression rank analyses

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.

Microbiome sequencing, assembly and analysis

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)<=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 Vibrio 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 Vibrio 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 Vibrio 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 Vibrio 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 Vibrio may affect the obtained microbiome sequencing signal of the respective treatment group (BE_V; i.e., if a large proportion of read depth was claimed by the experimental Vibrio, relatively fewer reads might have represented remaining Vibrio and non-Vibrio strains). Using Wilcoxon tests, we tested whether BE_V samples showed a higher read count than BE_C samples when all Vibrio strains were summed up per sample. Then, the same was tested after subtracting the experimental Vibrio counts in both groups. Subsequently, the sum of counts from samples that were exposed to the experimental bacteria (BE_V) was divided by the sum of all samples (BE_C+ BE_V) to estimate Vibrio treatment count bias in a proportion (the result indicated which proportion of reads came from BE_V 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 Vibrio was determined. We excluded genera that had counts in less than 5% of the 307 individuals considered, and we excluded the genus Vibrio, leading to 226 genera considered. Finally, a one-sample Wilcoxon test was used to test if the genus Vibrio was showing a significantly higher treatment bias than other genera.

For further analyses of the remaining microbiome, first, counts of the experimental Vibrio 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_C/VPat_C/VBE_CTemp₁₈), and one on samples with this exposure (Mat_C/VPat_C/VBE_VTemp₁₈; 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 Vibrio 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 Vibrio 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), Vibrio counts were suggested as a slightly better predictor than BE and it was therefore evaluated using the ANOVA (…,type=”III”) function.