Partitioning the environmental drivers of immunocompetence

sistance and tolerance and, in turn, to the vital rates of infectious agents and the progression of the disease they cause. Our results illuminate the environmental regulation of vertebrate immunity (given the evolutionary con- servation of the molecular pathways involved) and they identify mechanisms through which immunocompetence and host-parasite dynamics might be impacted by changing environments. In particular, we predict a dominant sensitivity of immunocompetence and immunocompetence-driven host-pathogen dynamics to host diet shifts. The


H I G H L I G H T S
• Diet and temperature are the main drivers of immune allocation in a wild vertebrate. • Immune allocation corresponds to immunocompetence (driving infection dynamics). • Diet shifts will be the dominant driver of immunocompetence under climate change. • Epidemiological models should incorporate environmentally-driven immunocompetence.

G R A P H I C A L A B S T R A C T
a b s t r a c t a r t i c l e i n f o

Introduction
Variation in immunity influences individual propensity to disease (Parkin and Cohen, 2001;Rook and Dalgleish, 2011), contributing to health in humans, to productivity and welfare in domesticated animals, and to fitness in wild animals. Moreover, it constrains, and may be necessary to understand, the population dynamics of infectious disease (Hedrick, 2017). Whilst immune variation is driven largely by the environment (Brodin et al., 2015;Beura et al., 2016) (and there are numerous piecemeal observations of such effects (Bowden, 2008)), the relative contribution of different environmental components is remarkably poorly understood. This can be said for vertebrates generally, even humans and the best studied laboratory models (Jackson, 2015). Our aim here is to achieve a quantitative understanding of how the total environment controls immune allocation and immunocompetence (the ability to produce an effective immune response) in a natural system. To this end, we employ a naturally-occurring piscine model (three-spined stickleback, Gasterosteus aculeatus) (Jones et al., 2012;Stewart et al., 2017), disentangling the environmental drivers of immunocompetence, and the consequences of this for infectious disease, through a combination of within-habitat field observation, experimentation and modelling.
Our approach here has been to decompose within-habitat seasonal trends in immunocompetence, as a tractable surrogate for environmental effects in general. The rationale for this is that seasonal variation in environmental variables in (north temperate) stickleback habitats is typically large and easily measurable. Moreover, such seasonal variation inherently tends to affect all individuals within a population synchronouslyallowing powerful analysis of population-level trends. And because seasonal variation typically encompasses all of the major abiotic and biotic variables affecting animals (e.g., temperature, food, infection pressures, breeding cues, etc.), it is expected to be an excellent proxy for environmental drivers in general, including those drivers acting idiosyncratically on individuals to produce individualistic (Arriero et al., 2017;Tinsley et al., 2020) variation.
We have previously been able to identify a regular (sinusoid) circannual population-level trend in immune allocation in wild sticklebacks, partitioning this from background variation by statistical modelling Stewart et al., 2018b). Recognizing the severe inherent limitations in extracting causality from observational datasets, here we adopt an experimental approach (e.g., Stewart et al., 2018a) to disentangle the environmental drivers of this trend. Controlled mesocosm and laboratory experiments, using acclimatized wild hosts to most closely represent wild phenotypes, have been employed to parameterize models and predict the trend in observational time series in the field. Our study thus stands in contrast from much work in natural systems, through its ability to go beyond mere correlation and provide clear, quantitative, causal inference.
Crucially, we selected our measurements of the immune system (reflecting the seasonal trend) in a data-driven way from independent pilot genome-wide transcriptomic studies (Jackson, 2015;Brown et al., 2016). In contrast to the more common practice of arbitrarily selecting a limited panel of measurements based on a subjective preconception of their importance -where the amount of natural variation left unrepresented is unknown -our approach allows us to be sure that our measurements represent a large and defined component of natural variation. Thus, we have previously shown that, in wild sticklebacks, circannual seasonality represents the major axis of genome-wide variation in the expression of conserved immune-associated genes Stewart et al., 2018b). Furthermore, we can combine a small number of quantitative real-time PCR (QPCR) gene expression measurements, originally selected from our transcriptomic data, into a single index (Seasonal Reporter Index, SRI) (Stewart et al., 2018b;Friberg et al., 2019) that accurately tracks the circannual (sinusoidlike) oscillation in the immune system. Through using standardized assaying methods, we are able to employ SRI as a common currency (representing the circannual oscillation) across different studies, permitting the combination of field and experimental datasets to quantify the environmental drivers of SRI.
Very importantly, we have linked SRI to functional phenotypes. Below we show it to predict resistance and tolerance to the ecologically relevant (Stewart et al., 2017) monogenean parasite Gyrodactylus gasterostei, and we have previously shown SRI to also predict resistance to the oomycete pathogen Saprolegnia parasitica (Stewart et al., 2018b). Moreover, through new and previously published field and experimental data for G. gasterostei, we have parameterized models that further link SRI, and the environmental variables that drive it, to the withinhost disease progression and the population-level dynamics of this infection.
In the analyses below we draw on new and previously published field and experimental studies (Table A.1; Materials and methods). The dataset comprises 1409 individuals with SRI measurements and a further 460 wild fish with G. gasterostei infection data. This encompasses 46 monthly spatiotemporal samples in the wild derived from two ecologically distinct habitats across two annual cycles, 5 laboratory experiments and 3 mesocosm experiments. To combine the published and new data we revisit and re-analyse gene expression data from the earlier studies to represent SRI for the subset of genes (cd8a, foxp3b, ighm, orai1, tbk1) employed here, where high SRI represents high lymphocyte activity (Stewart et al., 2018b;Masud et al., 2019).
In summary, we have attempted to decompose the major (seasonal) axis of variation in immune allocation in wild sticklebacks, employing a combination of field observations, experimental manipulations of environmental variables and modelling. We have developed a simple gene expression index, SRI, that accurately reports this seasonal variation and have parameterized models for SRI variation in the field, based on multiple experimental estimates and extensive field data. Furthermore, through challenge infection experiments, we have linked SRI to an infection resistance phenotype and we have parameterized further models of parasite population dynamics, simulating the effects of change in the main environmental drivers of SRI. Taken together, our results demonstrate that thermal conditions and diet alone are necessary and sufficient to explain most variation in immunocompetence, which in turn greatly influences infection dynamics. Our findings also emphasize that, in practice, diet shifts are likely to have the largest potential influence on immunocompetence in changing environments.

