Wild fish were sampled from an oligotrophic upland lake (FRN, 52.3599, −3.8776) and lowland river (RHD, 52.4052, −4.0372) in mid Wales in late (astronomical) summer 2012 and late winter 2013 (8/site in summer and 10/site in winter) (A diagrammatic summary of the overall study design is provided in Figure S1). RNA was extracted from whole fish and a genome-wide quantitative gene expression profile generated for each fish by RNAseq, as described in (21). Gene expression data utilized here are expressed in fragments per kilobase of exon per million reads mapped (FPKM), as previously reported. We additionally extracted reads of bacterial 16S rRNA origin from quality controlled Fastq run files for individual fish. Reads were classified against the Greengenes database (operational taxonomic units, OTUs, clustered at 99%) using Taxonomer (27) (for further details, including of controls and identification of possible contaminants, see Supplementary Materials and Methods and Table S1). There were no missing values in the wild fish dataset.

We carried out an experiment simulating season-specific diet at FRN under otherwise common garden conditions (A diagrammatic summary of the overall study design is provided in Figure S1). The seasonal diet treatments were motivated by the finding of Allen and Wootton (6), and our own casual observation (pers. obs.), that during colder periods of the year at FRN, stickleback gut contents can include a substantial proportion of plant-derived detritus and algae. The ingestion of these items may be due to less discriminate benthic foraging during periods of the year when preferred food resources are limited.

To illustrate this trend we re-analyzed the year-round monthly stomach content data recorded at FRN by Allen & Wootton (6), according to the following considerations. From Figure 1 in Allen & Wootton (6) we extracted monthly relative diet composition data, using Plot digitizer 2.6.8. An index constructed from these data, of the proportion of stomach contents made up by detritus and plant material as opposed to animals (detritus index, DI), varied seasonally when analyzed by cosinor regression (28) (overall deletion test of sinusoid terms, p = 0.05; amplitude parameter, p = 0.004; η² = 46.1%, n = 12), with lowest predicted values (highest animal content) in June and highest predicted values in January (see Figure 1). Detritus and algae constituted 0–39% of monthly stomach contents.

Taking into account seasonal representation of invertebrates in the study of Allen & Wootton (6) and in other seasonal diet studies of freshwater stickleback in the U.K. (29, 30) we then formulated a summer-like diet treatment consisting of common animal prey items (copepods, cladocerans and chironomid larvae) and a winter-like diet treatment consisting of plant detritus (more typical of stomach contents during winter months) and chironomid larvae (a prey item typically consumed year-round). For the winter-like diet, boiled (1 h), triturated spinach leaves, representing plant detritus, were mixed with triturated gamma irradiated chironomid larvae (Tropical Marine Centre 8153) at a ratio of 3.8 g per 100 ml sieved detritus. The summer-like diet consisted of equal amounts of chironomid larvae, cyclopoid copepods (Tropical marine centre 8171) and cladocerans (Tropical marine centre 8151) (all gamma irradiated). Each diet was provided in excess in each tank and was continuously available ad libitum. We confirmed visually (on dissection at the end of the experiment, see below) that fish fed the winter diet ingested substantial amounts of plant detritus during foraging for chironomid fragments.

Fish destined for the experiment (carried out at Aberystwyth University) were captured at FRN in winter (January) and acclimatized to outdoors tanks (exposed to natural photoperiod and temperature) for 1 week, during which they were treated twice with praziquatel, as previously described (24), to remove gyrodactylid and diplostomatid infections. Fish were then exposed to the diet conditions, in the same outdoors tanks (see below), for 21 days. Throughout the experiment, fish were maintained in outdoors 10 L tanks (n = 8 fish per tank; 2 tanks per treatment group), exposed to natural photoperiod and temperature variation (4 ± 3°C), in conditioned tap water buffered to a PH of 7.5. Water was changed every day, with food replaced after each water change. Ammonia, nitrite and nitrate levels were monitored daily and remained at negligible levels.

