BACKGROUND
seasonal immune function has often been observed in vertebrates  <cit> , including humans  <cit> , but is relatively poorly understood.
as with more studied circadian rhythms, though, there are fundamental implications for health  <cit>  .
taking a comparative approach  <cit>  and considering the conserved genes of the vertebrate immune system, here we use transcriptomic measurements to reveal seasonal re-adjustments of immunity in a naturally-occurring teleost model.
crucially, through focussing on wild organisms exposed to real-world environmental extremes, we expected to discover more measurable variation than would be the case in a study of domesticated animals, where seasonal variation may be muted by anthropogenic influences.
we chose the three-spined stickleback  as a subject because it has an annotated whole genome sequence  <cit>  and occurs accessibly in highly seasonal natural habitats.
also, a considerable knowledge base exists for this species: it is a highly studied model organism  <cit> , and there are particularly detailed ecological studies relating to our main study area, mid wales  <cit> .
we compared the transcriptomes of populations in late winter and late summer  in ecologically divergent natural populations, reasoning that a focus on common responses would provide a way to filter overriding seasonal trends from locality-specific variation.
we chose to primarily use global mrna extracts from individual whole fishes rather than from isolated cell populations or tissues.
this was because a fully reductionist approach to cell populations would be impractical, and because the majority of the teleost immune system is likely to be diffusely distributed in the gut, under the skin and mucosal surfaces and in association with the gills and liver  .
by considering global  samples, we were thus able to take a holistic view of which immune system pathways are differentially expressed at seasonal extremes.
we corroborated our transcriptomic analyses by targeted gene expression measurements in year-round samples of fishes from the original sites in a new annual cycle, and by tissue-specific analyses at a further site.
moreover, using network analyses we were able to ask what genes are important in regulating seasonal immune function and how do seasonally-biased immune networks interact with other seasonally-biased organismal processes?
RESULTS
seasonal expression bias of immune system genes occurs against a well-defined genome-wide seasonal signature
we analyzed the global  transcriptomes of g. aculeatus from two contrasting habitats in mid wales, river  rheidol  and lake  frongoch , in september  <dig> and march  <dig>  to begin our analysis we considered, genome-wide, which genes were associated with seasonal expression bias.
at frn,  <dig> genes were significantly differentially expressed from summer to winter with an individual cut-off  and  <dig> with a false discovery rate -adjusted cut-off.
at rhd,  <dig> genes were significantly differentially expressed with an individual cut-off and  <dig> with an fdr-adjusted cut-off.
genes that were seasonally differentially expressed at both localities tended, overwhelmingly, to show synchronous expression .
we hypothesized that these synchronously differentially expressed genes would also be those contributing to overarching seasonal responses.
thus, we categorized such genes on the basis that they were significantly  differentially expressed, in the same direction, at both sites at the individual error rate .
following this criterion,  <dig> genes were differentially expressed in a consistent direction ,  <dig> increasing expression during winter  and  <dig> increasing expression during summer .
we then considered genome-wide winter-summer expression changes in functional gene sets using gene set enrichment analyses  , focussing on predicted orthologues of human genes .
we compared ranked gene expression to kegg and reactome pathways, globally .
we also compared the ranked genes to a selected group of gene sets reflecting different organismal processes and immunological pathways .
finally, we considered hypergeometric overlap of the selected gene sets with the summer- and winter-biased genes defined above .
taken together, these analyses  indicated signatures of growth in summer  and physiological challenge in winter .
in all analyses, some immune pathways were seasonally-biased, with innate processes  emphasised in winter and adaptive  processes in summer .fig.
<dig> distinctive immunological and genome-wide gene expression signatures occurred at seasonal extremes.
a gene sets with significant summer  or winter  expression bias as indicated by gene set enrichment analysis .
ranked differential gene expression was compared, separately for the rhd and frn sites, to global kegg and reactome gene sets, and sets are shown where the combined fdr p value was significant ; gene set names are truncated but shown in full in additional file 2: table s2; stars indicate immunologically-relevant gene sets; the central dendrogram indicates the degree of overlap between gene sets.
b analyses of selected gene sets  representing immunological pathways and organismal signatures of stress, reproduction, growth and metabolism.
individual colour panels correspond, left to right, to the order of gene sets in additional file 3: table s <dig>  these sets were considered by gsea for rhd and frn separately, and by gene overlap  for the overall summer and winter-biased gene sets
the above analyses were carried out on expression data un-adjusted for individual size, as this variable was  approximately balanced across winter and summer samples.
however, as our sampling points bounded a non-recruiting population ageing in the interval between breeding seasons, we considered in more detail the potential influence of ontogenetic stage.
it is likely, given the months  in which we recorded reproductive activity in the field, and taking into account slower biological ageing at lower temperatures , that the 0+ cohort in our approximately size-matched summer and winter samples would have included individuals widely overlapping in effective age .
furthermore, as our sampling deliberately selected a wide range of fish sizes, it is probable that 0+ and 1+ cohorts  <cit>  were represented, resulting in a very extensive overlap of effective ages between summer and winter samples.
a close association between body size and age allows age to be partitioned from season in statistical models by the use of a size metric, such as body length, as a surrogate.
this is validated by data from experiments in artificial outdoor habitats, where we found that time explains at least 57 % of the variation in length .
in order to control for age  effects we applied general linear models  to each in turn of the  <dig> genes in the gsea dataset, including main effects for season, length, sex and site.
we found that season was the dominant predictor of gene expression .
consistent with this gene-by-gene analysis, a multivariate principal co-ordinates analysis  of the same data demonstrated marked differentiation across seasons against the two major axes , but none for length, sex or site .
we also re-ran “global” gsea analyses , first with genes ranked by confounder-adjusted seasonal effect, and then with genes ranked by confounder-adjusted length effects .
we found a similar outcome in the analysis ranked by confounder-adjusted seasonal effect to in the unadjusted analysis shown in fig.
1a, with the two analyses sharing  <dig> gene sets that were significantly seasonally enriched , including all of the immunological sets except leucocyte transendothelial migration .
in contrast, there was a distinctive outcome in the analysis ranked by confounder-adjusted length effect, where only  <dig> enriched gene sets were shared with the analysis shown in fig.
1a, including none of the immunological sets.
thus, the effect of season was a very dominant one, emerging clearly in analyses even without adjustment for ontogeny.
season is a dominant and consistent influence on immune gene expression
focussing on the immune system we considered genes with predicted orthology to those in the immport  <cit>  comprehensive list of immune-associated genes.
in total,  <dig> immune-associated genes out of  <dig> were seasonally-biased .
of these, greater absolute numbers of genes were winter-biased  than summer-biased .
as with the full gene list, immune-associated genes that were consistently seasonally differentially-expressed tended, very strongly, to be synchronously expressed across localities .fig.
<dig> extremes of season were the dominant predictor of immune gene expression.
a scatterplot of log <dig> winter-summer fold expression change  for all immune-associated  genes with significant seasonal difference  at both rhd and frn sites.
overwhelmingly such genes were regulated in the same direction across sites.
b season was overwhelmingly the most important predictor of immune gene expression, in comparison to site, sex and body size .
bar chart summarizes results from general linear models  fitted to each individual log <dig> immune-associated gene expression variable; bars are the mean observed f value  for each model term , expressed as a proportion of the critical value  and relate to models lacking interaction terms in the case of the main effects.
c principal co-ordinates  ordination of immune-associated gene expression , indicating strong divergence between summer and winter samples along similar site-specific trajectories; plot showing scatter of individual points against the  <dig> major axes  and concentration ellipsoids containing 50 % of points
we also asked how important seasonal influences on the expression of immune-associated genes were in comparison to other sources of variation .
to answer this we again considered lms fitted to expression data for all  <dig> immune-associated genes, initially with main effects for season, site, sex and body length and then with all 2-way interactions involving season.
as for the analyses of genome-wide expression above, the broad pattern in these models was for season to be the dominant influence on immune gene expression , compared to individual sex, body size or site.
also, the interactions of season with other terms tended to be small compared to the main effect of season, indicating consistent seasonal effects across site, sex and age.
moreover, pco ordination of all immport list genes  revealed clear differentiation between summer and winter samples along similar trajectories between sites .
these observations are consistent with overarching temporal environmental drivers acting similarly on the immune system across different habitat types and life-history stages.
adaptive immunity genes are summer-biased and innate immunity genes are winter-biased
the  <dig> consistently seasonally-biased genes from the immport list were individually evaluated to identify those with core immunological functions  .
such “core” genes in the summer-biased set included those involved in, or regulating, adaptive effector response pathways , corroborating the summer bias in lymphocyte responses suggested by gsea analyses.
one weakly expressed classical major histocompatibility class  iia locus  was also detected more strongly in summer, although this was not the case for other more highly expressed mhciia loci; the chromosome vii locus is hereafter referred to as mhciia.
there were also summer-biased genes involved in immunological cell adhesion  and toll-like receptor -mediated signalling .
in the winter-biased set there was a lack of genes clearly promoting adaptive immunity.
however, there were several genes involved in regulating or suppressing lymphocyte activity  .
in addition there were genes associated with innate immune cell activity  and interleukin   <dig> family signalling pathways , and genes like those up-regulated by type i interferons in mammals  or involved in tlr signalling pathways leading to the production of type i interferons .
key elements of non-classical complement pathways  were also winter-biased.
although selected on the basis of cuffdiff outputs  <cit> , all of the core immune genes were highly significantly seasonally-biased when analysed in confounder-adjusted lms with terms for season, site, sex and body length.
a set of highly co-expressed winter- and summer-biased immune genes can be identified that may lie close to a regulatory axis for seasonal immunity
as the regulatory influence of genes on other genes may be reflected in patterns of co-expression, we analyzed these patterns in our transcriptomic data to identify candidate seasonally-biased genes with high regulatory importance.
given the prominence of the seasonal signature characterized above, we used unadjusted data for these analyses in order to preserve the overall context of co-expression.
we selected  <dig> of the seasonally-biased genes with well-defined immune system functions  as hubs for co-expression networks constructed using an information theory based algorithm .
these hubs included most of the core genes, but excluded those with very strongly linked functions .
in an initial analysis of the hubs alone , summer and winter-biased genes overwhelmingly segregated to different regions of the network, interfacing principally through a relatively small number of genes  that were also highly interconnected within their respective summer or winter-biased sets.
this set of highly connected interfacing genes may lie close to the regulatory axis that controls seasonal changes in immune function and comprised orai <dig>  cd <dig>  tbk <dig> and il1r-like in the winter-biased set and colec <dig>  mhciia, foxp3b, cd8a, zap <dig>  and rag <dig> in the summer-biased set.
remarkably, some of these summer-biased genes  are involved in the t-cell – apc immunological synapse determining t-cell activation.
furthermore, the winter-biased orai <dig> codes for a calcium channel that is necessary for t-cell proliferation  <cit>  and that is recruited to the t-cell receptor complex during activation  <cit> .
loss-of-function mutations affecting mammalian orthologues of six of the genes cause either severe immunodeficiencies  or autoimmunity , reflecting their potential degree of influence upon the adaptive immune response.fig.
<dig>
aracne networks of seasonally-biased core immune genes.
winter- and summer-biased nodes segregate to different regions of networks and interface via a small set of central nodes that are highly connected amongst themselves and also within their respective winter- or summer-biased set.
a network  <dig>  nodes sized according to their betweenness-centrality ; network shown with a force-directed layout, modified to highlight edges  between summer and winter-biased genes .
b network  <dig> re-analyzed with summer-biased timeless added as an extra node; nodes sized according to eccentricity  and shown with an unmodified force-directed layout.
edges from timeless connect to winter-summer interface genes and further winter-biased genes
a different aracne network  including all summer- and winter-biased genes  was also constructed, using the same hub genes as above.
when gene modules from this analysis were considered in terms of the degree of gene overlap , modules associated with winter-biased hubs were largely segregated from those associated with summer-biased hubs.
the strongest similarity across winter- and summer-biased modules was between those associated with the hubs orai <dig> and cd8a, tending to confirm that genes interacting with orai <dig> and cd8a might be close to the axis of seasonal regulation.fig.
<dig>
aracne network  of all seasonally-biased genes, specifying core immune genes and timeless as hubs.
nodes, shown in a modified force-directed layout, represent gene module sizes associated with hubs and edges are jaccard similarity coefficients for module composition .
modules associated with winter and summer-biased immune hubs segregate to different regions of the network, with the strongest winter-summer module similarity between orai <dig> and cd8a.
the module associated with the summer-biased timeless is primarily similar to modules associated with winter-biased hubs.
for larger modules the heat map  shows significant odds ratio gene overlaps with gene sets representing organismal signatures of metabolism, growth, reproduction and stress
a third aracne network  was constructed from all immune associated genes  and, again, the same set of seasonally-biased immune hub genes used above.
in this network, there was a dominant cluster of large, overlapping modules .
most of these were enriched in both summer- and winter-biased genes and associated with hubs that were winter-summer interfacing nodes in network  <dig>  this dominant meta-module  is consistent with the existence of a coherent regulatory unit involved in seasonal immune function.fig.
<dig>
aracne network  of all immune-associated genes , specifying seasonally-biased core immune genes as hubs.
a nodes, shown in a modified force-directed network, represent gene module sizes associated with hubs and edges are jaccard similarity coefficients for module composition .
node colours indicate modules significantly enriched in winter-biased genes, summer-biased genes or both .
modules associated with hubs that were winter-summer interface nodes in network  <dig> tended to be large, to share a high degree of similarity in composition, and to contain significant enrichments of both winter- and summer-biased genes.
b form of simple structural equation model  used in assessing the influence of winter-summer interface  genes from network 1; w, main axis of covariation in winter-biased core immune genes ; s, main axis of covariation in summer-biased core immune genes; grey circle, expression of individual winter-summer interfacing gene
finally, we constructed a small  three-variable structural equation model  of the form shown in fig.
5b.
we used this to further assess the influence of individual winter-summer interfacing  genes from aracne network  <dig> on the seasonal transition in immune function.
in this analysis two of the variables were derived as the first components from separate principal components analyses  of summer and winter-biased core immune genes .
each component thus represented the major axis of covariation within the respective summer or winter-biased gene set.
the third variable was the expression of a key winter-summer interfacing gene, each of which, in turn, was evaluated in the model.
all of the winter-summer interfacing immune genes, except tbk <dig>  negated the direct effect of winter-biased on summer biased genes  and themselves showed significant associations, of opposite sign, with the summer and winter-biased genes.
this supports the linking role of these genes indicated by the aracne analyses.
consistent seasonality confirmed by year-round q-pcr measurement of key genes over a new annual cycle
arbitrarily selecting  <dig>  from the  <dig> key genes identified in network  <dig> above, we confirmed their seasonal expression bias by quantitative real-time pcr  in an independent sample set.
whole-fish gene expression measurements were made for  <dig> individuals sampled in a regular monthly design over a new annual cycle  in our original localities  and also in semi-natural artificial habitats populated by stock from frn.
all of the genes showed highly significant confounder-adjusted seasonal patterns , overall, with a clear peak in the predicted season.fig.
<dig> corroborative whole-fish and tissue-specific q-pcr gene expression measurements.
a temporal variation from october  <dig> to september  <dig> in whole-fish gene expression for winter-summer interface  genes from network  <dig>  n =  <dig>  relative expression   is indicated on the y-axis.
plots show thin-plate-spline smoothers for time fitted in a generalized additive model  with fixed effects for habitat, sex and length; shaded areas represent 95 % confidence regions.
samples were derived from frn, rhd and artificial outdoors habitats stocked from frn.
the other key gene examined, il1r-like, also demonstrated significant seasonal variation  with peak expression in december , if log <dig> transformed.
b tissue-specific expression of key genes at sto .
heat map showing relative gene expression across tissues; significant differences occurred for all genes .
c tissue-specific seasonal variation in key genes at sto.
mean relative expression  ±  <dig> se is shown on the y-axis.
p values  relate to directional  t-tests of seasonal shifts in the same direction as the whole-fish rnaseq study; 13/ <dig> of these tests were significant but, in comparison , 0/ <dig> of 1-tailed tests in the opposite direction were significant.
the calibrator sample for tissue-specific analyses was pooled whole-fish rna from  <dig> individuals in september.
the sto samples showed no significant difference in fish length between winter and summer and were balanced for sex ratio
tissue-specific expression of key genes suggests intense seasonality in the gill
we also considered seasonality of key genes  within specific tissues  at a new locality on the river stour in eastern england .
all of these genes, except for orai <dig>  were primarily expressed in organs with known concentrations of lymphoid tissue .
furthermore, all of the many instances of significant tissue-specific seasonal bias  occurred in the same direction as predicted by the whole-fish transcriptomic study .
outside of the thymus, expression of t-cell-associated genes  was highest in the gill, lower in head kidney, spleen and intestine and negligible in skeletal muscle .
this is consistent with a strong concentration of t-cell activity in the gill.
moreover, the summer bias of t-cell-associated genes was seen primarily in the gill .
in the case of orai <dig>  whose expression is important in mammalian t-cells  <cit>  but not narrowly characteristic of them, high expression occurred in skeletal muscle  .
this gene was, however, also robustly expressed  and winter-biased  in the organs with greatest t-cell-specific expression, thymus and gill, supporting a possible role in seasonal immunoregulation.
genes from innate signalling pathways  tended to be winter-biased in all tissues .
overall, gill most closely reflected the pattern of seasonal bias seen in whole-fish mrna pools .
seasonal immune gene expression links to wider life history signatures
some seasonal immune functions in vertebrates are controlled by photoperiodic time measurement  <cit>  and the circadian molecular clock may also have a role in co-ordinating circannual biological rhythms  <cit> .
a scan of the seasonally-biased genes for those involved in such processes revealed that timeless  occurred within the summer-biased set.
when added to the aracne networks above timeless was, remarkably, most strongly connected to key winter-summer interface genes, with more connections to winter-biased genes .
however, other genes involved in clock machinery or photoperiodism did not show the same tendency and timeless is known to have physiological functions in mammals that are independent of any role in biological clocks  <cit> .
in order to place the seasonal variation in immune pathways within an even wider organismal context, we next tested the larger gene modules from the aracne network of all seasonally-biased genes  for overlap with the gene sets representing other life history components and responses to stressors .
the modules formed by summer-biased hubs cd8a and zap <dig> contained genes associated with organismal growth.
in contrast, modules for winter-biased hubs contained genes associated with metabolism and responses to stress .
to explore these associations further, we then constructed an additional network  including seasonally-biased core immune genes, timeless, and representatives from the wider organismal gene sets.
the latter were selected on the basis that they were relatively differentially seasonally expressed in gsea analyses.
in the resulting network , in which all genes were specified as hubs, those representing growth processes were again associated primarily with summer-biased adaptive immune genes .
apart from some connections to growth-related genes , mostly the edges emanating from winter-biased genes were to timeless and to genes representing metabolism and oxidative stress.
notably, winter-summer interfacing  genes from network  <dig> showed especially high numbers of non-immune edges  in network  <dig>  .
as connectivity may reflect regulatory influence, this reinforces the potential regulatory importance of key winter-summer interface genes from network  <dig>  in fact, only two non-interface genes showed high numbers of non-immune edges in the network.
these were the innate winter-biased genes il1rap, which showed a pattern of edges similar to il1r-like, and, more distinctively, tnfrsf <dig>  which was the main immune node with connections to general stress responses.fig.
<dig> association between seasonal core immune genes and wider organismal signatures of growth, metabolism and stress.
a
aracne network  including seasonal core immune genes, timeless, and seasonal genes from curated sets representing growth, metabolism, and aspects of stress ; specifying all genes as hubs; nodes sized according to their betweenness-centrality.
network shown in an unmodified force-directed layout.
b bar chart showing, for the most highly connected core immune genes in network  <dig>  the distribution of edges with genes representing wider organismal signatures .
colour bar outside vertical axis indicates winter- or summer expression bias.
c scatterplot showing, for core immune genes in network  <dig>  numbers of edges to other immune genes vs numbers of edges to non-immune genes ; point sizes are proportional to the number of winter-summer edges for the gene in network  <dig>
tlr signalling pathways show seasonal modulation
finally, given the importance of tbk <dig> in gene co-expression networks  and also the contrasting seasonal expression of tirap  we examined tlr signalling pathways in more detail .
this used genes from the kegg tlr-signalling pathway supplemented by a few further a priori selected genes associated with tlr function in fishes .
there was a clear seasonal modulation.
most seasonally-biased genes  tended to be winter-biased, although some  were strongly summer-biased.
there was considerable discrimination between summer and winter samples when tlr signalling gene expression was ordinated using pco .fig.
<dig> seasonal bias in toll-like receptor  signalling pathway.
a differential winter-summer gene expression in pathway members , demonstrating winter bias in some cases and summer bias in others.
categorization of differential expression is based on overall significance levels in general linear models  with explanatory terms for site, sex and length .
b principal co-ordinates  ordination of gene expression in all pathway members, revealing considerable winter-summer discrimination; scatter of individual points against the  <dig> major axes  and concentration ellipsoids containing 50 % of points
discussion
we have demonstrated seasonal re-adjustments of immune system gene expression in naturally-occurring freshwater teleosts.
these occurred most intensely in the gill and were substantial  and over-arching .
in keeping with some previous suggestions about seasonal immune function in teleosts  <cit> , we found that genes marking adaptive immune processes were summer-biased , whilst certain innate immune genes were winter-biased.
however, as set out below, our observations provide considerable new insights into the control of seasonal immune responses.
transcriptomic analyses  indicated that summer-biased genes included many centrally involved in lymphocyte responses.
for example, the recombination activating genes  and genes associated with particular adaptive cell populations: t-cells , cytotoxic t-cells , helper t-cells  and b-cells .
in contrast, the set of winter-biased genes lacked those promoting adaptive effector responses.
in all cases, winter-biased genes associated with t- or b-cell responses were regulatory or even suppressive in nature.
this strongly suggests a regulatory control of adaptive immunity during winter, rather than, or additional to, a loss of function due to the kinetic consequences of low temperature in a cold-blooded organism.
furthermore, there were gene expression signatures of elevated innate immune functions in winter: including il- <dig> signalling and non-classical complement pathways.
a complex modulation of genes involved in innate tlr-mediated signalling occurred, with a predominant winter bias.
we designed the sampling for our transcriptomic study to, as far as possible, reduce correlation between season and ontogeny, and we also carefully considered, post-hoc, the possible role of ontogeny in generating apparent seasonal differences.
to ensure that an extensively overlapping range of effective ages was present in our winter and summer transcriptomic samples, we deliberately selected a wide range of fish sizes within samples .
through monitoring the growth of fishes in artificial outdoor habitats we confirmed that age predicted the majority of variation in length, and we adjusted for length, as a surrogate for age, in statistical models applied to transcriptomic data.
importantly, the much greater overall signature of season compared to length in statistical models applied to genome-wide and immune system-wide gene expression is not consistent with ontogeny being a major confounder in our study.
moreover, the balancing of size across seasonal samples, and the adjustment for length in our statistical modelling, also accounts for the possibility that growth allometries in different tissues  may have biased results for whole-fish samples.
we also considered whether the patterns of gene co-expression in our transcriptomic data could give insights into the regulation of seasonal immune function.
information theory-based network analyses  <cit>  of expression in seasonally-biased core immune system genes identified a small set of genes lying at the interface between summer- and winter-biased genes.
these were highly networked  amongst themselves and also each highly networked within the seasonally-biased group to which they respectively belonged.
remarkably, several of the interfacing genes have roles in apc-t-cell immunological synapses  and mutations leading to loss of function in their mammalian orthologues cause primary immunodeficiencies .
also amongst these interface genes, foxp <dig>  <cit>  is a master regulator of regulatory t-cell function and recombination-activating genes are central to the production of re-combined adaptive receptors  <cit> .
in mammals loss-of-function mutations in these genes respectively cause lethal autoimmunity and severe combined immunodeficiency  <cit> .
other interfacing genes are involved in innate processes that might precede antigen presentation: innate signalling pathways  and antigen internalization via phagocytosis or endocytosis .
when we added the summer-biased clock-associated gene timeless to these networks, it proved to be closely associated with interface genes, and especially with winter-biased  interface genes.
whilst this could reflect some co-ordination via a seasonal oscillator, though, other genes involved in photoperiodism or circadian rhythms did not enter the networks in corresponding ways.
it is also the case that timeless itself has a relatively poorly resolved role in the mammalian circadian clock and is known to have independent physiological functions  <cit>  .
during winter there was a genome-wide signature indicative of elevated metabolic processes and metabolite transfer and organismal stress, and in summer a signature of growth and developmental processes.
again using network analyses of our transcriptomic data, we finally asked how seasonal changes in immunity might be related to this background.
we found that the winter-summer interfacing  immune genes identified above were especially highly connected to genes involved in non-immune seasonal variation, further emphasizing their relevance in the seasonal control of immunity.
genes involved in metabolism and oxidative stress interconnected densely with winter-biased innate genes, and amongst these especially to the winter-summer interfacing genes tbk <dig> and il1r-like, and also to il1rap.
on the other hand, genes involved in general organismal stress responses linked differently to winter-biased immune responses: primarily via tnfrsf <dig>  a protein that triggers apoptotic pathways  <cit>  and restrains t-cell  <cit>  and b-cell  <cit>  responses.
in comparison, genes involved in non-immune summer signatures  networked primarily to summer-biased adaptive genes, especially to the summer-biased interface genes cd8a and zap <dig>  these observations suggest an unexpectedly strong link between growth processes and adaptive immunity and that one, or both, may favour permissive conditions for the other.
taken together, the above patterns indicate that multiple organismal processes are likely to interact with the seasonal regulation of immunity, additional to the possible influence of any “hard-wired” circannual oscillator.
it might be expected, then, that predictable seasonal influences will be modified by less predictable non-cyclical temporal variations in environmental stressors  <cit> .
to validate our transcriptomic analyses we returned to our original study localities, an upland lake and river in mid-wales, and also considered artificial outdoors habitats stocked from the lake site.
using q-pcr measurements we confirmed  seasonality in a panel of the key immunity genes predicted to be winter- or summer-biased.
this year-round monthly analysis considered whole-fish samples  and was carried out in a new annual cycle that lacked unusually cold winter or spring weather.
furthermore tissue-specific analyses  at an entirely new locality  found that all tissues showed seasonal expression changes and these changes all occurred in the same direction as in the whole-fish studies at our original sites.
thus, overall we considered  <dig> very divergent localities  across 2 years and found compelling evidence to support a general pattern such as that indicated in our initial transcriptomic measurements.
in addition to our analyses of whole-fish mrna pools, we also confirmed tissue-specific expression patterns through q-pcr measurements of key immunity genes at a new estuarine locality .
as indicated above we found many significant tissue-specific seasonal expression differences and all of these were in the direction predicted by our other whole-fish studies.
the most pronounced seasonal expression profile occurred in the gill .
furthermore, the gill contained the most intense concentration of t-cell activity outside of the thymus, with elevated expression of t-cell associated genes such as cd8a and foxp3b, and expression of these genes was seasonal in the gill but non-seasonal in the thymus.
these observations are consistent with the known responsiveness of immune gene expression in the teleost gill to environmental stimuli  <cit> , and also with the recent discovery and characterization of extensive, t-cell rich, interbranchial lymphoid tissue  in teleost fishes .
our results suggest the possibility that ilt may have an important role in seasonal immune function.
finally, and whilst the present study is intended to characterize the seasonal dynamics of gene expression, rather than identify environmental causation, we briefly consider what external agents may drive the responses that we observed.
in highly seasonal temperate zone habitats, such as the ones we consider here, each of temperature, diet, photoperiodic responses, pathogen exposures, or other biotic or abiotic manifestations of the environment, could be involved to unknown degrees.
in the future, by matching detailed field observations with mesocosm studies and laboratory experiments, we expect to dissect the relative contributions of these influences to seasonal immune variation and to the immune phenotype more generally.
CONCLUSIONS
our results suggest that in wild teleosts, during winter conditions, adaptive immune activity declines in a manner that involves the expression of regulatory genes affecting lymphocyte function.
this is indicative of a controlled, strategic response rather than a simple kinetic tracking of environmental temperature.
seasonal change is most prominent in the gill, suggesting ilt may be important in such responses.
further broad attention to seasonal immune function is certainly warranted, given the likely practical relevance – through effects on infectious disease susceptibility and inflammatory status – to health in humans and domesticated animals and to fitness in natural populations.

 seasonal immune function has often been observed in vertebrates, including humans, but is relatively poorly understood. as with more studied circadian rhythms, though, there are fundamental implications for health. taking a comparative approach and considering the conserved genes of the vertebrate immune system, here we use transcriptomic measurements to reveal seasonal re-adjustments of immunity in a naturally-occurring teleost model. crucially, through focussing on wild organisms exposed to real-world environmental extremes, we expected to discover more measurable variation than would be the case in a study of domesticated animals, where seasonal variation may be muted by anthropogenic influences. we chose the three-spined stickleback as a subject because it has an annotated whole genome sequence and occurs accessibly in highly seasonal natural habitats. also, a considerable knowledge base exists for this species: it is a highly studied model organism, and there are particularly detailed ecological studies relating to our main study area, mid wales. we compared the transcriptomes of populations in late winter and late summer in ecologically divergent natural populations, reasoning that a focus on common responses would provide a way to filter overriding seasonal trends from locality-specific variation. we chose to primarily use global mrna extracts from individual whole fishes rather than from isolated cell populations or tissues. this was because a fully reductionist approach to cell populations would be impractical, and because the majority of the teleost immune system is likely to be diffusely distributed in the gut, under the skin and mucosal surfaces and in association with the gills and liver. by considering global samples, we were thus able to take a holistic view of which immune system pathways are differentially expressed at seasonal extremes. we corroborated our transcriptomic analyses by targeted gene expression measurements in year-round samples of fishes from the original sites in a new annual cycle, and by tissue-specific analyses at a further site. moreover, using network analyses we were able to ask what genes are important in regulating seasonal immune function and how do seasonally-biased immune networks interact with other seasonally-biased organismal processes?  seasonal expression bias of immune system genes occurs against a well-defined genome-wide seasonal signature we analyzed the global transcriptomes of g. aculeatus from two contrasting habitats in mid wales, river rheidol and lake frongoch, in september and march to begin our analysis we considered, genome-wide, which genes were associated with seasonal expression bias. at frn, genes were significantly differentially expressed from summer to winter with an individual cut-off and with a false discovery rate -adjusted cut-off. at rhd, genes were significantly differentially expressed with an individual cut-off and with an fdr-adjusted cut-off. genes that were seasonally differentially expressed at both localities tended, overwhelmingly, to show synchronous expression. we hypothesized that these synchronously differentially expressed genes would also be those contributing to overarching seasonal responses. thus, we categorized such genes on the basis that they were significantly differentially expressed, in the same direction, at both sites at the individual error rate. following this criterion, genes were differentially expressed in a consistent direction, increasing expression during winter and increasing expression during summer. we then considered genome-wide winter-summer expression changes in functional gene sets using gene set enrichment analyses, focussing on predicted orthologues of human genes. we compared ranked gene expression to kegg and reactome pathways, globally. we also compared the ranked genes to a selected group of gene sets reflecting different organismal processes and immunological pathways. finally, we considered hypergeometric overlap of the selected gene sets with the summer- and winter-biased genes defined above. taken together, these analyses indicated signatures of growth in summer and physiological challenge in winter. in all analyses, some immune pathways were seasonally-biased, with innate processes emphasised in winter and adaptive processes in summer.fig. distinctive immunological and genome-wide gene expression signatures occurred at seasonal extremes. a gene sets with significant summer or winter expression bias as indicated by gene set enrichment analysis. ranked differential gene expression was compared, separately for the rhd and frn sites, to global kegg and reactome gene sets, and sets are shown where the combined fdr p value was significant; gene set names are truncated but shown in full in additional file 2: table s2; stars indicate immunologically-relevant gene sets; the central dendrogram indicates the degree of overlap between gene sets. b analyses of selected gene sets representing immunological pathways and organismal signatures of stress, reproduction, growth and metabolism. individual colour panels correspond, left to right, to the order of gene sets in additional file 3: table s these sets were considered by gsea for rhd and frn separately, and by gene overlap for the overall summer and winter-biased gene sets the above analyses were carried out on expression data un-adjusted for individual size, as this variable was approximately balanced across winter and summer samples. however, as our sampling points bounded a non-recruiting population ageing in the interval between breeding seasons, we considered in more detail the potential influence of ontogenetic stage. it is likely, given the months in which we recorded reproductive activity in the field, and taking into account slower biological ageing at lower temperatures, that the 0+ cohort in our approximately size-matched summer and winter samples would have included individuals widely overlapping in effective age. furthermore, as our sampling deliberately selected a wide range of fish sizes, it is probable that 0+ and 1+ cohorts were represented, resulting in a very extensive overlap of effective ages between summer and winter samples. a close association between body size and age allows age to be partitioned from season in statistical models by the use of a size metric, such as body length, as a surrogate. this is validated by data from experiments in artificial outdoor habitats, where we found that time explains at least 57 % of the variation in length. in order to control for age effects we applied general linear models to each in turn of the genes in the gsea dataset, including main effects for season, length, sex and site. we found that season was the dominant predictor of gene expression. consistent with this gene-by-gene analysis, a multivariate principal co-ordinates analysis of the same data demonstrated marked differentiation across seasons against the two major axes, but none for length, sex or site. we also re-ran “global” gsea analyses, first with genes ranked by confounder-adjusted seasonal effect, and then with genes ranked by confounder-adjusted length effects. we found a similar outcome in the analysis ranked by confounder-adjusted seasonal effect to in the unadjusted analysis shown in fig. 1a, with the two analyses sharing gene sets that were significantly seasonally enriched, including all of the immunological sets except leucocyte transendothelial migration. in contrast, there was a distinctive outcome in the analysis ranked by confounder-adjusted length effect, where only enriched gene sets were shared with the analysis shown in fig. 1a, including none of the immunological sets. thus, the effect of season was a very dominant one, emerging clearly in analyses even without adjustment for ontogeny. season is a dominant and consistent influence on immune gene expression focussing on the immune system we considered genes with predicted orthology to those in the immport comprehensive list of immune-associated genes. in total, immune-associated genes out of were seasonally-biased. of these, greater absolute numbers of genes were winter-biased than summer-biased. as with the full gene list, immune-associated genes that were consistently seasonally differentially-expressed tended, very strongly, to be synchronously expressed across localities.fig. extremes of season were the dominant predictor of immune gene expression. a scatterplot of log winter-summer fold expression change for all immune-associated genes with significant seasonal difference at both rhd and frn sites. overwhelmingly such genes were regulated in the same direction across sites. b season was overwhelmingly the most important predictor of immune gene expression, in comparison to site, sex and body size. bar chart summarizes results from general linear models fitted to each individual log immune-associated gene expression variable; bars are the mean observed f value for each model term, expressed as a proportion of the critical value and relate to models lacking interaction terms in the case of the main effects. c principal co-ordinates ordination of immune-associated gene expression, indicating strong divergence between summer and winter samples along similar site-specific trajectories; plot showing scatter of individual points against the major axes and concentration ellipsoids containing 50 % of points we also asked how important seasonal influences on the expression of immune-associated genes were in comparison to other sources of variation. to answer this we again considered lms fitted to expression data for all immune-associated genes, initially with main effects for season, site, sex and body length and then with all 2-way interactions involving season. as for the analyses of genome-wide expression above, the broad pattern in these models was for season to be the dominant influence on immune gene expression, compared to individual sex, body size or site. also, the interactions of season with other terms tended to be small compared to the main effect of season, indicating consistent seasonal effects across site, sex and age. moreover, pco ordination of all immport list genes revealed clear differentiation between summer and winter samples along similar trajectories between sites. these observations are consistent with overarching temporal environmental drivers acting similarly on the immune system across different habitat types and life-history stages. adaptive immunity genes are summer-biased and innate immunity genes are winter-biased the consistently seasonally-biased genes from the immport list were individually evaluated to identify those with core immunological functions. such “core” genes in the summer-biased set included those involved in, or regulating, adaptive effector response pathways, corroborating the summer bias in lymphocyte responses suggested by gsea analyses. one weakly expressed classical major histocompatibility class iia locus was also detected more strongly in summer, although this was not the case for other more highly expressed mhciia loci; the chromosome vii locus is hereafter referred to as mhciia. there were also summer-biased genes involved in immunological cell adhesion and toll-like receptor -mediated signalling. in the winter-biased set there was a lack of genes clearly promoting adaptive immunity. however, there were several genes involved in regulating or suppressing lymphocyte activity. in addition there were genes associated with innate immune cell activity and interleukin family signalling pathways, and genes like those up-regulated by type i interferons in mammals or involved in tlr signalling pathways leading to the production of type i interferons. key elements of non-classical complement pathways were also winter-biased. although selected on the basis of cuffdiff outputs, all of the core immune genes were highly significantly seasonally-biased when analysed in confounder-adjusted lms with terms for season, site, sex and body length. a set of highly co-expressed winter- and summer-biased immune genes can be identified that may lie close to a regulatory axis for seasonal immunity as the regulatory influence of genes on other genes may be reflected in patterns of co-expression, we analyzed these patterns in our transcriptomic data to identify candidate seasonally-biased genes with high regulatory importance. given the prominence of the seasonal signature characterized above, we used unadjusted data for these analyses in order to preserve the overall context of co-expression. we selected of the seasonally-biased genes with well-defined immune system functions as hubs for co-expression networks constructed using an information theory based algorithm. these hubs included most of the core genes, but excluded those with very strongly linked functions. in an initial analysis of the hubs alone, summer and winter-biased genes overwhelmingly segregated to different regions of the network, interfacing principally through a relatively small number of genes that were also highly interconnected within their respective summer or winter-biased sets. this set of highly connected interfacing genes may lie close to the regulatory axis that controls seasonal changes in immune function and comprised orai cd tbk and il1r-like in the winter-biased set and colec mhciia, foxp3b, cd8a, zap and rag in the summer-biased set. remarkably, some of these summer-biased genes are involved in the t-cell – apc immunological synapse determining t-cell activation. furthermore, the winter-biased orai codes for a calcium channel that is necessary for t-cell proliferation and that is recruited to the t-cell receptor complex during activation. loss-of-function mutations affecting mammalian orthologues of six of the genes cause either severe immunodeficiencies or autoimmunity, reflecting their potential degree of influence upon the adaptive immune response.fig. aracne networks of seasonally-biased core immune genes. winter- and summer-biased nodes segregate to different regions of networks and interface via a small set of central nodes that are highly connected amongst themselves and also within their respective winter- or summer-biased set. a network nodes sized according to their betweenness-centrality; network shown with a force-directed layout, modified to highlight edges between summer and winter-biased genes. b network re-analyzed with summer-biased timeless added as an extra node; nodes sized according to eccentricity and shown with an unmodified force-directed layout. edges from timeless connect to winter-summer interface genes and further winter-biased genes a different aracne network including all summer- and winter-biased genes was also constructed, using the same hub genes as above. when gene modules from this analysis were considered in terms of the degree of gene overlap, modules associated with winter-biased hubs were largely segregated from those associated with summer-biased hubs. the strongest similarity across winter- and summer-biased modules was between those associated with the hubs orai and cd8a, tending to confirm that genes interacting with orai and cd8a might be close to the axis of seasonal regulation.fig. aracne network of all seasonally-biased genes, specifying core immune genes and timeless as hubs. nodes, shown in a modified force-directed layout, represent gene module sizes associated with hubs and edges are jaccard similarity coefficients for module composition. modules associated with winter and summer-biased immune hubs segregate to different regions of the network, with the strongest winter-summer module similarity between orai and cd8a. the module associated with the summer-biased timeless is primarily similar to modules associated with winter-biased hubs. for larger modules the heat map shows significant odds ratio gene overlaps with gene sets representing organismal signatures of metabolism, growth, reproduction and stress a third aracne network was constructed from all immune associated genes and, again, the same set of seasonally-biased immune hub genes used above. in this network, there was a dominant cluster of large, overlapping modules. most of these were enriched in both summer- and winter-biased genes and associated with hubs that were winter-summer interfacing nodes in network this dominant meta-module is consistent with the existence of a coherent regulatory unit involved in seasonal immune function.fig. aracne network of all immune-associated genes, specifying seasonally-biased core immune genes as hubs. a nodes, shown in a modified force-directed network, represent gene module sizes associated with hubs and edges are jaccard similarity coefficients for module composition. node colours indicate modules significantly enriched in winter-biased genes, summer-biased genes or both. modules associated with hubs that were winter-summer interface nodes in network tended to be large, to share a high degree of similarity in composition, and to contain significant enrichments of both winter- and summer-biased genes. b form of simple structural equation model used in assessing the influence of winter-summer interface genes from network 1; w, main axis of covariation in winter-biased core immune genes; s, main axis of covariation in summer-biased core immune genes; grey circle, expression of individual winter-summer interfacing gene finally, we constructed a small three-variable structural equation model of the form shown in fig. 5b. we used this to further assess the influence of individual winter-summer interfacing genes from aracne network on the seasonal transition in immune function. in this analysis two of the variables were derived as the first components from separate principal components analyses of summer and winter-biased core immune genes. each component thus represented the major axis of covariation within the respective summer or winter-biased gene set. the third variable was the expression of a key winter-summer interfacing gene, each of which, in turn, was evaluated in the model. all of the winter-summer interfacing immune genes, except tbk negated the direct effect of winter-biased on summer biased genes and themselves showed significant associations, of opposite sign, with the summer and winter-biased genes. this supports the linking role of these genes indicated by the aracne analyses. consistent seasonality confirmed by year-round q-pcr measurement of key genes over a new annual cycle arbitrarily selecting from the key genes identified in network above, we confirmed their seasonal expression bias by quantitative real-time pcr in an independent sample set. whole-fish gene expression measurements were made for individuals sampled in a regular monthly design over a new annual cycle in our original localities and also in semi-natural artificial habitats populated by stock from frn. all of the genes showed highly significant confounder-adjusted seasonal patterns, overall, with a clear peak in the predicted season.fig. corroborative whole-fish and tissue-specific q-pcr gene expression measurements. a temporal variation from october to september in whole-fish gene expression for winter-summer interface genes from network n =  relative expression is indicated on the y-axis. plots show thin-plate-spline smoothers for time fitted in a generalized additive model with fixed effects for habitat, sex and length; shaded areas represent 95 % confidence regions. samples were derived from frn, rhd and artificial outdoors habitats stocked from frn. the other key gene examined, il1r-like, also demonstrated significant seasonal variation with peak expression in december, if log transformed. b tissue-specific expression of key genes at sto. heat map showing relative gene expression across tissues; significant differences occurred for all genes. c tissue-specific seasonal variation in key genes at sto. mean relative expression  ±  se is shown on the y-axis. p values relate to directional t-tests of seasonal shifts in the same direction as the whole-fish rnaseq study; 13/ of these tests were significant but, in comparison, 0/ of 1-tailed tests in the opposite direction were significant. the calibrator sample for tissue-specific analyses was pooled whole-fish rna from individuals in september. the sto samples showed no significant difference in fish length between winter and summer and were balanced for sex ratio tissue-specific expression of key genes suggests intense seasonality in the gill we also considered seasonality of key genes within specific tissues at a new locality on the river stour in eastern england. all of these genes, except for orai were primarily expressed in organs with known concentrations of lymphoid tissue. furthermore, all of the many instances of significant tissue-specific seasonal bias occurred in the same direction as predicted by the whole-fish transcriptomic study. outside of the thymus, expression of t-cell-associated genes was highest in the gill, lower in head kidney, spleen and intestine and negligible in skeletal muscle. this is consistent with a strong concentration of t-cell activity in the gill. moreover, the summer bias of t-cell-associated genes was seen primarily in the gill. in the case of orai whose expression is important in mammalian t-cells but not narrowly characteristic of them, high expression occurred in skeletal muscle. this gene was, however, also robustly expressed and winter-biased in the organs with greatest t-cell-specific expression, thymus and gill, supporting a possible role in seasonal immunoregulation. genes from innate signalling pathways tended to be winter-biased in all tissues. overall, gill most closely reflected the pattern of seasonal bias seen in whole-fish mrna pools. seasonal immune gene expression links to wider life history signatures some seasonal immune functions in vertebrates are controlled by photoperiodic time measurement and the circadian molecular clock may also have a role in co-ordinating circannual biological rhythms. a scan of the seasonally-biased genes for those involved in such processes revealed that timeless occurred within the summer-biased set. when added to the aracne networks above timeless was, remarkably, most strongly connected to key winter-summer interface genes, with more connections to winter-biased genes. however, other genes involved in clock machinery or photoperiodism did not show the same tendency and timeless is known to have physiological functions in mammals that are independent of any role in biological clocks. in order to place the seasonal variation in immune pathways within an even wider organismal context, we next tested the larger gene modules from the aracne network of all seasonally-biased genes for overlap with the gene sets representing other life history components and responses to stressors. the modules formed by summer-biased hubs cd8a and zap contained genes associated with organismal growth. in contrast, modules for winter-biased hubs contained genes associated with metabolism and responses to stress. to explore these associations further, we then constructed an additional network including seasonally-biased core immune genes, timeless, and representatives from the wider organismal gene sets. the latter were selected on the basis that they were relatively differentially seasonally expressed in gsea analyses. in the resulting network, in which all genes were specified as hubs, those representing growth processes were again associated primarily with summer-biased adaptive immune genes. apart from some connections to growth-related genes, mostly the edges emanating from winter-biased genes were to timeless and to genes representing metabolism and oxidative stress. notably, winter-summer interfacing genes from network showed especially high numbers of non-immune edges in network. as connectivity may reflect regulatory influence, this reinforces the potential regulatory importance of key winter-summer interface genes from network in fact, only two non-interface genes showed high numbers of non-immune edges in the network. these were the innate winter-biased genes il1rap, which showed a pattern of edges similar to il1r-like, and, more distinctively, tnfrsf which was the main immune node with connections to general stress responses.fig. association between seasonal core immune genes and wider organismal signatures of growth, metabolism and stress. a aracne network including seasonal core immune genes, timeless, and seasonal genes from curated sets representing growth, metabolism, and aspects of stress; specifying all genes as hubs; nodes sized according to their betweenness-centrality. network shown in an unmodified force-directed layout. b bar chart showing, for the most highly connected core immune genes in network the distribution of edges with genes representing wider organismal signatures. colour bar outside vertical axis indicates winter- or summer expression bias. c scatterplot showing, for core immune genes in network numbers of edges to other immune genes vs numbers of edges to non-immune genes; point sizes are proportional to the number of winter-summer edges for the gene in network tlr signalling pathways show seasonal modulation finally, given the importance of tbk in gene co-expression networks and also the contrasting seasonal expression of tirap we examined tlr signalling pathways in more detail. this used genes from the kegg tlr-signalling pathway supplemented by a few further a priori selected genes associated with tlr function in fishes. there was a clear seasonal modulation. most seasonally-biased genes tended to be winter-biased, although some were strongly summer-biased. there was considerable discrimination between summer and winter samples when tlr signalling gene expression was ordinated using pco.fig. seasonal bias in toll-like receptor signalling pathway. a differential winter-summer gene expression in pathway members, demonstrating winter bias in some cases and summer bias in others. categorization of differential expression is based on overall significance levels in general linear models with explanatory terms for site, sex and length. b principal co-ordinates ordination of gene expression in all pathway members, revealing considerable winter-summer discrimination; scatter of individual points against the major axes and concentration ellipsoids containing 50 % of points discussion we have demonstrated seasonal re-adjustments of immune system gene expression in naturally-occurring freshwater teleosts. these occurred most intensely in the gill and were substantial and over-arching. in keeping with some previous suggestions about seasonal immune function in teleosts, we found that genes marking adaptive immune processes were summer-biased, whilst certain innate immune genes were winter-biased. however, as set out below, our observations provide considerable new insights into the control of seasonal immune responses. transcriptomic analyses indicated that summer-biased genes included many centrally involved in lymphocyte responses. for example, the recombination activating genes and genes associated with particular adaptive cell populations: t-cells, cytotoxic t-cells, helper t-cells and b-cells. in contrast, the set of winter-biased genes lacked those promoting adaptive effector responses. in all cases, winter-biased genes associated with t- or b-cell responses were regulatory or even suppressive in nature. this strongly suggests a regulatory control of adaptive immunity during winter, rather than, or additional to, a loss of function due to the kinetic consequences of low temperature in a cold-blooded organism. furthermore, there were gene expression signatures of elevated innate immune functions in winter: including il- signalling and non-classical complement pathways. a complex modulation of genes involved in innate tlr-mediated signalling occurred, with a predominant winter bias. we designed the sampling for our transcriptomic study to, as far as possible, reduce correlation between season and ontogeny, and we also carefully considered, post-hoc, the possible role of ontogeny in generating apparent seasonal differences. to ensure that an extensively overlapping range of effective ages was present in our winter and summer transcriptomic samples, we deliberately selected a wide range of fish sizes within samples. through monitoring the growth of fishes in artificial outdoor habitats we confirmed that age predicted the majority of variation in length, and we adjusted for length, as a surrogate for age, in statistical models applied to transcriptomic data. importantly, the much greater overall signature of season compared to length in statistical models applied to genome-wide and immune system-wide gene expression is not consistent with ontogeny being a major confounder in our study. moreover, the balancing of size across seasonal samples, and the adjustment for length in our statistical modelling, also accounts for the possibility that growth allometries in different tissues may have biased results for whole-fish samples. we also considered whether the patterns of gene co-expression in our transcriptomic data could give insights into the regulation of seasonal immune function. information theory-based network analyses of expression in seasonally-biased core immune system genes identified a small set of genes lying at the interface between summer- and winter-biased genes. these were highly networked amongst themselves and also each highly networked within the seasonally-biased group to which they respectively belonged. remarkably, several of the interfacing genes have roles in apc-t-cell immunological synapses and mutations leading to loss of function in their mammalian orthologues cause primary immunodeficiencies. also amongst these interface genes, foxp is a master regulator of regulatory t-cell function and recombination-activating genes are central to the production of re-combined adaptive receptors. in mammals loss-of-function mutations in these genes respectively cause lethal autoimmunity and severe combined immunodeficiency. other interfacing genes are involved in innate processes that might precede antigen presentation: innate signalling pathways and antigen internalization via phagocytosis or endocytosis. when we added the summer-biased clock-associated gene timeless to these networks, it proved to be closely associated with interface genes, and especially with winter-biased interface genes. whilst this could reflect some co-ordination via a seasonal oscillator, though, other genes involved in photoperiodism or circadian rhythms did not enter the networks in corresponding ways. it is also the case that timeless itself has a relatively poorly resolved role in the mammalian circadian clock and is known to have independent physiological functions. during winter there was a genome-wide signature indicative of elevated metabolic processes and metabolite transfer and organismal stress, and in summer a signature of growth and developmental processes. again using network analyses of our transcriptomic data, we finally asked how seasonal changes in immunity might be related to this background. we found that the winter-summer interfacing immune genes identified above were especially highly connected to genes involved in non-immune seasonal variation, further emphasizing their relevance in the seasonal control of immunity. genes involved in metabolism and oxidative stress interconnected densely with winter-biased innate genes, and amongst these especially to the winter-summer interfacing genes tbk and il1r-like, and also to il1rap. on the other hand, genes involved in general organismal stress responses linked differently to winter-biased immune responses: primarily via tnfrsf a protein that triggers apoptotic pathways and restrains t-cell and b-cell responses. in comparison, genes involved in non-immune summer signatures networked primarily to summer-biased adaptive genes, especially to the summer-biased interface genes cd8a and zap these observations suggest an unexpectedly strong link between growth processes and adaptive immunity and that one, or both, may favour permissive conditions for the other. taken together, the above patterns indicate that multiple organismal processes are likely to interact with the seasonal regulation of immunity, additional to the possible influence of any “hard-wired” circannual oscillator. it might be expected, then, that predictable seasonal influences will be modified by less predictable non-cyclical temporal variations in environmental stressors. to validate our transcriptomic analyses we returned to our original study localities, an upland lake and river in mid-wales, and also considered artificial outdoors habitats stocked from the lake site. using q-pcr measurements we confirmed seasonality in a panel of the key immunity genes predicted to be winter- or summer-biased. this year-round monthly analysis considered whole-fish samples and was carried out in a new annual cycle that lacked unusually cold winter or spring weather. furthermore tissue-specific analyses at an entirely new locality found that all tissues showed seasonal expression changes and these changes all occurred in the same direction as in the whole-fish studies at our original sites. thus, overall we considered very divergent localities across 2 years and found compelling evidence to support a general pattern such as that indicated in our initial transcriptomic measurements. in addition to our analyses of whole-fish mrna pools, we also confirmed tissue-specific expression patterns through q-pcr measurements of key immunity genes at a new estuarine locality. as indicated above we found many significant tissue-specific seasonal expression differences and all of these were in the direction predicted by our other whole-fish studies. the most pronounced seasonal expression profile occurred in the gill. furthermore, the gill contained the most intense concentration of t-cell activity outside of the thymus, with elevated expression of t-cell associated genes such as cd8a and foxp3b, and expression of these genes was seasonal in the gill but non-seasonal in the thymus. these observations are consistent with the known responsiveness of immune gene expression in the teleost gill to environmental stimuli, and also with the recent discovery and characterization of extensive, t-cell rich, interbranchial lymphoid tissue in teleost fishes. our results suggest the possibility that ilt may have an important role in seasonal immune function. finally, and whilst the present study is intended to characterize the seasonal dynamics of gene expression, rather than identify environmental causation, we briefly consider what external agents may drive the responses that we observed. in highly seasonal temperate zone habitats, such as the ones we consider here, each of temperature, diet, photoperiodic responses, pathogen exposures, or other biotic or abiotic manifestations of the environment, could be involved to unknown degrees. in the future, by matching detailed field observations with mesocosm studies and laboratory experiments, we expect to dissect the relative contributions of these influences to seasonal immune variation and to the immune phenotype more generally.  our results suggest that in wild teleosts, during winter conditions, adaptive immune activity declines in a manner that involves the expression of regulatory genes affecting lymphocyte function. this is indicative of a controlled, strategic response rather than a simple kinetic tracking of environmental temperature. seasonal change is most prominent in the gill, suggesting ilt may be important in such responses. further broad attention to seasonal immune function is certainly warranted, given the likely practical relevance – through effects on infectious disease susceptibility and inflammatory status – to health in humans and domesticated animals and to fitness in natural populations. 