Data sources and study sites
The experiments and observational field datasets upon which analyses of gene expression were based are summarized in Table A.1. In all experimental studies, wild fish were acclimatized to mesocosm or laboratory conditions following the application of anti-parasitic treatments, as previously described (Stewart et al., 2018a(Stewart et al., , 2018bFriberg et al., 2019;Masud et al., 2019). Most gene expression datasets were based on the analysis of whole-fish samples, but we also considered experiments based on gill samples, as gill undergoes similar SRI gene expression responses to environmental variation to those seen in whole-fish samples . We used fish from different localities (mostly peri-coastal Welsh sites) for experiments, as consistent seasonal environmental responses (in the SRI gene expression metrics we consider) occur across localities in England and Wales . Studies were carried out at one of two independent centres: either Cardiff University or Aberystwyth University. Importantly, we replicated the effects of the key experimental variables identified here (diet and temperature), in both whole-fish and gill samples, in fish from different localities, in mesocosm and laboratory experiments and in experiments executed at the two study centres. As previously described, the field study sites were an upland lake (Llyn Frongoch,FRN,52.3599,) and a lowland river (Afon Rheidol, RHD, 52.4052, −4.0372)  forming parts of different catchments in mid Wales.

Analysis software
All analyses were carried out in R version 3.5.2 (R Core Team, 2018).

Parasitology data
Previously unreported population data for Gyrodactylus gasterostei infection were collected for wild fish at FRN over a 25-month period, October 2013-October 2015, with missing monthly values at October 2014 and December 2014. Monthly samples (20 fish/month; total n = 460) were taken at the same times as the samples of fish taken for gene expression analyses (10 fish/month, see above, Table A.1). Fish within each monthly sample were captured randomly with a hand net, except for the final 5 months of the series, for which sampling was stratified across young-of-the-year (0 + ) and older fish (1 + ). Following sampling, fish were placed in plastic aquaria filled with lake water and returned to the laboratory, where they were killed by concussion and decerebration. All body surfaces were scanned under a dissecting microscope, whilst immersed in lake water, for the presence of G. gasterostei, which were individually counted. The two missing monthly data (October 2014 andDecember 2014) were included in the SIS model fitting procedures below as the average of the flanking values.
For all samples contributing data to the present study (see Table A.1), the body cavity of fish used for gene expression analyses was searched, under a dissecting microscope, for Schistocephalus solidus plerocercoids. This examination occurred following storage in RNA stabilization solution and prior to homogenization of the fish or fish tissues and nucleic acid extraction (as previously described; Hablutzel et al., 2016).

Assessment of breeding status
For all samples contributing data to the present study (see Table A.1), we checked the fish used for gene expression analysis for signs of reproductive maturity. Males were recorded as being in a breeding-related state ("breeding") if showing evidence of red breeding coloration. Females were recorded as breeding if a ripening or ripe ovary was present.

Gene expression measurements
Gene expression measurements from pre-existing datasets (see Table A.1) were made using standardized methods described in the original publications. Equivalent methods were used to generate the gene expression data for the new experiment described below (experiment 1). Briefly, RNA was extracted from whole fish samples preserved in RNA stabilization solution using the Isolate II RNA mini kit (Bioline). Whole individual fishes were homogenized in kit lysis buffer using a 5 mm stainless steel bead (Qiagen, 69989) in a Qiagen TissueLyser LT system and a standard aliquot of the homogenate passed through the manufacturer-recommended protocol. RNA extracts were DNAse treated and converted to cDNA using the High-Capacity RNA-to-cDNA™ Kit (ThermoFisher), according to manufacturer's instructions, including reverse transcription negative (RT-) controls for a subsample. Assays were pipetted onto 384 well plates by a robot (Pipetmax, Gilson) using a custom programme and run on a QuantStudio 6-flex Real-Time PCR System (ThermoFisher) at the machine manufacturers default realtime PCR cycling conditions. Reaction size was 10 μl, incorporating 1 μl of template and Applied Biosystems™ Fast SYBR™ Green Master Mix (ThermoFisher) and primers at the machine manufacturer's recommended concentrations. Samples from different experimental treatment groups were dispersed across 3 plates. Each plate contained all target gene expression assays and two endogenous control gene assays, for samples (in duplicate) and a calibrator sample (in triplicate). Endogenous control genes (yipf4, acvlr1) were previously validated , as a pairing, for stability under seasonal variation. Primers used were reported in Brown et al. (2016) and Hablutzel et al. (2016). In addition, no template controls for each gene were included on each plate. Template cDNA (see above) was diluted 1/20 prior to assay. The calibrator sample (identical on each plate) was created by pooling cDNA derived from whole fish RNA extracts from wild sticklebacks captured in summer. Relative gene expression values used in analyses are RQ values calculated by the QuantStudio 6-flex machine software according to the ΔΔCt method, indexed to the calibrator sample. Melting curves and amplification plots were individually inspected for each well replicate to confirm specific amplification.

Seasonal reporter index (SRI)
For the studies listed in Table A.1, we measured expression in a set of 5 genes (seasonal reporter, SR, genes) whose expression we have shown to reflect a major immune-system-wide seasonal expression signature Stewart et al., 2018b). These include three genes expressed highly in summer: cd8a (Ensembl gene identifier: ENSGACG00000008945), which codes for the CD8 T-cell receptor coreceptor in cytotoxic T cells, foxp3b (ENSGACG00000012777), which codes for a transcription factor that regulates the development and function of regulatory T cells, and ighm (ENSGACG00000012799), which codes for the constant region of the immunoglobulin heavy chain μ; and two genes expressed highly in winter: orai1 (ENSGACG00000011865), which codes for a calcium selective ion channel involved in T-cell activation and apoptosis, and tbk1 (ENSGACG00000000607), which codes for an enzyme with kinase activity involved in innate immune signaling pathways (see Brown et al., 2016). As previously described, we combined these gene expression measurements into an additive index (seasonal reporter index, SRI) that we have shown to reflect a major pattern of seasonality in the expression of stickleback immune-associated genes. For this, each raw relative gene expression variable was first log 10 transformed and standardized. The values for each gene variable were then summed, assigning negative or positive values to genes according to whether they were most expressed in winter (negative) or in summer (positive) in the transcriptomic study of Brown et al. (2016) . High values of SRI correspond to the high expression of the adaptive arm of the immune system .