At the end of the experiment fish were individually netted and immediately killed by concussion and decerebration and immersed in RNA stabilization solution. A small longitudinal slit was made in the ventrum of each animal prior to immersion to assist penetration of the stabilization solution to the body cavity. Storage of preserved samples was at −70°C. Weight (mg), body length (snout to tail fork, mm) and sex (based on gonad form and pigmentation) were recorded upon thawing prior to nucleic acids extraction (see below). For each fish, we generated a quantitative gene expression profile by quantitative real-time PCR (QPCR, see below).

The anti-microbial treatments applied to the diets (boiling or gamma-irradiation), and also the conduct of the tank experiment under non-sterile conditions, entail that microbes in the foods would not be completely eliminated but would be much reduced and greatly altered in composition. We would thus expect the dietary microbe-driven effect on host immunity, if any, to be very unlike that seen in the wild (where microbial exposures would be very different). Equally, an effect on host immunity through different microbial colonization of the (initially identical) tank waters due to the different foods, if any, would also likely be very different to any exposure in the wild. Moreover, the latter possibility would be limited by regular (24 h) water changes and the low environmental temperature (4 ± 3°C), curtailing microbial growth in the tank water. A similar argument would apply to differential toxicity of the experiment diets (i.e., if these were unrepresentative of the wild diets then different effects on immunity would be expected). As the immunophenotypic responses to diet treatment measured below were identical (rather than different) in form to those occurring in the wild under intentionally similar diet changes, we consider the above possibilities unlikely. We additionally note that a design that totally eliminates notional differential bacterial colonization of the food or environment is only possible in a fully microbe-free system, and that such a system (necessarily including captive-reared hosts) would be so artificial as to offer limited insight in the present case.

For experiment fish we focussed on a set of 5 genes (seasonal reporter, SR, genes) whose expression we have shown to report a major immunome-wide seasonal expression signature (21, 22). These include two genes expressed highly in summer: cd8a (Emsembl gene identifier: ENSGACG00000008945) and foxp3b (ENSGACG00000012777); and three genes expressed highly in winter: orai1 (ENSGACG00000011865), tbk1 (ENSGACG00000000607), and an IL-1 receptor genomic cluster member termed here il1r-like (ENSGACG00000001328). For wild fish we extracted expression estimates from RNAseq data and there were no missing values (see above). For experimental fish we measured gene expression in spleen, gill, liver and skin (fin) by QPCR using methods similar to those previously described (22, 24) (see Supplementary Materials and Methods). All 5 genes above were measured in all 4 tissues, with the exception of orai1, which was not measured in fin due to a technical omission. QPCR measurements considered below are relative expression values, normalized to two endogenous controls genes and indexed to a common calibrator sample by the ΔΔCt method (see Supplementary Materials and Methods). A limited number of missing values occurred for some tissue-specific expression variables and these are summarized in Table S2.

DNA was extracted from intestinal tracts preserved in RNA stabilization solution using the MoBio PowerSoil DNA Isolation Kit and samples (randomized) were passed to a sequencing service (Macrogen) for 16S rRNA sequencing (paired end) on a MiSeq machine (Illumina). DNA quality was assessed by electrophoresis and quantified using Quant-iT™ PicoGreen™ dsDNA Assay Kit by fluorimetry on a Victor³ reader (Perkin-Elmer). Libraries targetted the 16S V3 and V4 region (amplicon ~460 bp) with Illumina-recommended primers based on (31) (forward 5′ TCGTCGGCAGCGTCAGATGTGTATAAGAGACAGCCTACGGGNGGCWGCAG; reverse 5′ GTCTCGTGGGCTCGGAGATGTGTATAAGAGACAGGACTACHVGGGTATCTAATCC). Library size distributions were assessed with an Agilent Technologies 2100 Bioanalyzer using a DNA 1000 chip. Following sequencing, Fastq files supplied by the sequencing service were merged using the vsearch script -fastq mergepairs (32), then using usearch v9 (33) quality filtered (expected error < 0.5 removed) and truncated to 440 bp. Sequences were de-replicated, sorted by cluster size and OTUs clustered at 97% using the UPARSE algorithm (34). Suspected chimeric sequences were identified using UCHIME (35) against the fullSilva 97% formatteddatabase (36) and removed. Taxonomy was assigned against the full SILVA 97% OTUs reference database. OTU clusters with <5 reads, or of archaeal, eukaryotic, unassigned or chloroplast origin, were removed from the final OTU table. There were no missing values in this dataset (For further details, including of controls and identification of contaminants, see Supplementary Materials and Methods and Table S3).

