- Research article
- Open Access
Changing expression of vertebrate immunity genes in an anthropogenic environment: a controlled experiment
BMC Evolutionary Biologyvolume 16, Article number: 175 (2016)
The effect of anthropogenic environments on the function of the vertebrate immune system is a problem of general importance. For example, it relates to the increasing rates of immunologically-based disease in modern human populations and to the desirability of identifying optimal immune function in domesticated animals. Despite this importance, our present understanding is compromised by a deficit of experimental studies that make adequately matched comparisons between wild and captive vertebrates.
We transferred post-larval fishes (three-spined sticklebacks), collected in the wild, to an anthropogenic (captive) environment. We then monitored, over 11 months, how the systemic expression of immunity genes changed in comparison to cohort-matched wild individuals in the originator population (total n = 299). We found that a range of innate (lyz, defbl2, il1r-like, tbk1) and adaptive (cd8a, igmh) immunity genes were up-regulated in captivity, accompanied by an increase in expression of the antioxidant enzyme, gpx4a. For some genes previously known to show seasonality in the wild, this appeared to be reduced in captive fishes. Captive fishes tended to express immunity genes, including igzh, foxp3b, lyz, defbl2, and il1r-like, more variably. Furthermore, although gene co-expression patterns (analyzed through gene-by-gene correlations and mutual information theory based networks) shared common structure in wild and captive fishes, there was also significant divergence. For one gene in particular, defbl2, high expression was associated with adverse health outcomes in captive fishes.
Taken together, these results demonstrate widespread regulatory changes in the immune system in captive populations, and that the expression of immunity genes is more constrained in the wild. An increase in constitutive systemic immune activity, such as we observed here, may alter the risk of immunopathology and contribute to variance in health in vertebrate populations exposed to anthropogenic environments.
During the transition between natural and anthropogenic environments the vertebrate immune system faces combinations of conditions unlike those it evolved to deal with. This is known to result in functional changes [1, 2] and what these changes are, and how and why they occur, is a key problem. There is a direct parallel to health in humans inhabiting relatively anthropogenic settings (for example, higher income countries), where an increasing burden of illness results from non-infectious diseases with inflammatory origins [3–5]. There is an equal relevance to domesticated animals or wild animals occupying urbanized habitats, where the immune system also functions (or malfunctions) under environmental conditions very different to those in nature.
Despite this background, research comparing immune function in the wild to in artificial habitats is in its early stages. The main body of existing work, comparing wild rodents with laboratory counterparts, suggests increased immunological activation in wild animals [1, 6–8], which may result from a greater exposure to infection in nature. Some responses to stimulation may be more intense and variably expressed in wild rodents [6, 7, 9], whilst other responses may be attenuated . Although these results are of great interest, the laboratory rodent models used as the basis for comparison bring with them aspects that may be unrepresentative of the real-world problem. Thus inbred mouse lines, unlike humans and domesticated animals, are genetically homogenous  and even outbred stocks may show restricted genetic variability . Furthermore laboratory rodents are maintained under extremely benign and pathogen-free conditions, whereas humans and domesticated animals still encounter many infections and environmental insults, albeit that these are different to those occurring in nature . Moreover, the singular genealogies of laboratory rodents make them difficult to compare directly with wild counterparts, even leaving the effects of inbreeding aside. Thus most laboratory stocks and lines have been in captivity for very many generations and are thus distant from the originator population, if this is identifiable at all. And they will often have been generated through arbitrary crosses [10, 11, 13, 14], resulting in haplotypes unrepresentative of those seen in nature. Hence complex genetic influences confound any comparison made to wild animals, leading to a basic lack of experimental control and uncertainty in interpretation.
In order to understand how loss of natural environment modifies the immune system it will be informative to study the immunophenotypic trajectory of wild animals newly acclimatized to anthropogenic conditions, with matched in situ controls in the wild [3, 12]. Importantly, this allows effects due to plasticity and to loss of on-going natural selection in the wild to be studied, unconfounded by long term effects of selection and breeding patterns within the anthropogenic environment. Whilst the latter processes are important in animal domestication, they are a separate issue that is not considered here. Furthermore, it should be noted that selection within the anthropogenic environment is unlikely to explain recent upwards trends in human immunopathologies, given their historical context [3, 12]. These are more likely driven by relatively recent environmental changes, making the focus of the present study of particular relevance.
To provide one case study of the type of acclimatization described above we focussed on the 3-spined stickleback (Gasterosteus aculeatus), a species that is accessible and much-studied in the wild , that easily acclimates to captivity, and that has an annotated whole genome , facilitating post-genomic studies. In the same way that other teleosts, such as zebrafish and medaka, are increasingly used to study disease processes relevant to mammalian health , the 3-spined stickleback – because it contains all of the central elements of adaptive immunity  - has a general comparative relevance for immunity in other vertebrates.
We transplanted post-larval fishes from a natural habitat to replicated artificial mesocosm habitats and, following anthelmintic treatment of the transplanted individuals, synchronously monitored both wild and transplanted (captive) cohorts through time. This study design embodies a general scenario typical of anthropogenic environments: where parasite exposure is reduced through anthelmintic treatment and curtailment of transmission , where bacterial exposures are altered due to artificial diets  and substrates , and where environmental stressors are different to in the wild due to plentiful food, altered density and social interaction, confined spatial ranges, absence of predation, and altered microclimate and chemical exposures. Our central aim here, though, is not necessarily to dissect the relative contributions of all these influences, but rather to generate a representative scenario and consider the immunological consequences and their health correlates.
As immunological readouts from our experiment we consider changes in the expression of a representative panel of conserved vertebrate immunity genes  in whole-fish mRNA pools . In using this whole-organism measurement approach we kept in mind the wide dispersal of the teleost immune system in different tissues  and aimed to achieve a holistic metric of immune activity – averaging across the entire immune system and all the tissues of the body. Such a metric is arguably more relevant to the general risk of systemic immunopathology, in comparison to a narrow, arbitrarily chosen focus on a single tissue or cell population within a tissue. Reductionist measurements of the latter type could be unrepresentative at the organism level (as so much is necessarily left unmeasured) and are in danger of reporting cellular trafficking between anatomical compartments, rather than overall levels of systemic activity.
This is the first study that we know of to carry out a closely matched immunological comparison of wild and captive vertebrates (i.e., where the individuals are ontogenetically matched and the results not clouded by differing genealogical histories in the study groups). We use our measurements to ask how do individual immune expression profiles vary between the wild and captivity and do certain immune expression profiles in captivity lead to adverse individual health outcomes?
Our experiment synchronously compared gene expression in wild lacustrine fishes with that in age-class-matched fishes from the same population that had been transferred to a representative anthropogenic habitat. This comparison was extended through time from the point when juvenile wild-caught fish had acclimatized to the anthropogenic habitat to a time when the wild cohort became scarce in nature.
For the anthropogenic habitat, twelve 300 L mesocosm tanks (52.4151°, −4.0670°), arranged in a 3 × 4 array, were evenly stocked with 480 postlarval sticklebacks taken from a lake habitat in mid Wales (52.3599°, −3.8773°) in July-August 2013. The mesocosms constituted two re-circulating systems, with 6 tanks in each. Within each re-circulating system water was re-circulated (3310 L h−1) by a pump (Blagdon, MDP3500) via a 90 L biological filter in a header tank. A manipulation of temperature was carried out within the mesocosm array that allowed the effects of thermal variation to be assessed independent of the comparison between wild and mesocosm habitats (as temperature variation might sometimes be confounded with these habitats due to natural climatic variation). Thus, one system was run at ambient temperature and the other heated to 2 °C above ambient temperature from the first experimental sampling point. Temperature in the heated system was maintained by 300 W heaters controlled by digital thermostats (1 per heated tank) with 0.1 °C sensitivity. For temperature control purposes each heated tank was paired to an adjacent unheated tank, with both providing thermistor feeds to the associated digital thermostat. Trials within the tank microenvironments showed that flow rates were sufficient to disperse temperature gradients around heaters within tanks. Tanks from each system were interspersed (alternating rows of 3 heated or 3 un-heated tanks across the array), to reduce positional effects as far as possible. Each tank contained standardized environmental enrichment (plastic aquarium plants) and a layer of light coloured gravel.
Following capture, wild postlarval fishes were subjected to 2 consecutive anthelmintic praziquantel treatments (24 h at 4 mg l−1; FlukeSolve, Fish Treatment Limited), separated by four days, following manufacturer’s recommendations. Of the common infections present in the wild population (based on 510 fish monitored between July 2013 and October 2015), these treatments completely removed Gyrodactylus sp. and a diplostomatid digenean species (the latter infecting the retinal layer of the eye), but had no measurable effect on Schistocephalus solidus worms or on ectocommensal trichodinids and epistylids. Prior to the commencement of the experiment in October, fishes were acclimatized for 4–6 weeks within the mesocosm system. Salinity in the system was routinely maintained at approximately 1 % (10 g L−1) as a prophylactic measure to suppress opportunistic microbial infections. Nitrite and nitrate levels (Tropic Marin Nitrite-Nitrate test) were continuously monitored throughout the experiment and remedial water changes carried out when nitrite levels rose above 0.02 mg L−1. Animals were fed daily on standard (per mesocosm) rations of frozen chironomid larvae (“bloodworm”), occasionally supplemented with frozen cladocerans (Tropical Marine Centre). Temperature in each tank was logged every 5 min, to a reading resolution ≤ 0.05 °C, throughout the experiment by Tinytag radio temperature loggers (TGRF-3024) networked through a Tinytag Radio system. In the field, temperatures were logged every 5 min by a Tinytag Aquatic 2 (TG-4100) data logger to a reading resolution ≤ 0.01 °C. Mean temperatures tended to be slightly higher (average difference < 1 °C) in the un-heated mesocosms than in the natural lake habitat (Additional file 1: Fig. S1a). Due to monthly sampling, individual density in the mesocosms fell from ~ 0.13 individuals L−1 at the start of the study period to ~ 0.028 individuals L−1 at the end. Reduction in individual density though, tended to be compensated by individual growth, with biomass density ranging within relatively narrow limits: from 0.013 to 0.027 g L−1 (Additional file 1: Fig. S1b). As the densities that we employed were relatively very low, they would have reduced the negative biological effects of crowding. At the same time, sufficient numbers of fish remained in each tank to allow these to aggregate in large groups and undergo social interactions. This scenario may have limited any influence of density variation on immunity (via crowding [24, 25] or social  mechanisms), although such effects cannot be eliminated. Notwithstanding, as density effects would not be biologically comparable in the field and mesocosms (wild individuals being able to range over much larger distances and interact with a much larger total population), we did not include density estimates in the analyses below. Instead we expected any between-habitat density-related differences to emerge in the habitat and habitat x time terms of statistical models. As considered above, altered density and constraints on movement and social interaction are just one of the many environmental components varying between natural and anthropogenic habitats - and it is beyond the scope of the present study to identify and fully dissect the influence of each of these components.
Animals were sampled monthly, on the same day, between October 2013 and August 2014: 10 animals month−1 in the originator lake habitat and 20 animals month−1 in the mesocosms (10 animals each from the heated and un-heated systems, drawing animals evenly from amongst different tanks). For the mesocosms, some additional monthly sampling (at the same sampling points as above) was carried out to provide (unused) spare capacity for the present study. Monthly samples from the wild and from mesocosms would have been approximately matched for age because they originated in the same annual recruitment cohort and this would have ceased to recruit new individuals soon after the last mesocosm stock were collected in the field. Although sticklebacks may breed across the spring and summer in the lake habitat, leading to a range of individual ages within the year cohort, this variation would likely have been distributed similarly within individual monthly samples (of 10 or 20 fish). Furthermore, the inclusion of length (a partial surrogate for age ) in statistical models (see below) will have additionally adjusted for age variation between individual fishes.
Animal handling and nucleic acids preparation
All animal maintenance and sampling of animals in the field followed U.K. Home Office regulations and local (Aberystwyth University) ethical approval procedures. Mesocosm work involved only manipulations fully compatible with routine captive fish husbandry (where the aim is to maintain a healthy stock) and was carried out in consultation with the HO inspectorate. Sampling occurred at standardized times of day, 09.00–13.00 h (UTC), and within a 2 h window at each sampling occasion. Sticklebacks were captured individually using a dip net and immediately killed by concussion and de-cerebration and stored in RNAlater™ RNA stabilization solution (ThermoFisher Scientific) according to manufacturer’s instructions (>5 volumes of stabilization solution; ≤ 0.5 cm solid tissue thickness along smallest dimension; a ventral incision was made to the abdominal cavity of each specimen to aid penetrance of stabilization solution). Samples were transferred to 4 °C overnight and then to −80 °C for long-term storage. Immediately prior to RNA extraction, sticklebacks were thawed at 4 °C, dabbed dry with tissue, and examined for Schistocephalus infection under a dissecting microscope (via a ventral incision). Fish weight (mg, minus the weight of any Schistocephalus infection) and standard length (mm) were recorded. RNA from whole fishes was extracted using the Isolate II RNA mini kit (Bioline): whole individual fishes were homogenized in lysis buffer using a 5 mm stainless steel bead (Qiagen, 69989) in a Qiagen TissueLyser LT system and a standard aliquot of the homogenate (diluted in lysis buffer to be equivalent to 11 mg of tissue) passed through the manufacturer-recommended protocol. Trials with whole fish preserved as above indicated that high RNA yields of good quality were consistently obtained from internal organs (spleen and head kidney), providing evidence of effective penetrance of RNA stabilization solution to the deep tissues.
Choice of target genes
Target genes were partly selected to represent a diversity of immunological pathways and partly to include those immune-associated genes determined to have high environmentally-related variability in a previous transcriptomic (RNAseq) study of sticklebacks in mid Wales . Adaptive responses were represented by the T-cell receptor co-receptor alpha chain from cytotoxic T-cells (cd8a) , the heavy chains of immunoglobulins M (igmh) and Z (igzh, see Gambón-Deza et al. ), the pro-inflammatory T-helper cell type 1 cytokine interleukin 12 (il12ba) , the pro-inflammatory T-helper cell type 17 cytokine interleukin 17 , the T-helper cell type 2 cytokine interleukin 4 (il4, see Ohtani et al., 2008 ), the regulatory T-cell transcription factor forkhead box P3 (foxp3b)  and a gene encoding a calcium channel known to be necessary for T-cell activation and proliferation in mammals (orai1) . Innate defences with direct activities against microbes, that may be expressed constitutively or inducibly following the stimulation of innate cells, were represented by the antimicrobial peptide beta defensin (defbl2,) [34, 35] and the bacteriolytic enzyme lysozyme C (lyz) . Innate inflammatory responses were represented by genes involved in toll-like receptor (TLR) (tirap, tbk1) [37, 38] and interleukin 1 (IL-1) family signalling (il1r-like) . We also measured the expression of the anti-oxidative enzyme, glutathione peroxidise 4a (gpx4a) , to reflect responses to oxidative stress that might drive, or be driven by, inflammation . Importantly, because of the elevated salt concentrations maintained in the mesocosm systems, none of these loci are present in a previously reported list of 1844 stickleback salt responsive genes .
Q-PCR gene expression measurements
Quantitative real-time PCR (Q-PCR) gene expression measurements were carried out in a 2-step format with SYBR Green chemistry. Primer sets (Table 1) used to assay target and endogenous control genes featured exon-spanning primers (except in the case of tirap, a single exon gene) and were determined to be 100 ± 10 % efficient under reaction conditions. Endogenous control genes, yipf4 and acvr1l, were previously selected as an optimally stable pairing for whole-fish analyses by the NormFinder algorithm  applied to RNAseq expression data for 36 wild fishes sampled at different sites in summer and winter . As for the target genes listed above, both endogenous control genes are absent from a reported list of stickleback salt-responsive genes . RNA samples were DNase treated and quantified spectrophotometrically (NanoDrop 2000) prior to conversion to cDNA with the High Capacity RNA-to-cDNA Kit (Life technologies) at the kit’s maximum capacity. A proportion of samples were also processed with reverse-transcriptase negative controls. Each cDNA product was diluted 1:20 prior to assaying. Samples were run in a 384-well plate format on a Life technologies QuantStudio 12 K flex real-time PCR system using SensiFAST SYBR® Lo-ROX mix (Bioline) and with the instrument manufacturer’s default fast cycle settings and recommended reaction concentrations scaled to 10 μl reactions. All assay plates were pipetted by an Eppendorf epMotion M5073 robot using a standard programme and layout. In addition to unknown samples (run in duplicate), each assay plate included a calibrator sample (run in triplicate) and no-template control wells (for every gene). For each sample, all target and endogenous control genes were run within a single plate and samples from different sampling units (time points and origins) were balanced across plates. Melting curves (run at instrument default settings) were examined for each assay well to assess the possibility of non-specific amplification. Relative expression (indexed to the calibrator sample and based on mean values for the technical replicates) was calculated using the ΔΔCt method. The calibrator sample was created by pooling cDNA aliquots from all individual samples.
Only wild fishes from the same year class as the mesocosm fishes were included in analyses. Our analyzed dataset consisted of 299 fishes, 84 from the wild and 215 from the mesocosms. Of these, 4.6 % overall (1 wild fish and 13 mesocosm fishes) had one or more missing measurements and thus are omitted from some or all analyses. An initial test of the hypothesis that habitat affected the expression of immunity genes was carried out using multivariate analysis of variance (MANOVA), with all immune gene expression variables as the responses, and with habitat, sex, Schistocephalus infection (present/absent) and time (month) as factors and body length, body condition and mean temperature in the 14 days before sampling as covariates. Inclusion of body length in the model allows for variation in immune expression due to age (for which body length is a substantial surrogate ) and general size, including the possibility of tissue allometry affecting whole fish measurements . Body condition was derived as the residuals from a quadratic regression of body weight on body length  (and there were no linear associations of the weight residuals with sampling time in either habitat, suggesting that the use of the residuals as a condition index was not biased by ontogenetic stage). Following a significant result in the MANOVA, the effects of habitat on the expression of individual immunity genes, and also gpx4a, were considered using general linear models (LMs) of analogous structure to the MANOVA (and allowing variation in the set of explanatory terms used for each gene expression response; see next). Because temperature was partially confounded with habitat (Additional file 1: Fig. S1a), and did not affect many of the genes examined across the relatively large (2 °C) thermal manipulation in the mesocosms, the MANOVA described above may have been conservative in detecting habitat effects. Thus for the LMs we only included a temperature covariate for those genes where we found a significant effect of temperature in the mesocosms. The effect of Schistocephalus infection and temperature within the mesocosm environment was tested in LMs containing factors for sex, Schistocephalus infection, time and temperature treatment (ambient/+ 2 °C) and covariates for body length and body condition. The association of individual gene expression variables with body condition (weight residuals, see above) in mesocosm fishes was also tested in LMs, in this case additionally containing factors for sex, Schistocephalus infection, time and temperature treatment. In all cases gene expression data were log10 transformed before analysis by MANOVA or LM. Principal components analysis (PCA) was used to ordinate gene expression measurements from individual fishes in order to identify habitat-specific patterns of variation. Differences in the overall variability of untransformed individual gene expression variables between the field and the mesocosms were tested using Levene’s test for equality of variances. MANOVA, LMs, PCA and Levene’s test were implemented in Minitab version 16.2.2. For co-expression analyses we used gene expression residuals from LMs (with terms for sex, length, body condition and Schistocephalus infection) in order to adjust for confounding variation. We initially constructed Pearson correlation matrices from the residuals and tested for significant structure in individual matrices (null = an identity matrix) using a Steiger test , and for differences between matrices with a Jennrich test  (using the R package psych). We tested for similarity in structure between matrices with a Mantel test [47, 48] (using the R package ape). We also used the residuals to construct co-expression networks with the information-theory (mutual information, MI) based algorithm ARACNe2 (Algorithm for the construction of accurate cellular networks) [49, 50], which takes non-linear associations into account. For this analysis we set all genes as hubs and constructed networks using the adaptive partitioning algorithm, estimating the MI threshold by a pre-processing run. P threshold was set at 1 × 10−7 and the data processing inequality at zero. Networks were bootstrapped (2000 resamples; significance cut-off for reported edges, P = 1.0 × 10−6). Cytoscape 2.8 was employed to visualize the networks (using arbitrary force-directed layouts) and to calculate network intersections (Network Analyzer plugin).
A range of innate and adaptive immunity genes are up-regulated in captivity and this is, in part, associated with a loss of seasonality
Fishes reared in captivity or taken from the wild showed similar, and extensively overlapping, biometric characteristics (Additional file 1: Fig. S1c, d). An initial overall confounder-adjusted MANOVA test for expression in 13 immunity genes revealed a highly significant difference (main effect) between wild and mesocosm fishes (F 14,256 = 3.72, P <0.0005). To dissect this result further we then examined each of the individual variables in confounder-adjusted general linear models (LMs) (Fig. 1a; Additional file 2: Table S1), finding significant main effects for igmh, cd8a, defbl2, lyz, il1r-like and tbk1. In all of these cases, expression was very significantly up-regulated in the mesocosms. We also considered the antioxidant enzyme, gpx4a, and found that this too was very strongly up-regulated in the mesocosm fishes. Furthermore, this gene co-varied positively with most immunity genes, consistent with a link to oxidative stress that might drive, or be driven by, immune responses. Thus log-transformed gpx4a was positively associated with the major axis of a log scale PCA of immunity genes in both habitats (lake r = 0.59, P <0.0005; mesocosms r = 0.55, P <0.0001), where this axis predominantly represented positive covariation, and it was individually positively associated (P <0.05) with 6 out of 13 log-transformed immunity genes in lake fishes and 12 out of 13 genes in mesocosm fishes.
For genes with significant main effects for habitat we also examined all two-way interactions involving habitat (except for the interaction with temperature, where this was included in the model). This revealed that in most cases there was a complex interaction with time (Fig. 1b; Additional file 2: Table S1), though the form of his interaction was not consistent across genes, with relatively elevated expression in the mesocosms occurring in the autumn and winter (igmh, cd8a, lyz), the winter and spring (defbl2) or without a clear pattern (il1r-like). In some cases, if data for July and August are put aside due to the small number of wild fishes sampled, peak average expression in the mesocosms reached higher than it ever did in the wild (lyz, il1r-like, defbl2, Fig. 1b). On the other hand, for cd8a and igmh (genes previously known to show summer-biased expression in the wild ) mean expression was more consistent in mesocosm fishes, diminishing any seasonal pattern. Thus similar high expression levels of cd8a and igmh occurred in the mesocosms and the wild during summer, but mesocosm fishes maintained high expression during winter periods when this was reduced in wild fishes. Additional to the interactions with time, there was an interaction of habitat with length (F 1,259 = 10.09, P = 0.002) for defbl2, with a trend for increased expression in smaller mesocosm fishes.
Captive animals have altered patterns of gene co-expression
As the regulatory influence of genes on other genes may be reflected in patterns of co-expression, we searched for differences in such patterns that might reveal regulatory changes in the immune systems of wild and captive fishes. We based this analysis on fishes from un-heated mesocosms, so that temperature conditions were similar between the two samples (<1 °C average difference) and there were approximately equal sample sizes in the two groups. In order to remove across-site biases due to individual host variables, we took residuals from LMs including terms for sex, length, body condition and Schistocephalus infection.
Initially considering pair-wise (residual) gene expression correlation matrices for the two sites, we found that both were strongly structured (Steiger test; wild, P = 2.7 × 10−62; mesocosm, P = 1.2 × 10−79). The matrices differed significantly (Jennrich test, P = 6.89 × 10−13), but retained some similarity (Mantel test, P = 0.001) (Fig. 2a), with correlation coefficients from wild fishes explaining 27 % of the total variation in coefficients from captive fishes in a linear regression (Fig. 2b).
We then compared (residual) gene co-expression networks for wild and captive fish constructed using the information theory-based algorithm ARACNe, which takes into account non-linear associations between the individual variables (Fig. 2c ). Some edges (statistical links between genes) were identical in both networks (43 % shared in wild fishes and 38 % in mesocosm fishes) but there were also a large number of differences. These included the absence of orai1, a gene predicted to be involved in seasonal immune function , from the mesocosm network. Also the antioxidant enzyme gpx4a showed co-expression with more genes in the mesocosm network.
Only defbl2 expression was associated with individual condition in captive fishes
We considered how the genes up-regulated in captivity might be associated with individual condition (weight residual, see Material and methods) specifically within the mescosm environment. We found only a (negative) effect for defbl2 (F 1,192 = 9.49, P = 0.002).
The expression of immunity genes is more variable in captive fishes
We also considered the comparative amount of variability in gene expression in the wild and in captivity. Individual scores along the major axis of variation in a PCA of raw (untransformed) data, including all specimens from both habitats, were more variable in the mesocosms (Fig. 3a) (Levene’s test, P = 0.024). This major axis (PC1) accounted for 42 % of total variation and principally represented positive covariation in gene expression (i.e., with component loadings mainly of the same sign; see Fig. 3a). Mesocosm fishes extended much farther than wild fish along PC1 in the direction of higher expression values (but not in the other direction), indicating that fishes with high levels of overall immune activity were responsible for the greater variability in the mesocosms. A similar result was obtained (with greater statistical significance), if only individuals from un-heated mesocosm tanks were included in the comparison (P = 0.016). For individual genes (Fig. 3b), untransformed expression of igzh, foxp3b, lyz, defbl2, il1r-like and gpx4a was significantly more variable in the mesocosms, and this remained the case if only un-heated tanks were included in the comparison (Additional file 3: Table S2). Across this set of heteroscedastic genes (especially dfbl2 and igzh), considerable numbers of mesocosm fishes demonstrated expression values well in excess of that ever seen in wild fishes (Fig. 3b).
There is no effect of Schistocephalus infection on genes up-regulated in the mesocosms
Finally we asked whether a reduced burden of macroparasitic infections might contribute to the immunological changes seen (above) in the mesocosms. Although differences in infection pressures were generally confounded with other mesocosm characteristics, we were able to cleanly assess the effect of one infection, Schistocephalus, amongst mesocosm fishes. This infection was refractory to the anthelmintic treatment used at the start of the study and infected a substantial proportion of mesocosm fishes (36 %), providing the basis for comparison on a level playing field. We found that there was no effect of Schistocephalus infection on any of the genes up-regulated in captivity.
We found significant up-regulation of aspects of both innate and adaptive immune gene expression in captive (mesocosm) compared to wild fishes. When host and environmental confounders were accounted for, 6 out of 13 (46 %) of the immunity genes examined were differentially expressed, in all cases up-regulated in the mesocosm environment. Of the six up-regulated genes, two are adaptive immunity genes pivotal in cytotoxic T-cell (cd8a) and antibody (igmh) responses. Of the other up-regulated genes, the bacteriolytic enzyme lysozyme (lyz) and the antimicrobial peptide beta defensin (defbl2), code for molecules with direct antimicrobial activity and that may be deployed by the host as standing or inducible innate defences. The remaining up-regulated genes (il1r-like, tbk1) are components of signalling pathways involved in innate inflammation. Although we did not directly observe how these gene expression responses corresponded to functional immune responses, their relevance to functional immune status is supported by the frequently observed correlation between gene expression and infection in sticklebacks [51–54] and many other vertebrates. Furthermore, although complex post-transcriptional kinetics might occur for some individual genes that lead to mRNA and bioactive protein responses of differing sign , the concordance in the direction of the above responses is indicative of a broad up-regulation of diverse immune pathways in captivity.
For some of the genes up-regulated in captivity, there was evidence that this was due to a diminution of seasonality. Thus, for two adaptive immunity genes previously known to show summer-biased expression in the wild (cd8a, igmh) , there was a time × habitat interaction in statistical analyses. The form of this interaction was such that maximum average expression values were similar in both habitats, but captive fish tended to maintain high expression during winter periods when expression was seasonally reduced in the wild. In these cases more challenging environmental conditions in winter that lead to a reduction of adaptive immune activity may be the origin of the overall difference . In contrast, for other genes (defbl2, lyz, il1r-like), differing temporal responses occurred that involved higher average expression in captive fishes than was ever recorded in wild fishes. Importantly, as such responses may extend (in some individuals) outside the range seen in nature, and thus outside the parameter space previously operated upon by natural selection, they may uncover maladaptations in immunoregulatory networks and increase the risk of immunopathology .
Analysis of inter-individual variability confirmed that the expression of immunity genes was generally less constrained in captive than in wild fishes. A number of genes showed significant differences in the variability of expression in the wild and captivity (igzh, foxp3b, lyz, defbl2, il1r-like and gpx4a), and in all cases these were more variable in captivity. This increased variability was largely due to some captive individuals with high expression levels. These observations are consistent with a relaxed control of immunological activity in captive fishes, perhaps due to a lack of strong and overriding environmental pressures that tend to affect all individuals in the population (for example seasonal change, or greater motor activity required in foraging and predator avoidance).
When we considered the pattern of co-expression amongst genes, via gene-by-gene correlation or information theory based networks, our analyses revealed a similar core structure in the wild and in captivity, but also significant divergence between the two. For example, orai1 - a gene we have already implicated in seasonal immune function  - was (statistically) linked to other genes within the co-expression network for wild fishes, but not within the network for mesocosm fishes. This further emphasizes a disruption of seasonal responses (see above) as one of the consequences of captivity. Interestingly, expression of the antioxidant enzyme, glutathione peroxidise 4 (gpx4a)  tended to correlate positively with that of immunity genes in both captive and wild fishes, and was very strongly up-regulated in captive fishes in comparison to wild fishes, suggesting a propensity to track the general level of immune activity. In networks and correlation analyses this gene demonstrated more statistical links to immunity genes in captive than in wild fishes. Taken together, these observations are consistent with captive fishes generating greater anti-oxidant activity to counteract oxidative pressures that might drive, and be driven by, increased inflammation . Alternatively, there may be mechanisms that limit immunological activity according to anti-oxidative capacity  (i.e., where anti-oxidative potential is low immune responses are curtailed in order to prevent oxidative self-damage). In this latter regard, one possible link to environmental variation might be through selenium intake in the diet, as this element can be present in limiting quantities in nature and is a cofactor for selenoprotein enzymes (including glutathione peroxidases) involved in anti-oxidative physiology [58, 59].
We also considered whether genes that were up-regulated in captivity were associated with the health status of individual captive fishes, using body condition (weight adjusted for length) as a proxy for health. We found that only defbl2 was associated, negatively, with body condition. Interestingly this gene also tended to be expressed more strongly in smaller captive fishes (adjusting for body condition), but did not show a similar trend in wild fishes, consistent with a link to growth patterns in captivity. As maximum individual expression levels of defbl2 in captivity extended well outside those ever seen in the wild, this may be an example of a maladaptive immune response - previously unrefined by natural selection - that results in immunopathology. It is worth noting, furthermore, that high beta defensin gene copy number (which can increase beta defensin expression levels) is known to be associated with immunologically based disease in humans (Chron’s disease and psoriasis) .
Our design compared the overall effects of a “complex” of environmental changes, associated with transfer to artificial environments, to the situation in the wild. Although the loss of most parasites was confounded with these general environmental changes, the refractoriness of Schistocephalus to anthelmintic treatment meant that the effect of this infection could be considered against the background of domestication conditions. Schistocephalus infection, though, did not explain variation in the expression of the immunity genes that were up-regulated in captivity. Thus, if parasite loss is an important agent of immunological change in captivity, then this might only be through a subset of the naturally-occurring assemblage not including Schistocephalus. Interestingly, Schistocephalus, which occupies a sterile site of infection (the abdominal cavity), may not be able to influence the immune system indirectly through interactions with symbiotic microbiota .
A further possibility that should be considered is that selection in the wild resulted in some of the expression differences that we observed. Thus the mesocosms, where mortality was negligible during the course of the experiments, would have experienced a very different selective landscape to in the wild, where there is steady attrition of year cohorts. The present results could in part be consistent, for example, with selection in the wild purging individuals whose expression profiles drift outside a given parameter space. On the other hand, cryptic selection at the point of stocking mesocosm fishes is an unlikely explanation for the patterns that we saw because, rather than occupying a subset of the phenotypic space spanned by wild fishes, mesocosm fishes in fact occupy a much wider phenotypic space.
In summary, we have conducted a large, carefully controlled experiment comparing immune system expression in the wild and in a representative anthropogenic setting, using ontogenetically synchronized subjects with matched genealogical and environmental origins. Adjusting for the effects of confounder variables, a broad set of vertebrate immunity genes, representing diverse pathways, were up-regulated and expressed more variably in populations in the anthropogenic environment. This was accompanied by changes in patterns of co-expression and an increased signature of anti-oxidative activity (expression of gpx4a). The latter trend and the generally strong patterns of co-expression involving gpx4a, suggest that immune activity levels are tightly linked to oxidative status. It should be noted that observed changes from the wild state were coincident with the complex of environmental shifts and the altered selective landscape experienced by individuals acclimatising to the anthropogenic environment. However, it was possible to rule out immunological effects of the larval pseudophyllidean cestode Schistocephalus as a source of changes and to pinpoint a dampening of seasonal responses as an important effect for some genes. Taken together, our results demonstrate widespread regulatory modulation across the immune system of individuals maintained in an anthropogenic environment. Importantly, our whole-organism measurement approach allows us to infer an increase in systemic immunological activity – extending, in some individuals, into areas of phenotypic space that are rarely occupied by wild animals. The occupation of such exotic phenotypic spaces (as, for example, in the case of beta defensin) may uncover features in immunoregulatory networks that have never been refined by selection and that may result in immunopathology.
General linear model
Multivariate analysis of variance
Principal components analysis
Quantitative real-time PCR
High throughput RNA sequencing
Devalapalli AP, Lesher A, Shieh K, Solow JS, Everett ML, et al. Increased levels of IgE and autoreactive, polyreactive IgG in wild rodents: implications for the hygiene hypothesis. Scand J Immunol. 2006;64:125–36.
Ledin A, Arnemo JM, Liberg O, Hellman L. High plasma IgE levels within the Scandinavian wolf population, and its implications for mammalian IgE homeostasis. Mol Immunol. 2008;45:1976–80.
Friberg IM, Bradley JE, Jackson JA. Macroparasites, innate immunity and immunoregulation: developing natural models. Trends Parasitol. 2010;26:540–9.
Jackson JA, Friberg IM, Little S, Bradley JE. Review series on helminths, immune modulation and the hygiene hypothesis: immunity against helminths and immunological phenomena in modern human populations: coevolutionary legacies? Immunology. 2009;126:18–27.
Rook GA. Regulation of the immune system by biodiversity from the natural environment: an ecosystem service essential to health. Proc Natl Acad Sci U S A. 2013;110:18360–7.
Abolins SR, Pocock MJ, Hafalla JC, Riley EM, Viney ME. Measures of immune function of wild mice, Mus musculus. Mol Ecol. 2011;20:881–92.
Boysen P, Eide DM, Storset AK. Natural killer cells in free-living Mus musculus have a primed phenotype. Mol Ecol. 2011;20:5103–10.
Trama AM, Holzknecht ZE, Thomas AD, Su K-Y, Lee SM, et al. Lymphocyte phenotypes in wild-caught rats suggest potential mechanisms underlying increased immune sensitivity in post-industrial environments. Cell Mol Immunol. 2012;9:163–74.
Lesher A, Li B, Whitt P, Newton N, Devalapalli AP, et al. Increased IL-4 production and attenuated proliferative and pro-inflammatory responses of splenocytes from wild-caught rats (Rattus norvegicus). Immunol Cell Biol. 2006;84:374–82.
Beck JA, Lloyd S, Hafezparast M, Lennon-Pierce M, Eppig JT, et al. Genealogies of mouse inbred strains. Nat Genet. 2000;24:23–5.
Chia R, Achilli F, Festing MFW, Fisher EMC. The origins and uses of mouse outbred stocks. Nat Genet. 2005;37:1181–6.
Jackson JA. Immunology in wild nonmodel rodents: an ecological context for studies of health and disease. Parasite Immunol. 2015;37:220–32.
Sakai T, Kikkawa Y, Miura I, Inoue T, Moriwaki K, et al. Origins of mouse inbred strains deduced from whole-genome scanning by polymorphic microsatellite loci. Mamm Genome. 2005;16:11–9.
Takada T, Ebata T, Noguchi H, Keane TM, Adams DJ, et al. The ancestor of extant Japanese fancy mice contributed to the mosaic genomes of classical inbred strains. Genome Res. 2013;23:1329–38.
Wootton RJ. The Darwinian stickleback Gasterosteus aculeatus: a history of evolutionary studies. J Fish Biol. 2009;75:1919–42.
Jones FC, Grabherr MG, Chan YF, Russell P, Mauceli E, et al. The genomic basis of adaptive evolution in threespine sticklebacks. Nature. 2012;484:55–61.
Schartl M. Beyond the zebrafish: diverse fish species for modeling human disease. Dis Model Mech. 2014;7:181–92.
Boehm T, Swann JB. Origin and evolution of adaptive immunity. Annu Rev Anim Biosci. 2014;2:259–83.
Borkow G, Leng Q, Weisman Z, Stein M, Galai N, et al. Chronic immune activation associated with intestinal helminth infections results in impaired signal transduction and anergy. J Clin Invest. 2000;106:1053–60.
Wang J, Linnenbrink M, Künzel S, Fernandes R, Nadeau M-J, et al. Dietary history contributes to enterotype-like clustering and functional metagenomic content in the intestinal microbiome of wild mice. Proc Natl Acad Sci U S A. 2014;111:E2703–10.
Hanski I, von Hertzen L, Fyhrquist N, Koskinen K, Torppa K, et al. Environmental biodiversity, human microbiota, and allergy are interrelated. Proc Natl Acad Sci U S A. 2012;109:8334–9.
Jackson JA, Begon M, Birtles R, Paterson S, Friberg IM, et al. The analysis of immunological profiles in wild animals: a case study on immunodynamics in the field vole, Microtus agrestis. Mol Ecol. 2011;20:893–909.
Brown M, Hablützel P, Friberg IM, Thomason AG, Stewart A, et al. Seasonal immunoregulation in a naturally-occurring vertebrate. BMC Genomics. 2016;17:1–18.
Jia R, Liu B-L, Feng W-R, Han C, Huang B, et al. Stress and immune responses in skin of turbot (Scophthalmus maximus) under different stocking densities. Fish Shellfish Immunol. 2016;55:131–9.
Yarahmadi P, Miandare HK, Fayaz S, Caipang CMA. Increased stocking density causes changes in expression of selected stress- and immune-related genes, humoral innate immune parameters and stress responses of rainbow trout (Oncorhynchus mykiss). Fish Shellfish Immunol. 2016;48:43–53.
Filby AL, Paull GC, Bartlett EJ, Van Look KJW, Tyler CR. Physiological and health consequences of social status in zebrafish (Danio rerio). Physiol Behav. 2010;101:576–87.
Takizawa F, Dijkstra JM, Kotterba P, Korytář TT, Kock H, et al. The expression of CD8α discriminates distinct T cell subsets in teleost fish. Dev Comp Immunol. 2011;35:752–63.
Gambón-Deza F, Sánchez-Espinel C, Magadán-Mompó S. Presence of an unique IgT on the IGH locus in three-spined stickleback fish (Gasterosteus aculeatus) and the very recent generation of a repertoire of VH genes. Dev Comp Immunol. 2010;34:114–22.
Zhang L, Zhang BC, Hu YH. Rock bream (Oplegnathus fasciatus) IL-12p40: identification, expression, and effect on bacterial infection. Fish Shellfish Immunol. 2014;39:312–20.
Zou J, Secombes C. The function of fish cytokines. Biology. 2016;5:23.
Ohtani M, Hayashi N, Hashimoto K, Nakanishi T, Dijkstra J. Comprehensive clarification of two paralogous interleukin 4/13 loci in teleost fish. Immunogenetics. 2008;60:383–97.
Wen Y, Fang W, Xiang LX, Pan RL, Shao JZ. Identification of Treg-like cells in Tetraodon: insight into the origin of regulatory T subsets during early vertebrate evolution. Cell Mol Life Sci. 2011;68:2615–26.
Qu B, Al-Ansary D, Kummerow C, Hoth M, Schwarz EC. ORAI-mediated calcium influx in T cell proliferation, apoptosis and tolerance. Cell Calcium. 2011;50:261–9.
Zou J, Mercier C, Koussounadis A, Secombes C. Discovery of multiple beta-defensin like homologues in teleost fish. Mol Immunol. 2007;44:638–47.
Cuesta A, Meseguer J, Esteban MÁ. Molecular and functional characterization of the gilthead seabream β-defensin demonstrate its chemotactic and antimicrobial activity. Mol Immunol. 2011;48:1432–8.
Wang R, Feng J, Li C, Liu S, Zhang Y, et al. Four lysozymes (one c-type and three g-type) in catfish are drastically but differentially induced after bacterial infection. Fish Shellfish Immunol. 2013;35:136–45.
Purcell MK, Smith KD, Hood L, Winton JR, Roach JC. Conservation of toll-like receptor signaling pathways in teleost fish. Comp Biochem Physiol Part D Genomics Proteomics. 2006;1:77–88.
Valero Y, Morcillo P, Meseguer J, Buonocore F, Esteban MA, et al. Characterization of the IFN pathway in the teleost fish gonad against vertically transmitted viral nervous necrosis virus. J Gen Virol. 2015;96:2176–87.
Gu YF, Fang Y, Jin Y, Dong WR, Xiang LX, et al. Discovery of the DIGIRR gene from teleost fish: a novel Toll-IL-1 receptor family member serving as a negative regulator of IL-1 signaling. J Immunol. 2011;187:2514–30.
Grim JM, Hyndman KA, Kriska T, Girotti AW, Crockett EL. Relationship between oxidizable fatty acid content and level of antioxidant glutathione peroxidases in marine fish. J Exp Biol. 2011;214:3751–9.
Sneddon AA, Wu HC, Farquharson A, Grant I, Arthur JR, et al. Regulation of selenoprotein GPx4 expression and activity in human endothelial cells by fatty acids, cytokines and antioxidants. Atherosclerosis. 2003;171:57–65.
Wang G, Yang E, Smith KJ, Zeng Y, Ji G, et al. Gene expression responses of threespine stickleback to salinity: implications for salt-sensitive hypertension. Front Genet. 2014;5:312.
Andersen CL, Jensen JL, Ørntoft TF. Normalization of real-time quantitative reverse transcription-PCR data: a model-based variance estimation approach to identify genes suited for normalization, applied to bladder and colon cancer data sets. Cancer Res. 2004;64:5245–50.
Friberg IM, Lowe A, Ralli C, Bradley JE, Jackson JA. Temporal anomalies in immunological gene expression in a time series of wild mice: signature of an epidemic? PLoS One. 2011;6:e20070.
Steiger JH. Testing pattern hypotheses on correlation matrices: alternative statistics and some empirical results. Multivar Behav Res. 1980;15:335–52.
Jennrich RI. An asymptotic X 2 test for equality of 2 correlation matrices. J Am Stat Assoc. 1970;65:904–12.
Mantel N. The detection of disease clustering and a generalized regression approach. Cancer Res. 1967;27:209–20.
Manly BFG. Multivariate statistical methods: a primer. London: Chapman & Hall; 1986.
Margolin AA, Wang K, Lim WK, Kustagi M, Nemenman I, et al. Reverse engineering cellular networks. Nat Protocols. 2006;1:662–71.
Williams TD, Turan N, Diab AM, Wu H, Mackenzie C, et al. Towards a system level understanding of non-model organisms sampled from the environment: a network biology approach. PLoS Comput Biol. 2011;7, e1002126.
Lenz TL, Eizaguirre C, Rotter B, Kalbe M, Milinski M. Exploring local immunological adaptation of two stickleback ecotypes by experimental infection and transcriptome-wide digital gene expression analysis. Mol Ecol. 2013;22:774–86.
Huang Y, Chain FJJ, Panchal M, Eizaguirre C, Kalbe M, et al. Transcriptome profiling of immune tissues reveals habitat‐specific gene expression between lake and river sticklebacks. Mol Ecol. 2016;25:943–58.
Haase D, Rieger JK, Witten A, Stoll M, Bornberg-Bauer E, et al. Immunity comes first: the effect of parasite genotypes on adaptive immunity and immunization in three-spined sticklebacks. Dev Comp Immunol. 2015;54:137–44.
Haase D, Rieger JK, Witten A, Stoll M, Bornberg-Bauer E, et al. Specific gene expression responses to parasite genotypes reveal redundancy of innate immunity in vertebrates. PLoS One. 2014;9:e108001.
Dunne DW, Cooke A. A worm’s eye view of the immune system: consequences for evolution of human autoimmune disease. Nat Rev Immunol. 2005;5:420–6.
Mittal M, Siddiqui MR, Tran K, Reddy SP, Malik AB. Reactive oxygen species in inflammation and tissue injury. Antioxid Redox Signal. 2014;20:1126–67.
Cram DL, Blount JD, York JE, Young AJ. Immune response in a wild bird is predicted by oxidative status, but does not cause oxidative stress. PLoS One. 2015;10:e0122421.
Pacitti D, Lawan MM, Feldmann J, Sweetman J, Wang T, et al. Impact of selenium supplementation on fish antiviral responses: a whole transcriptomic analysis in rainbow trout (Oncorhynchus mykiss) fed supranutritional levels of Sel-Plex®. BMC Genomics. 2016;17:1–26.
Biller-Takahashi JD, Takahashi LS, Mingatto FE, Urbinati EC. The immune system is limited by oxidative stress: dietary selenium promotes optimal antioxidative status and greatest immune defense in pacu piaractus mesopotamicus. Fish Shellfish Immunol. 2015;47:360–7.
Machado LR, Ottolini B. An evolutionary history of defensins: a role for copy number variation in maximizing host innate and adaptive immune responses. Front Immunol. 2015;6:115.
Friberg IM, Little S, Ralli C, Lowe A, Hall A, et al. Macroparasites at peripheral sites of infection are major and dynamic modifiers of systemic antimicrobial pattern recognition responses. Mol Ecol. 2013;22:2810–26.
We are very grateful to Rob Darby, Rory Geohagen and Gareth Owen (IBERS, Aberystwyth University) for assistance.
This work was funded by research grant RPG-301 from the Leverhulme Trust.
Availability of data and materials
Datasets will be available in the European Nucleotide Archive under primary accession number PRJEB13319. Other supporting data are available as additional files.
JAJ conceived of and managed the study. MB, PH and JAJ carried out the fieldwork and mesocosm work. MB and PH carried out the gene expression assays. PH and JAJ analysed the data and JAJ wrote the manuscript with input from PH and IMF. All authors read and approved the final manuscript.
The authors declare that they have no competing interests.
Consent for publication
Ethics approval and consent to participate
Studies were approved by the animal welfare committee of the Institute of Biological, Environmental and Rural Sciences (IBERS), Aberystwyth University.
Biometric and environmental variation in the study populations and sites. (PDF 42 kb)
Significant effects of habitat and habitat × time on the expression of immune-associated genes. (PDF 322 kb)
Differing gene expression variance in wild and mesocosm fishes. (PDF 136 kb)