SRI variability in field populations
To assess how well SRI represented genome-wide immune-associated gene expression in the transcriptomic (RNAseq) data of Brown et al. (2016) (expressed as log-transformed fragments per kilobase of exon per million fragments mapped, FPKM) we carried out a Mantel correlation of genome-wide immune-associated gene expression with SRI (based on Manhattan distance matrices). We also used a principal coordinates analysis (PCO) to represent genome-wide immune-associated gene expression, correlating the major axes with SRI (Pearson, r). Immune-associated genes in these analyses (n = 3648 genes) were defined as in Brown et al. (2016) based on the ImmPort comprehensive list of immune-related genes from the ImmPort database (Bhattacharya et al., 2014).
We analyzed SRI time series at FRN and RHD (October 2013 -October 2015) (F1-F4, Table A.1), employing the VariancePartition package (calcVarPart function) (Hoffman and Schadt, 2016) to calculate the proportion of variation associated with each explanatory term in a linear mixed model (LMM) fitted using the lmer function in lme4 (Bates et al., 2015). Explanatory terms included length (snout to end of tail in midline, i.e., fork length; mm), year (2013-2014/2014-2015), site, sex (male/female), breeding condition (breeding/non-breeding), S. solidus infection (infected/non-infected), assay plate and season. The season effect was represented by a sinusoid function of time (Stewart et al., 2018b) with separate slopes for each site × year combination.

Statistical modelling to assess effects on SRI in multiple experiments and field studies
To quantify effects that might be important drivers of seasonal SRI we re-analyzed previously published field, mesocosm and laboratory datasets, alongside previously unpublished data relating to these datasets (for parasite infection and host breeding condition) and the data from our new experiment (Table A.1). Field parasite records consisted of data for S. solidus (a large body cavity cestode) (Stewart et al., 2017) only, as the recovery of other parasites was incompatible with our methods for gene expression measurement (see Hablutzel et al., 2016). Overall prevalence for S. solidus in our field samples was 21%. Laboratory parasite records were for controlled challenge infections with the oomycete parasite Saprolegnia parasitica (Stewart et al., 2017) and G. gasterostei (Masud et al., 2019;present study). We also considered effects where subjects in mesocosm experiments remained infected with naturally-acquired S. solidus infections (that were refractory to the anti-parasite treatments we employed during acclimation of wild fish, see Hablutzel et al. (2016)); overall prevalence of S. solidus in mesocosm experiments was 23%. In datasets where some fish were sampled in breeding condition, this was also considered as an explanatory variable. Depending on which variables were represented in a given dataset, statistical models for SRI contained combinations of the following explanatory terms: sex (male/female), length (snout to end of tail in midline, i.e., fork length; mm), breeding condition (breeding/ non-breeding), S. solidus infection (infected/non-infected), G. gasterostei infection (challenged and infected/not challenged), Saprolegnia infection (overtly infected/non-infected), season (sinusoid function of time), temperature (°C) and other thermal treatments (factor), and diet (summer-like diet/winter-like diet).
For these analyses, the original raw datasets were stripped of any data rows with missing values (<0.5%, 7/1409; with respect to the sample sizes described in Table A.1) and the RQ values were converted to the scale of the calibrator sample used in Brown et al. (2016). The RQ values were then log-transformed (log 10 [x + 1]) and standardized (centered at the overall mean and divided by the overall standard deviation) and used to calculate SRI (as above). For simplicity and consistency each study was analyzed in a separate linear model (LM, lm function), not considering random terms for assay plate or other batch effects as we have not found these to be large or consistent effects on SRI previously. Due to the standard assay procedures employed, the parameters from different statistical models were then comparable on the same scale (Stewart et al., 2018b) and could also be used to parameterize the sinusoid SRI model below (see Table A.2).

A model for circannual SRI variation with sinusoid temperature and diet drivers
We modelled the sinusoid circannual variation in SRI, taking a similar approach to in (Stewart et al., 2018a), with two separate sinusoid functions (1.1), one representing seasonal thermal effects on SRI (thermal driver, D T ) and the other seasonal diet effects on SRI (diet driver, D D ).
D T depended on a cosine function (1.2), where t = time, τ = period and the parameters A T (amplitude) and Φ T (acrophase) were estimated through cosinor regression (Tong, 1976;Stolwijk et al., 1999;Cornelissen, 2014) of field temperature data. The output of this sinusoid function was then converted into SRI units by a coefficient (a) that was set as the thermal effect on SRI measured in different experiments manipulating temperature. For a description of how a was parameterized see Table A.2.
D D also depended on a cosine function (1.3), with parameters A D (amplitude) and Φ D (acrophase) that was converted into SRI units by a coefficient (d).
Parameterization of d. A D from experiments manipulating diet is described in Table A.2. The parameters x and Φ D were estimated through inverse fitting the partially parameterized overall function (1.1) to our field data for SRI using the FME package  (see Stewart et al., 2018a).

2.
10. An experiment (experiment 1) to calibrate diet-driven immunocompetence to intrinsic rate of increase in G. gasterostei 2.10.1. Rationale for experiment We carried out a challenge infection experiment, manipulating diet, to calibrate SRI variation to variation in an infection phenotype. For the infection system we employed the viviparous monogenean G. gasterostei (Stewart et al., 2017), a species that proliferates in situ on the host.