Analyses were carried out in R Version 3.4.4. We employed permutational multivariate analysis of variance using distance matrices (PERMANOVA-DM), implemented in the vegan package, to partition variance in bacterial community composition (OTU relative abundance) between host and environmental variables. For consistency, we also used PERMANOVA-DM to partition variation in immunophenotypic measures between the same sets of explanatory variables. Analyses of wild site and experiment datasets were conducted separately. PERMANOVA-DM models of the wild site data reported below analyzed bacterial community composition or immune gene expression variables as the response and contained season (summer/winter), site (river/lake), sex (male/female), and body length (continuous) as explanatory terms. PERMANOVA-DM models of the experiment data reported below analyzed bacterial community composition or immune gene expression variables as the response and contained diet treatment (summer-like/winter-like), sex (male/female), and body length (continuous) as explanatory terms. We additionally considered body condition as an explanatory term in the wild site models, and body condition and tank as explanatory terms in the experiment models, but omitted these from final models as they were always non-significant. Where it was desirable to adjust for the effects of one variable that was colinear with another (diet and condition in the experiment), we limited permutations within strata of the adjusted variable. To consider the association between bacterial community composition and immunophenotypic profile we used Mantel correlations of distance matrices to assess the crude (unadjusted) association and partial Mantel correlations to assess the association adjusted for the effect of season or diet treatment (vegan package). Principal co-ordinates analysis (PCO; labdsv package) was used to ordinate both the bacterial and immunophenotypic datasets, visualizing the scatter of individuals against the major axes of variation. Manhattan distance was used to construct all distance matrices. Gene expression variables in all of the above analyses were first log transformed (log₁₀ x or log₁₀ [x + 1]). Post-hoc (to the PERMANOVA-DMs) gene-by-tissue confounder-adjusted analyses of season or diet effects on gene expression were carried out in general linear models (LMs) using the core lm command. Prior to these analyses, individual gene expression variables were subject to an optimal power transformation (via a Box-Cox procedure, MASS package) and then standardized (zero mean, unit standard deviation) for comparability of parameters.

Consistency of gene expression responses across different tissues (in fish from the experiment) was initially assessed by a pairwise correlation analysis (Pearson coefficients) and by ordination of tissue specific gene expression by PCO (see above). We additionally implemented a multivariate generalized linear mixed model (MGLMM) (37, 38) of multiple gene expression responses to diet in different tissues, with individual identity as a random term (MCMCglmm package). This allowed estimates of repeatability (across tissues) for individual genes and the estimation of a diet × tissue interaction term to assess the existence of a consistent (systemic) response (i.e., where interaction is important, this is evidence for a lack of consistency across tissues). Gene expression data were log-transformed and standardized for this analysis.

Comparing our two datasets (wild fish at different seasons and experimental fish with different season-specific diet) we also asked whether there was any consistency in the bacterial OTUs associating with season and diet treatment. In each datasest, we ran individual confounder adjusted permutation tests (based on LMs; lmPerm package) for each OTU (relative abundance) against season or season-specific diet treatment, identifying genes that were significant at a false discovery rate (FDR) cut-off of P = 0.05. We additionally ran random forests machine learning analyses (39) to predict season or diet treatment from the bacterial datasets, identifying variables with high importance.

Whilst our analyses of bacterial variation concentrated on community composition, we did secondarily consider rarefied species richness (in LMs with explanatory terms equivalent to those in the PERMANOVA-DMs used to analyse community composition). The only significant trends were modest associations with body length, which might possibly be due to the dynamics of colonization in wild fish.