these genes acted as hubs linking summer-biased adaptive pathways, winter-biased innate pathways and other organismal processes, including growth, metabolic dynamics and responses to stress and temperature. extremes of season were the dominant predictor of immune expression. further consideration should be given to the effects of such responses in seasonally-occurring disease. in contrast, late winter was accompanied by signatures of innate immunity and modulated toll-like receptor signalling. signatures of adaptive immunity were elevated in late summer. our results suggest major and predictable seasonal re-adjustments of immunity. fishes show seasonal patterns of immunity, but such phenomena are imperfectly understood in vertebrates generally, even in humans and mice. as these seasonal patterns may link to infectious disease risk and individual condition, the nature of their control has real practical implications. network analyses identified a small set of immune genes that might lie close to a regulatory axis. negative regulators of t-cell activity were prominent amongst winter-biased genes, suggesting that adaptive immunity is actively down-regulated during winter rather than passively tracking ambient temperature. seasonal change was most pronounced in the gill, which contains a considerable concentration of t-cell activity in the stickleback. we made genome-wide measurements of whole-fish mrna pools at the end of summer and winter in contrasting habitats and focussed on common trends to filter habitat-specific from overarching temporal responses. here we characterize seasonal dynamics in the expression of conserved vertebrate immunity genes in a naturally-occurring piscine model, the three-spined stickleback.
 fishes show seasonal patterns of immunity, but such phenomena are imperfectly understood in vertebrates generally, even in humans and mice. as these seasonal patterns may link to infectious disease risk and individual condition, the nature of their control has real practical implications. here we characterize seasonal dynamics in the expression of conserved vertebrate immunity genes in a naturally-occurring piscine model, the three-spined stickleback.  we made genome-wide measurements of whole-fish mrna pools at the end of summer and winter in contrasting habitats and focussed on common trends to filter habitat-specific from overarching temporal responses. we corroborated this analysis with targeted year-round whole-fish gene expression studies in a different year. we also considered seasonal tissue-specific expression at a third contrasting locality by q-pcr, further validating the generality of the patterns seen in whole fish analyses. extremes of season were the dominant predictor of immune expression. signatures of adaptive immunity were elevated in late summer. in contrast, late winter was accompanied by signatures of innate immunity and modulated toll-like receptor signalling. negative regulators of t-cell activity were prominent amongst winter-biased genes, suggesting that adaptive immunity is actively down-regulated during winter rather than passively tracking ambient temperature. network analyses identified a small set of immune genes that might lie close to a regulatory axis. these genes acted as hubs linking summer-biased adaptive pathways, winter-biased innate pathways and other organismal processes, including growth, metabolic dynamics and responses to stress and temperature. seasonal change was most pronounced in the gill, which contains a considerable concentration of t-cell activity in the stickleback.  our results suggest major and predictable seasonal re-adjustments of immunity. further consideration should be given to the effects of such responses in seasonally-occurring disease. electronic supplementary material the online version of this article contains supplementary material, which is available to authorized users. keywords seasonalityrnaseqteleostthree-spined sticklebackimmunityimmunoregulationecoimmunologywildlifehttp://dx.doi.org/ /501100000275leverhulme trustrpg-301issue-copyright-statement© the author 2016 