Hosts
Sticklebacks (n = 75) were caught using hand nets from Roath Brook stream (RBR), Cardiff, UK (51.499858°, −3.168780°) in January 2018 and transported to the aquarium at Cardiff University. As the source population was exposed to a natural community of parasites, the fish were subject to anti-parasitic treatment. All individuals were first submerged in 0.004% formaldehyde solution for 1 h, with a 0.5 h rest period in between, and then maintained in 1% aquarium salt and 0.002 g/l methylene blue for 48 h. Subsequently, each fish was visually screened (Schelkle et al., 2009) three times, involving anesthetization in 0.02% MS-222 and manual removal of any remaining gyrodactylids under a dissecting microscope with fibre-optic illumination. Following treatment, sticklebacks were maintained in 70l tanks with dechlorinated water (12 ± 1°C; 12 h light/12 h dark). All temperature and lighting conditions were maintained throughout the experimental trials.

Experimental diets and challenge infection
Two host dietary treatments were applied that aimed to mimic the extremes of seasonal diet quality observed in the wild (following Friberg et al. (2019)). A winter-like diet was made up of boiled, triturated spinach leaves, representing plant detritus, mixed with triturated gamma irradiated chironomid larvae (Tropical Marine Centre 8153) at a ratio of 2.5 g per 100 ml sieved detritus. A summer-like diet (identical to the summer-like diet of Friberg et al. (2019)) consisted of equal amounts of chironomid larvae, cyclopoid copepods (Tropical marine centre 8171) and cladocerans (Tropical marine centre 8151) (all gamma irradiated). Fish destined for the experiment were approximately size matched in pairs and the members of each pair were assigned randomly to different diet treatments. Each diet was provided in excess in each tank and was continuously available ad libitum. We note that wild stickleback may have lower food consumption in colder months of the year Wootton, 1983, 1984), and provision of unlimited access to food may thus underestimate the full adverse dietary effects of winter (i.e., diet effects could be even more substantial than established below). Nonetheless, we also note that the winter diet treatment in both the present feeding experiment and that of Friberg et al. (2019) produced body condition phenotypes like those of wild winter fish (see below), suggesting that these treatments were reasonably representative.
For each diet treatment, the fish were further (randomly) assigned to two groups. One group were infected with G. gasterostei after 3 weeks and the other group were sham-infected serving as negative controls. To perform controlled infections, heavily infected wildcaught donor fish were killed and the caudal fin was placed just touching an experimental fish until three worms were observed to transfer to the caudal fin of the recipient (donor-recipient contact lasting <5 min). This procedure was carried out using a dissecting microscope with fibre-optic illumination. All fish were screened the next day to ensure transferred worms were still present (this was the case for all fish). Infected fish and sham-infected control fish were isolated in 1 l pots and screened (King and Cable, 2007) every 4 days over a 25-day infection trajectory. Screening consisted of anesthetization in 0.02% MS-222 and enumeration of worms on each fish under a dissecting microscope with fibre-optic illumination.

Feeding response in experiment
The feeding response to the summer-like diet in the present experiment was not as strong as in the experiment reported by Friberg et al. (2019). In both experiments, fish entered the experiment with a winter-like phenotype (captured in January) and the winter-like diet sustained fish at a similar body condition (scaled mass index, SMI (Peig and Green, 2009)) to that in the wild . In contrast, the summer-like diet drove a + 86.9 mg SMI response in the experiment of Friberg et al. (2019) (M1, see Table A.1), compared to only a + 26.0 mg response in the present experiment (E5) (see Fig. A.1).

Analysis of experiment 1 data
To quantify the 25-day G. gasterostei infection trajectories on individual hosts in experiment 1, we estimated the intrinsic rate of increase, r (Birch, 1948), as the maximum value observed over any of the 4-day observation periods, calculated from the relationship: where N is the worm count, t is time in days and e is the base of the natural logarithm.
SRI and G. gasterostei r data from experiment 1 were analyzed, respectively, as the response in linear models (LMs; lm function) with sex (male/female), length (snout to end of tail in midline, i.e., fork length; mm), reproduction (breeding/non-breeding), assay plate and treatment group (winter-like diet/summer-like diet) as starting explanatory terms. LM results for SRI and r reported below are based on a minimal model selected by AICc (and in neither case were the results presented sensitive to different model selection procedures). Survival analysis was carried out employing a Cox proportional hazards model (coxph function) (Andersen and Gill, 1982;Therneau and Grambsch, 2000), assessing treatment group and length as explanatory variables.

A model for circannual r variation with sinusoid temperature and diet drivers
A sinusoid model similar to that constructed for SRI (Eq. (1.1), above) was constructed for G. gasterostei r, by combining Eqs. (2.4)-(2.6) below.