Where body condition measures are considered below we analyzed both the residuals from a quadratic regression of body weight on length (RBW) (40) and the scaled mass index (SMI) (41) (SMI and RBW being highly correlated). Term-specific effect size measures presented below are partial R² for PERMANOVA-DM and classical η² for LMs.

The basic data from this study will be available in the European Nucleotide Archive (primary accession number PRJEB13319).

We note that inherently different nucleic acid extraction, sequencing and bioinformatics methodologies in our RNAseq and 16S amplicon datasets, likely leading to different biases, make inference from differences between these untenable (i.e., methodology is confounded with dataset). However, we a priori determined a strategy of only searching for, and deriving inference from, commonalities between the datasets. Such commonalities, emerging despite methodological differences, are likely to be more robust even than the case of a common methodology (where shared methodological biases might lead to systematic error).

We first asked whether wild sticklebacks differed in winter and summer. We found a modest effect of season on bacterial community composition [RNAseq data; PERMANOVA-DM, F_{(1, 33)} = 3.32, P < 0.001, R² = 8.7%], a weak effect of site [F_{(1, 33)} = 1.91, P = 0.012, R² = 5.0%] and no effect of sex or body length (Figure 2A). Season was a substantial source of variation in immunome-wide [RNAseq data; n = 3,648 genes; PERMANOVA-DM, F_{(1, 31)} = 6.29, P < 0.001, R² = 13.1%] and SR gene expression [RNAseq data; n = 5 selected seasonal immunome genes; F_{(1, 34)} = 63.79, P < 0.001, R² = 65.2%] (Figure 2A; see also Figure S2), as previously reported (21). Amongst the SR genes, orai1, tbk1, and il1r-like were expressed relatively more in winter (winter-biased) and cd8a and foxp3b were expressed more in summer (summer-biased), also as previously reported (21). Bacterial community composition was uncorrelated with immunome-wide expression but significantly correlated with SR gene expression (Mantel r = 0.21, P = 0.001), though this association disappeared when season was adjusted for. Body condition (RBW and SMI) did not vary between the winter and summer fish [see also (24)] and had no association with bacterial community composition or immune gene expression when added to the PERMANOVA-DM models above.

To assess the causal effect of diet in driving the seasonal changes seen in the wild, we carried out an experiment in which we provided season-specific diets under otherwise common garden conditions. For this experiment, we acclimatized wild fish (captured in January) to outdoors aquaria and applied the diet treatments under ambient temperature and photoperiod conditions. For representativeness, we monitored gene expression in spleen, liver, gill, and fin.

We found a dominant overall effect of diet on SR gene expression [QPCR data; PERMANOVA-DM, F_{(1, 15)} = 11.80, P < 0.001, R² = 44.8%] (Figures 2B, 3A,B; see also Figure S3) and diet was also the dominant source of variation in all individual tissues (compared to sex and body size) (Figure 2B). For all tissues the summer-like diet forced SR gene expression changes in the same direction as winter to summer changes in the wild (Figure 4). In post-hoc analysis of individual genes, those that have previously been established to be winter-biased in the wild (orai1, tbk1, and il1r-like) [see (21, 22)] tended to be down-regulated, whilst known summer-biased genes (cd8a and foxp3b) tended to be up-regulated, by the summer-like diet (Table 1). Stronger and more consistent diet effects were observed for T-cell-associated genes (especially cd8a and foxp3b), but significant trends, in the expected season-specific direction, were observed for all SR genes in one or more organs. There were no effects of sex or body length on immune expression (and no sex by diet interaction) (PERMANOVA-DM) and no effects of body condition (RBW and SMI; accounting for diet effects in a stratified PERMANOVA-DM).