2.
12. An SIS model for G. gasterostei dynamics at FRN fitted to monthly prevalence data 2.12.1. Background We wished to develop a model that represented the dynamics of an ecologically relevant infection in terms of a transmission parameter with defined contributions from environmentally-determined host immunocompetence. Our aim was to parameterize this model using experimentally-derived estimates of the environmental effects and field observations and, ultimately, to simulate the effects of changing environmental conditions through perturbing the estimated model parameters.
We considered G. gasterostei as a model infection as this species was prevalent at our field study locality (FRN), had an existing background of information on its thermal responses (Harris, 1983) and was amenable to further laboratory experimentation. Moreover, although most helminth species at FRN showed limited seasonality in occurrence , seasonality was very marked in G. gasterostei, providing a strong signal to be explained (Fig. A.2). In our monthly field records for 2013-2015, the population dynamics of the species was dominated by a seasonal dilution effect due to host recruitment (Fig. A.2). Thus, there was a precipitous drop in prevalence and abundance to low levels as the (annual) host population turned over from 1 + to 0 + age cohorts as part of its usual circannual cycle (Wootton et al., 1978). Following host recruitment, prevalence (<1%) and abundance increased from low levels to a maximum (100% prevalence) in the month before the next host recruitment (Fig. A.2).
2.12.2. SIS model G. gasterostei, as for other Gyrodactylus species, typically produces self-limiting infections (Harris, 1983;Bakke et al., 2007) on isolated laboratory fish but stimulates only a transient refractory period in the host (Harris, 1983). This refractory period is likely even less important in the wild, where parasite emigration and host-to-host transmission may blur the infrapopulation kinetics observed in isolated laboratory hosts. Therefore, we chose an SIS model, with frequency-dependent transmission (Johnson et al., 2011), assuming any refractory period to be negligible and describing the numbers of susceptible (S) and infected (I) fish in the total population (N).
In this model (parameters and starting conditions defined in Table A.3) we included host population dynamics as a fixed daily recruitment during April-September (μ) and a fixed per capita mortality rate (ν). For a given starting October population size, to mirror the typical 0 + /1 + cohort turnover in the study population, we determined the ν that would deplete the 1 + cohort by the following September and μ that would restore the total October population to the same level as the previous year. We did not incorporate Gyrodactylus-induced host mortality as there is no evidence for this in wild stickleback. For the transmission function, we took the transmission rate β to be proportional to the parasite on-host intrinsic rate of increase, r, which was determined by additive contributions from a thermal effect (r T ) and a diet effect (r D ). The intrinsic rate of increase was then converted to the transmission rate via a coefficient, β′, that was estimated by fitting the model to field data (see below).
The respective thermal and diet contributions to r were derived from experimental data. In the case of r T we used data for thermal experiments on G. gasterostei published by Harris (1983), fitting a polynomial function to these (and forcing the intercept through zero). Given its seasonal nature, temperature (T) was represented in the model by a sinusoid function fitted to annual temperature data records at FRN (parameters defined in Table A.3).
Diet, likely affecting r through altering host immunocompetence, is similarly seasonal in nature and was thus also represented by a sinusoid function. To obtain the amplitude (A D ) of this function, we first estimated a diet effect on G. gasterostei r from experiment 1 in the present study. We then calibrated this effect to the observed SRI response to diet in the experiment and predicted the annual diet-driven range of r in the field as that corresponding to the diet-driven SRI range required, alongside thermal effects, to explain SRI circannual variation (based on analysis of Eq. (1.1), above). We set A D as half of this annual range. The mesor for the sinusoid function (M D ) was set at the same value as A D , zeroing the diet effect on r at the maximum food status (maximum immunocompetence). The acrophase (Φ D ) for the sinusoid function was determined by cosinor regression of the detritus index calculated by  based on the data of Allen and Wootton (1984) for FRN.
For the function representing temperature dependent per capita host recovery, we used data from (Harris, 1983), setting recovery below 5°C at zero and at 5°C and above based on a polynomial fit. Because these data were based on isolated laboratory hosts and might not reflect the situation in the wild, where parasites can continuously transmit from fish to fish without building up to high population sizes so rapidly, we included a coefficient (γ′) converting γ into the final recovery rate and estimated this by inverse fitting of the model. In the event, γ′ was small when fitted to field data, tending to confirm the above possibility and also supporting our earlier assumption to neglect any refractory period.

Fitting the model to field data
Partially parameterizing the SIS model as described above (summarized in Table A.3) we then fitted it to a field dataset for monthly G. gasterostei prevalence at FRN (R 2 = 86%; Fig. A.3, Table A.3) to estimate β′ and γ′. For this we used the R packages deSolve  (ode function) to solve models and FME  (modcost and modFit functions) to carry out constrained fitting to the data using a pseudorandom search algorithm.

Simulations based on the fitted model, varying temperature and diet conditions
We then took the fully parameterized model (including the fitted estimates of β′ and γ′) and varied the acrophase of the diet function to 180°in either direction (at which points it would be fully out-ofphase with the temperature function). Against this background we also varied the amplitude and mesor of the diet function. Amplitude was varied between zero and its original estimated value whilst maintaining the mesor equal to the amplitude (zeroing the diet effect on r at the highest diet-driven immunocompetence). Additionally, against a background of changing diet function amplitude (as above), we increased the temperature (mesor) between 0 and 2°C above the function based on field data (T). Model outputs (G. gasterostei prevalence, averaged across one year) were plotted using the R package plot3D (Surf3D function) (Soetaert, 2013).

Animal welfare
Experiment 1 was approved by the Cardiff University animal ethics committee and conducted under UK Home Office license PPL 303424.

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

Confirmation that SRI reflects a dominant, environmentally-driven, circannual axis of variation in immune allocation in wild sticklebacks
We have previously demonstrated, in sticklebacks from our main study sites (FRN and RHD), that seasonality is the major source of variation in the genome-wide expression of immune-associated genes, being of greater importance than host sex, ontogeny or habitat Friberg et al., 2019). The SRI measure we employ here strongly correlates with the genome-wide expression of immune-associated genes (Fig. 1a) in the transcriptomic dataset of Brown et al. (2016) (Mantel r = 0.46, based on 3648 immuneassociated genes). Furthermore, SRI undergoes a well-defined seasonal oscillation in our independent QPCR datasets for FRN and RHD Fig. 1. A simple gene expression index (Seasonal Reporter Index, SRI) reports a dominant seasonal genome-wide fluctuation in the expression of immune-associated genes. (a) Association of SRI with the first (PCO1) and second (PCO2) axes from a principal co-ordinates analysis of genome-wide immune-associated gene expression (n = 3648 genes); PCO1 and PCO2 respectively represent 43% and 14% of total variation; based on the transcriptomes of individual wild fish sampled in winter or summer at FRN and RHD (data from Brown et al., 2016); Pearson correlations (r) shown on panel. Overall Mantel correlation between genome-wide immune-associated gene expression and SRI was 0.46. (b) Non-parametric smoothers (from generalized additive models, GAMs) representing temporal variation in SRI, based on a 25-month time series at two localities, FRN (upland lake) and RHD (lowland river) (data from F1-F4, Table A.1). (c) SRI variance in wild sticklebacks at FRN and RHD partitioned by a linear mixed model (LMM) (data from F1-F4, Table A.1). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) (Fig. 1b). LMM analysis suggests that this seasonal oscillation accounts for 36% of total SRI variation, as much as the residual variation, and greatly exceeding other sources of variation (Fig. 1c). Interestingly, as habitat (FRN, lake vs RHD, stream) can be expected to be associated with genetic divergence in sticklebacks (e.g., Berner et al., 2009;Ravinet et al., 2013;Lohman et al., 2017;Huang et al., 2019;Rennison et al., 2019), as well as with habitat-specific environmental effects, its small effect on SRI suggests the genetic contribution to interindividual SRI variation in the wild is small (or possibly in opposition to habitat-specific environmental effects).

Diet and thermal variation are sufficient to explain circannual SRI oscillation in the field
Focussing on the SRI genes as a precise reporter of the seasonal (environmentally-driven) oscillation, we have previously been able to measure SRI responses to thermal and diet perturbations in multiple experiments (Table A.1), where higher temperatures (Stewart et al., 2018b) and summer-like (animal protein-rich) diets  drive high SRI. When we parameterized sinusoid models, containing separate sinusoid drivers for diet and temperature, using averaged estimates from these experiments (Materials and methods, Table A.2) and environmental observations in the field, we found that the combined effects of temperature and diet were sufficient to fully explain observed circannual SRI variation in the field, even at our site, FRN, with the highest observed circannual amplitude (Fig. 2a). At sites with lower observed circannual amplitude (RHD and mesocosms, see Stewart et al., 2018b) displacement of the diet driver peak towards the autumn, or reduced diet range, is required to prevent overestimation of the SRI oscillation (i.e., diet and thermal effects are even more than sufficient if in phase).

Diet and thermal variation are necessary to explain circannual SRI oscillation in the field
Further to the above, our data strongly indicate that diet and temperature are not only sufficient, but also necessary, for the circannual SRI oscillations in the field, as we can rule out other plausible drivers. In earlier experimental studies we have directly excluded the effects of ontogeny, photoperiodic or clock responses (Stewart et al., 2018b), or locomotory activity (Masud et al., 2019). Hydrochemical or external microbiological drivers are very unlikely explanations because they would be expected to show considerable habitat-specific temporal variation, whereas we have found broadly synchronous SRI oscillation across ecologically disparate sites (an upland lake, a lowland river and mesocosms filled from the tap-water supply) and modest betweensite differences Stewart et al., 2018b). Other important influences might include parasite exposures (Behnke et al., 1992) or effects due to breeding (Hanssen et al., 2005), and we are able to address these directly from our extensive experimental and field data (Fig. 2b). In the case of parasites, the effect of active S. solidus (larval cestode) (Stewart et al., 2017) or G. gasterostei infection on SRI is typically small or undetectable (Fig. 2b). The only large and consistent infection effect on SRI we have observed was due to overt S. parasitica infestation in laboratory experiments, but such phenotypes (representing an advanced acute infection that would result in death within hours or a few days in the laboratory) are rarely seen in the wild populations (<<1% prevalence in our datasets) and cannot account for observed seasonal SRI variation. Equally, through many observations on fish in breeding condition in the field and in multiple experiments, we find that the effect of this on SRI is generally negligible (Fig. 2b). Where we have observed a strong association between breeding and SRI in the field, this was only at one (riverine) site, suggesting the association was likely due to a habitat-specific environmental interaction. Furthermore, the direction of the latter effect (low SRI in breeding fish in the spring-summer breeding season, Fig. 2b) was inconsistent with the seasonal progression in SRI (rising in spring and peaking in summer). Thus, we can strongly infer that diet and temperature are the main drivers of circannual immune variation through the exclusion of other candidate drivers.

Calibration of SRI to an infection phenotype
We also evaluated SRI variation as a marker of functional immunocompetence, calibrating it to a relevant infection phenotype (G. gasterostei infestation (Stewart et al., 2017)). For this we fed acclimated wild fish winter-and summer-specific diets (following Friberg et al. (2019)), producing an uplift in SRI (+1.20 ± 0.04, P = 0.005) and  Table A.1). Significant experimental thermal effects, experimental diet effects and length effects are re-scaled to the seasonal range for these variables (see Table A.2). Significant breeding and infection effects are not rescaled so that they remain visible. body condition (scaled mass index (Peig and Green, 2009), +26.0 ± 10.6 mg, P = 0.018) in the summer-like diet group. We then challenged these fish with G. gasterostei (Experiment 1; details in Materials and methods). Parasite intrinsic rate of increase (r) was higher in the winter-like diet group (+0.10 ± 0.04 parasites/host/day, P = 0.038). Despite being unaffected by infection, SRI was a strong predictor of parasite r (−0.05 ± 0.016 parasites/host/day, P = 0.006, η 2 = 39%) (Fig. 3a) in infected fish, with diet treatment itself becoming a nonsignificant predictor when included in the same model. These results thus strongly support SRI as a marker of constitutive immunocompetence, as we have previously shown in S. parasitica infection (Stewart et al., 2018b). Strikingly, mortality was zero in control fish for both diet treatments, but when exposed to infection, fish were 10.2 × more likely to die if fed the winter diet (driving low SRI) rather than summer diet (Cox proportional hazards analysis, summer diet 2.368 ± 1.062, P = 0.004) (Fig. 3b). In terms of host response strategy, on-host r may be taken to reflect a resistance response (limiting the presence of parasite individuals) by the host and mortality to reflect a large component of tolerance (in the ecological sense (Ayres and Schneider, 2012;Jackson et al., 2014)). This experiment thus allows us to see that in its summer-like state (marked by high SRI) the stickleback's defensive system contributes very substantially to tolerance, as well as to resistance.