There was a substantial element of individual consistency in the gene expression responses across different tissues, especially for those genes showing strong responses. Thus, T-cell associated genes (cd8a, foxp3, orai1) had large across-tissue repeatabilities (0.4–0.5), but other genes showed non-significant repeatability overall. Nonetheless, tissue-specific responses of all individual genes tended to be positively correlated amongst different tissues and to respond to diet treatment along a consistent multivariate trajectory (in a MGLMM, tissue × treatment interaction was significant overall, p = 0.005, but not a comparatively large effect) (Figure 5). Taken together, there was thus strong evidence of a prominent systemic trend in the immunophenotypic response to diet.

The winter-like diet maintained fish in similar body condition to wild fish in January and February (Figures 6A,B). In great contrast to the lack of winter-summer variation in condition seen in wild fish (see above), the summer-like diet promoted a substantial increase in body condition compared to both the winter-like diet (RBW, LM parameter = 48.0 ± 13.2 mg, tukey P = 0.002; SMI, LM parameter = 48.3 ± 10.3 mg, tukey P < 0.0005) and wild winter fish (RBW, LM parameter = 57.4 ± 12.9 mg, tukey P < 0.0005; SMI, LM parameter = 60.0 ± 10.0 mg, tukey P < 0.0005) (Figures 6A,B). There was a slight, but non-significant, increase in mean length in the summer-like diet group, compared to the winter-like diet group (which were size matched at the start of the experiment) (Figure 6C).

We measured bacterial community composition in the gut and gill in the above experiment (Figures 2B, 7). There was a substantial effect of diet on bacterial community composition in the gut [16S amplicon sequencing data; PERMANOVA-DM, F_{(1, 25)} = 4.66, P = 0.002, R² = 15.7%] and gill [F_{(1, 16)} = 6.87, P < 0.001, R² = 30.0%], but no effects of sex or body length (and no sex by diet interaction). Body condition (RBW and SMI) was not associated with bacterial community composition in gut or gill if the effect of diet treatment was accounted for (stratified PERMANOVA-DM). We also tested our prediction (above) that the correlation of microbiome composition with immune expression would be stronger in the gill than the gut. This prediction was confirmed when local microbiome composition was compared, in a planned comparison (unadjusted for diet effect), to systemic (splenic) SR gene expression: Mantel r for the gill was 0.64 (p = 0.001, 95% CI = 0.51–0.72, n = 28) and only 0.17 (p = 0.014, 95% CI = 0.09–0.27, n = 28) for the gut (Splenic expression was chosen for the above comparison because spleen was considered likeliest to reflect systemic responses and also showed the largest diet effect size). Adjustment for diet effect tended to negate the association between the microbiome and splenic SR gene expression, which became non-significant for the gut microbiome and only marginally significant for the gill microbiome.

There was no genus-level consistency in the bacterial lineages that best predicted season in the wild, compared to the lineages that best predicted season-specific diet treatment in our experiment (Figure 8). Strikingly, however, some of the OTUs most predictive of season and diet treatment were mycolic acid-producing Corynebacteriales that increased in abundance under summer conditions. More specifically, a Tsukamurella lineage was 4th most predictive of season out of 1,816 lineages (individually significant with FDR adjustment) (Figure 8A). At the same time a Corynebacterium lineage was by far the most predictive of diet treatment out of 3,398 gut lineages (Figure 8B), and the only lineage in which the diet effect was significant after FDR error rate adjustment. Post hoc analyses indicated that Corynebacteriales as a whole, and the frequently pathogenic corynebacterial genus Mycobacterium individually, increased significantly (p < 0.05) in relative abundance in the gill and gut of experimental fish fed the summer-specific diet (Figure 8C). Although the corresponding trends were not observed in wild fish, the relative abundance of a further corynebacterial genus, Gordonia, did increase significantly in summer. We note that, whilst less frequently pathogenic than Mycobacterium, Corynebacterium, Tsukamurella, and Gordonia have all sometimes been associated with disease in fish (42–44). We also note that a lineage of the potentially pathogenic firmicute genus Streptococcus was the 5th most predictive of the summer condition in wild fish (Figure 8A), increasing in relative abundance. Post hoc analysis revealed this genus also increased in both the gut and gill of experimental fish fed a summer-like diet (P < 0.05; Figure 8C).

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.