Immunocompetence and infection dynamics will be sensitive to shifts in seasonal thermal and diet profiles
Crucially, the large effect sizes and partly independent seasonal nature of temperature and diet mean that immunocompetence (SRI) and infectious disease progression (r) can be predicted to be extremely sensitive to shifts in their seasonal timings, as illustrated in Fig. 4a, b. In simulations with the sinusoid model derived above for SRI, in which the phase of seasonal diet variation is altered, seasonality in SRI is maximized under typical circannual conditions (high temperature and diet quality driving a peak in summer) (Fig. 4a). Conversely, if a similar model is constructed for parasite r (based on effects on r in our diet experiment and previously reported thermal effects (Harris, 1983)), seasonality in r tends to be minimized under typical conditions (r is promoted by high summer temperature but depressed by the higher quality summer diet) (Fig. 4b). Consequently, phase shift between temperature and diet (departing from the typical in-phase state, with peaks for both in summer) tends to reduce seasonality in immunocompetence (SRI) but increase seasonality in r, altering their seasonal relationship and the chance of seasonal outbreaks.
To further illustrate the epidemiological importance of these effects, we constructed an SIS (susceptibleinfectioussusceptible) model capturing the seasonal dynamics of G. gasterostei infection at our most seasonal field site (Fig. 4c, d; Materials and methods; Table A.3). These dynamics are consistent with a dominant circannual dilution effect ( Fig. A.2) (Harris, 1983) due to seasonal host recruitment superimposed upon a tendency for progressive accumulation of infection in the exposed population. The SIS model included fixed host recruitment dynamics and assumed the transmission (β) parameter would be proportional to the on-host parasite r (the ultimate determinant of the number of transmissible agents in the environment) which, in turn, would be driven by temperature and diet. Using parameter estimates on r from infection experiments manipulating temperature and diet, and converting the overall r estimate into a transmission rate via a coefficient β r , we then fitted the model to an observed prevalence time series (Fig. A.3), estimating β r (Table A.3). By altering the effect sizes on r within the parameterized model, and their circannual phases, we were then able to project the effect of diet, temperature and their relative seasonal timing on mean parasite prevalence. As for SRI and r, mean parasite prevalence is considerably sensitive to the seasonal timing of diet and temperature (Fig. 4c). In terms of longer term variation, moreover, the importance of diet effects are emphasised (Fig. 4c,  d), because where thermal variation might be limited by climatic constraints, fluctuation in food resources is not necessarily so constrained. For example, a difference in circannual mean (mesor) temperature of up to 2°C exerts a relatively small effect on infection burden compared to shifts in the diet within the natural circannual range (Fig. 4d).

Discussion
It is increasingly recognised that the environment drives much immunophenotypic variation (as supported by this study) but that we lack a systematic knowledge of what environmental drivers are quantitatively important (Jackson, 2015;Brodin et al., 2015;Beura et al., 2016). For the first time we have provided a quantitative and precisely qualified answer to this question: diet and thermal variation fully explain the dominant, seasonal environment-driven axis of immune variation in a wild vertebrate model.
To arrive at this result, we deliberately sought to disentangle and identify the environmental drivers of seasonal immunophenotypic variation in a highly seasonal, north temperate, system. We reasoned that these drivers would be a proxy for environmental drivers in general, given the tendency of most key environmental variables to vary seasonally in such a system. Moreover, a large seasonal environmental fluctuation would be expected to produce a clear, sinusoid-like seasonal immunophenotypic signal that could be explained analytically by the superimposition of sinuosoid functions for candidate environmental driver variables (Stewart et al., 2018a) (although this did not exclude the additional assessment of non-periodic drivers). We employed a gene expression index (SRI) that reflected a quantitatively dominant genome-wide immunological response to season in the wild and then we made controlled experimental estimates of different causal drivers of this index. We then sought to establish which of these drivers could quantitatively predict the SRI variation seen in replicated circannual field datasets. Only the effects of temperature and diet, in combination, were necessary and sufficient to explain the seasonal SRI pattern in the field. Although many environmental variables (e.g., infections) or intrinsic host variables might drive detectable signatures in the expression of individual immune-associated genes, our present results suggest these would be a small contribution in terms of the major axis of genome-wide variation seen in the wild populations over time.
Whilst there are naturally limits to how far this result can immediately be generalized, nonetheless the genes involved in seasonal immune expression in sticklebacks are highly conserved amongst all vertebrates , and thus might be controlled by conserved regulatory circuits. Furthermore, whilst G. aculeatus is ectothermic, and could be considered unrepresentative of endothermic organisms from the point of view of thermal biology, we note that dependencies on temperature are still likely relevant in endotherms. For example, endotherms may be unable to narrowly control temperature in peripheral tissues, affecting immunocompetence (Evans et al., 2015;Foxman et al., 2015), (a) Individual scatter of Gyrodactylus r against SRI. (b) Compared to a summer-specific diet, a winter-specific diet, that drove a decrease in SRI, resulted in a 10× increase in the probability of death following Gyrodactylus challenge. Predicted survival curves and 95% confidence intervals from a Cox Proportional Hazards analysis. and they might also use temperature as a strategic cue (Stewart et al., 2018b) or suffer trade-offs against the cost of heat generation (Cross et al., 2016).
Our present experiment (experiment 1) demonstrated that dietdriven components of immunophenotypic variation, reflected by SRI, influenced infection resistance and tolerance, affecting the vital rates of Gyrodactylus and the progression of the disease it causes (altering parasite r and parasite-induced host mortality). As for Gyrodactylus in the present study, we have previously demonstrated that increasing SRI constrains infection and disease progression in the oomycete pathogen Saprolegnia (Stewart et al., 2018a(Stewart et al., , 2018b. Whilst a vast additional range of pathogens exists in the environment, and each of these might respond differently to the same host immunophenotypic shifts, nonetheless our existing results, in quite disparate pathogens, suggest that widespread diet-driven perturbation of infection is possible. Importantly, the partially independent nature of the distinct thermal and dietary drivers we have identified entails that when these change phase, or become mismatched, there are likely to be major consequences for immunompetence and the dynamics of disease. This is important in the context of global climate change, and environmental change in general, where rapid ecological re-adjustment may promote temporal shifts in food availability (e.g., Burger et al., 2012;Doiron et al., 2015;Régnier et al., 2019). Indeed, whilst infectious disease outbreaks in wildlife are often attributed to a new parasite or pathogen strain, or to transmission phenomena, the size of the effects we have demonstrated suggest that disease outbreaks that are fundamentally due to environmentallydriven immunocompetence might also be frequent.
It is important to emphasize that our aim in this study has been to go beyond the mere description of seasonal patterns and to employ these patterns as a device to disentangle the environmental drivers of immunocompetence more generally. For the purposes of our study we deliberately sought out a system with intense circannual immunophenotypic variation. To produce such intense variation, it is unsurprising that the main environmental drivers would be seasonal and in-phase, as we, in fact, observed. Importantly, in other systems, exactly the same environmental drivers might be aperiodic, or out-of-phase, leading to a lack of immunophenotypic seasonality. This was illustrated for SRI by the simulations we carried out, parameterized from independent experiments, that varied the seasonal phase relationship of the diet and thermal drivers. In these simulations, seasonality in SRI is maximized under typical observed field conditions for our study site (high quality food and high temperature in summer) but declines as the food driver changes phase towards lower food quality in summer (Fig. 4a). It is important to note, therefore, that environmental drivers, such as we identified here, may operate just as influentially in systems with no or less seasonal immunophenotypic variation, if their respective temporal patterns do not happen to combine into a seasonal pattern.
Equally, parasite population dynamics may very well be aperiodic, or fluctuate apparently independently of variation in immunocompetence, even in the presence of dynamically important environmental drivers of immunocompetence. Often this might be due, for example, to the independent effects of temperature on parasite life processes (Stewart et al., 2018a). This was illustrated by our simulation of Gyrodactylus r (again, Fig. 4. Diet is expected to be a dominant driver of host immunocompetence during environmental change, exerting a large influence on within-host infectious disease progression and population-level infection dynamics. Upper panels show changes in circannual seasonal reporter index (SRI) (a) and Gyrodactylus intrinsic rate of increase, r (b) as the phase of seasonal diet oscillation (Phase Δ, radians) is varied from in-phase to 180°out-of-phase in either direction. Based on sinusoid models containing separate sinusoid driving functions for temperature and diet, parameterized from experimental results and field observations. Lower panels (c-d) show simulated average Gyrodactylus prevalence based on an SIS model parameterized from experimental results and field observations. In (c) the phase of the diet-driven oscillation is varied from in-phase to 180°out-of-phase in either direction, whilst the diet effect size is varied from a constant effect (maximum observed immunocompetence) to a seasonal effect like that observed in the field (progressing sinusoidally between the maximum and minimum observed immunocompetence). In (d) the diet effect size is varied as in (c), whilst the circannual mean temperature is increased across an interval of +2°C. For (c) and (d) the diet effect is determined by varying the amplitude (A) from 0.5× the estimated maximum range, to zero, with the mesor (M) set as M = A, so that decreasing A corresponds to increasing overall immunocompetence. parameterized from independent experiments) in which low seasonality in r is associated with typical observed field conditions (high quality food and high temperature in summer) (Fig. 4b). Here, positive thermal driving of r is counterbalanced by diet driving of SRI that tends to force r downward with high food quality. As the diet driver moves out of phase with the thermal driver, towards lower food quality in summer, however, the seasonality in r increases, as the seasonality in SRI decreases (see above). Aspects of transmission dynamics, for example the dramatic circannual dilution effect that we observed for Gyrodactylus, may also be superimposed upon and mask the effects of environmentally-driven immunocompetence. Nonetheless, important environmental drivers of immunocompetence may still exist, shifting the parameters of the overall host-parasite system from where they would otherwise be.
The latter point was illustrated by our simulations of Gyrodactylus population dynamics, employing an SIS model in which we included a transmission parameter that was a function of diet-and thermallydriven host immunocompetence. These simulations, parameterized using a combination of experimental and field observations, highlighted the large potential influence of environmentally-driven immunocompetence on levels of parasitism. This was even in a system where purely transmission phenomena (a dilution effect) generated profound circannual dynamics that did not map clearly onto circannual immunophenotypic variation. Furthermore, the outputs of the model provided some insight in terms of the likely sensitivity to thermal and diet change (Fig. 4c, d). Importantly, if thermal effects on immunocompetence primarily track the prevailing environmental temperature, as we have previously found in sticklebacks (Stewart et al., 2018a), this constrains them within certain limits and magnifies the potential importance of diet effects. For example, mean within-habitat temperatures are only likely to vary by one or two degrees at most year-on-year (O'Reilly, 2015;Sharma et al., 2015) and shorter-term fluctuations and shifts in circannual thermal profiles are also likely to be subject to significant climatic limits. Diet, on the other hand, might be even less constrained (potentially varying from plenty to famine) and might drive dynamics across a much wider outcome space. Notably, the availability of food organisms would be expected to respond to complex ecological drivers that might include the effect of thermal variation itself.
A further noteworthy feature of the present system, that may also be relevant in other systems where the thermal reaction norms for parasite host exploitation and host defence are positively correlated (Jackson and Tinsley, 2002;Stewart et al., 2018a), is that changing temperature may produce matched (compensating) changes in parasite virulence and host immunocompetence that wholly or partly cancel each other out. For example, both Gyrodactylus r (Harris, 1983) and stickleback immunocompetence (SRI) (Stewart et al., 2018a(Stewart et al., , 2018b increase with temperature, over a broad range, and may have antagonistic effects on disease outcome (Stewart et al., 2018a). In contrast, increasing diet quality might preponderantly increase immunocompetence; although a direct positive nutritional effect on the parasite is theoretically possible, there is relatively little evidence that this is important in vertebrate hosts (Pike et al., 2019). These fundamental constraints support that food availability, rather than temperature per se, will often be a key driver of immune-mediated climate change effects on host-parasite dynamics.
Taken together, our results argue for the incorporation of environmentally-driven immunocompetence into the parameters of epidemiological models, such as the transmission function β (McCallum et al., 2017), to attain a real understanding of the origins of infectious disease dynamics. Furthermore, given the large influence and relative simplicity of the key environmental drivers we have observed, this raises the prospect of better predicting disease dynamics through relatively simple, indirect environmental measurements (e.g., of temperature or food variability), and even of manipulating them through environmental interventions (e.g., food supplementation).

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.