Abstract
Individuals vary in their immune response and, as a result, some are more susceptible to infectious disease than others. Little is known about the nature of this individual variation in natural populations, or which components of immune pathways are most responsible, but defining this underlying landscape of variation is an essential first step to understanding the drivers of this variation and, ultimately, predicting the outcome of infection. We describe transcriptome-wide variation in response to a standardised immune challenge in wild field voles. We find that genes (hereafter 'markers') can be categorised into a limited number of types. For the majority of markers, the response of an individual is dependent on its baseline expression level, with significant enrichment in this category for conventional immune pathways. Another, moderately sized, category contains markers for which the responses of different individuals are also variable but independent of their baseline expression levels. This category lacks any enrichment for conventional immune pathways. We further identify markers which display particularly high individual variability in response, and could be used as markers of immune response in larger studies. Our work shows how a standardised challenge performed on a natural population can reveal the patterns of natural variation in immune response.
Similar content being viewed by others
Introduction
Individuals vary in their immune response. Within a population, some individuals may fail to make protective immune responses following either natural infection or vaccination and so are especially vulnerable to infectious disease1,2,3,4. Defining the patterns of such variability will enhance our ability to manage the health of individuals – especially those that are most susceptible to infectious disease in human, livestock or wildlife populations.
Studies in laboratory mice are the cornerstone of immunology and have provided a detailed understanding of the molecular and cellular pathways by which immune responses are effected. This impressive mechanistic understanding, however, has only been achieved by minimising genetic and environmental variation within a laboratory setting. Where laboratory studies have examined the effects of variability – in genetics, microbiota or diet – both qualitative and quantitative differences in immune responses have been observed, with consequent effects on infection5,6,7. Nevertheless, natural variability cannot be fully reproduced in the laboratory, which has led to a recent effort to characterise the immune response in wild populations of mice or other rodents. Recent work in mice from agricultural and other anthropogenic settings is consistent with the expectation that exposure to complex environments greatly alters immune function8. New populations of memory T cells, present only in non-laboratory mice, have also been identified9.
One commonly used measure of an immune response is to assess the amount of one or more markers (e.g. transcripts or proteins) produced by a population of cells following stimulation by an immune agonist. From this ex vivo assay, one can gain insight into the types of immune response that could be made to a pathogen in vivo. Such responses will of course depend on the cell type, the time point and the immune agonist used. Nevertheless, for any molecular marker, individuals, in natural populations especially, could exhibit different marker abundances prior to and/or following stimulation (hereafter ‘baseline’ and ‘stimulated’ abundances respectively), leading to differences in their response to stimulation, and likely differences in their ability to control infection. An individual’s response to stimulation is here defined as the difference between its baseline and stimulated abundance of a particular marker. Furthermore, we argue that the most useful (and interesting) markers, in terms of understanding why individuals vary in their ability to mount a successful immune response, are likely to be those for which response to stimulation is most variable among individuals.
We use a wild population of field voles (Microtus agrestis) to examine naturally occurring patterns of individual variation in response to stimulation, across the transcriptome, as an essential first step towards understanding the processes driving these patterns. Given that we study a mixed population of cells (splenocytes), variation in the responses we measure could result from both the relative proportion and activation state of different cell types in different individual animals. Each of these things could, potentially, be driven independently by different causal drivers in individual animals and, furthermore, both could affect the strength and phenotypic direction of any response. Nonetheless, such composite responses are still likely to map onto functional responses and to represent a useful marker of individual variability.
We describe three main categories of response to stimulation: (i) uncorrelated response, (ii) constant response and (iii) baseline-dependent response (depicted in Fig. 1). We also identify markers, across these three categories, which show particularly high individual variability in response. We suggest that such categorisation is useful in organising natural variation in response, since little is known about which components of immune pathways are responsible for natural variability in immune response, or about the nature of such variability. Indeed, this categorisation is not limited to the components of conventional immune pathways. The ability of an immune response to effect protection against infection, for example, will be supported by a variety of non-immune functions, that will also be activated following stimulation by an agonist, and vary to a greater or lesser extent among individuals within a natural population. By identifying the components (whether conventionally immunological or not) that are likely to be responsible for natural variability in immune response, and by describing the nature of their variability, we are laying the groundwork for exploring the processes, whether genetic or environmental, which drive individual variation in immune response.
Methods
Field methods
Sixty-two field voles were collected between July and October 2015 to assay expression by RNASeq. Voles came from four sites in Kielder Forest, Northumberland (55° 13′ N, 2° 33′ W). Each site contained a trapping grid of regularly spaced traps (at approximately 5 m intervals) and was also used for other components of a larger field study (for more details see10).
Ethics statement
All animal procedures carried out as part of this larger field study were performed with approval from the University of Liverpool Welfare Committee and the UK Home Office (PPL 70/8210), and in accordance with the Animals (Scientific Procedures) Act 1986.
Laboratory methods
Following collection, voles were transported to the laboratory where they were weighed, aged (as either mature or immature) and sexed. Maturity (or immaturity) was determined ultimately following dissection. Males with scrotal testes and females with perforate vaginas were classified as mature. Non-scrotal males and non-perforate females can either be immature, or mature but sexually inactive. In these cases a consensus classification based on the internal reproductive organs, coat colour and size was used. Sampled voles varied in body weight (Median = 22.5 g; 25% centile = 18.8 g; 75% centile = 30.2 g) and included both mature (n = 43) and immature (n = 19) individuals. As our UK Home Office licence did not allow us to collect overtly pregnant females, there are fewer females (n = 18) than males (n = 44) in our sample. Voles were killed by a rising concentration of CO2 followed by exsanguination.
Splenocyte cultures
The whole spleen was removed from each vole, and a splenocyte culture was established from this. The splenocyte culture from each vole was split into two populations, one of which was stimulated with anti-CD3 antibodies (Hamster Anti-Mouse CD3e, Clone 500A2 from BD Pharmingen) and anti-CD28 antibodies (Hamster Anti-Mouse CD28, Clone 37.51 from BD Tombo Biosciences) at concentrations of 2 μg/ml and of 1 μg/ml respectively for 24 hr, and the other was left as an unstimulated control to act as a reference level. We refer to this reference level as the 'baseline', and control samples as 'baseline' samples. However, it is important to note that this level will vary for an individual, not only on a day to day basis, but throughout its life. Cell culture conditions were otherwise equivalent to those used in Jackson et al.11. Costimulation with anti-CD3 and anti-CD28 antibodies was used to selectively promote the proliferation of T cells12. We assumed that this would reflect the potential response of T-cell populations, a central process contributing to complex immune responses against invading pathogens, in vivo. We did not count the numbers of cells of different types (e.g. T cells, B cells or monocytes) in each spleen, but expect this to vary between spleens and therefore between splenocyte cultures.
RNASeq preparation and mapping
RNA was extracted from splenocyte cultures using Invitrogen PureLink kits. The amount of total RNA, and A260/A280 and A260/A230 nm ratios were assessed using the NanoDrop 1000 (Thermo Fisher Scientific). RNA integrity was assessed using the Agilent 2100 Bioanalyzer (Agilent Technologies). In each sample the RNA integrity number was measured (RINs > 7). For each sample, 100 ng of total RNA were used to prepare cDNA libraries using Illumina RiboZero kits and to construct libraries with NEBNext Ultra directional RNA library prep kit according to the manufacturer's protocols. Samples were sequenced to produce 2 × 75 bp paired-end reads on an Illumina HiSeq. 4000 platform. Adaptor sequences were removed with CUTADAPT version 1.2. and further trimmed with SICKLE version 1.200 (minimum window quality score of 20). This resulted in a mean library size of 18 million (range = 5–50 million) paired-end reads.
High-quality reads were mapped against a draft genome for M. agrestis (GenBank Accession no. LIQJ00000000) using TOPHAT version 2.1.013, and a set of predicted gene models was generated using BRAKER214. Mapped reads were counted using FEATURECOUNTS15. Further analysis was performed on counts of mapped reads for each gene in R version 3.4.216. These count data were initially filtered to remove unexpressed genes (those genes with fewer than three counts per million across all samples; n = 13). Following filtering, library sizes were recalculated and data were normalised to represent counts per million (cpm). No correction for gene length was necessary as all analyses were based on comparisons across (rather than within) samples. Transcript abundance for a particular gene here represents a single, functional measure of its activity across some, mixed, cell population.
In order to validate these RNASeq data, two-step reverse transcription quantitative PCR (Q-PCR) was performed on the same (unstimulated) samples and the expression of four example genes was measured: Arg1, Foxp3, Cd8a and Il1b. Foxp3 and Arg1 primer sequences were designed de novo and validated in house. Arg1 primer sequences were as follows GCAATTGGAAGCATCTCTGG (forward primer sequence) and TGATATCGGTGTGAGGATCC (reverse primer sequence) and Foxp3 primer sequences were CAGCAGCTGGAGCTGGAAAA (forward primer sequence) and GGACCAGGCTGGGAGAACAC (reverse primer sequence). Cd8a and Il1b primers were designed de novo and supplied by Primer Design (Chandler’s Ford, UK; see10 for detailed Q-PCR methodology). Expression levels estimated by RNASeq and Q-PCR were found to be positively correlated for three out of the four example genes (Arg1: rho = 0.61, p < 0.0001; Cd8a: rho = 0.52, p < 0.0001; Il1b: rho = 0.63, p < 0.0001; Supplementary Fig. S1a–c).
In order to maximise the power of our analysis to identify biologically relevant patterns, we focussed on those genes which were expressed at an ‘informative’ level in the spleen prior to and/or following stimulation. Genes expressed at a mean level greater than 200 cpm were considered ‘informative’ (n = 1350 or 6%). Of these, 1159 were annotated (see Supplemetary Table S1 for a full list of these genes). A relationship between gene-wise mean expression level and variance is typically found in RNASeq data, with low mean transcript abundance being strongly associated with high variability17 and therefore low repeatability (e.g. Foxp3; Supplementary Fig. S1d). We chose a threshold of 200 cpm here, because this is the region in which the mean-variance relationship asymptotes for our data i.e. a gene’s variance becomes independent of its mean expression level (Supplementary Fig. S2). This threshold is of a similar magnitude to that reported for a similar dataset, composed of unrelated human individuals, with high levels of biological variation17. As weakly expressed genes were removed (minimising heteroscedasticity), log-transformation of data was unnecessary. Both baseline and stimulated samples varied in their expression of informative genes (Supplementary Fig. S3).
Statistical methods
Genes for which a response to stimulation was observed across individuals were identified, and, as elaborated in the Results, categorised on the basis of (i) the dependence of an individual’s response on its baseline abundance, and (ii) the degree of variability in response across individuals.
Baseline-dependence of response
The dependence of an individual’s response on its baseline abundance was quantified by testing the relationship between that individual’s baseline abundance (cpmbase) and its stimulated abundance (cpmstim) using a linear regression (also known as an ordinary least squares or OLS regression), taking the form
as well as a quadratic regression, taking the form
Thus evaluating if stimulated abundance showed a curvilinear relationship with differences in baseline abundance. For approximately one third of genes (n = 466), the residuals from both of these regressions deviated significantly from the assumptions of normality and/or homoscedasticity, and a non-parametric Kendall–Theil linear regression was fitted instead.
Regression fits varied from gene to gene (R2 ranging from <0.001 to 0.85). As well as evaluating the significance of the quadratic term, the slopes of linear regression fits were tested for whether they were (i) significantly different from zero, and (ii) significantly different from one (we report whether slopes were greater or less than one in the text). Intercepts of linear regression fits were also tested for whether they were significantly different from zero. All resulting p-values, used to categorise genes, were adjusted for multiple testing using the Benjamini-Hochberg method18.
OLS regression assumes that the variable on the x-axis is measured without measurement error. Although this is not true for baseline abundance, by excluding those genes with low mean abundance, and focussing instead on high abundance genes, we were able to minimise this measurement error. As a result, OLS regression was still deemed appropriate, despite a small bias, especially as our goal was not precise slope estimation, but rather, assigning slopes to broad categories, and statistically testing whether slopes were different from reference values19.
Individual variability in response
Individual variability in response was quantified by comparing the coefficient of variation (CV) for baseline abundances across individuals (CVbase), and the CV for stimulated abundances across individuals (CVstim). Given that an individual’s response is defined as the difference between its baseline and stimulated abundance, a significant difference in the individual variability present at these two levels, either
or
is indicative of high individual variability in response. As we restricted our analysis to informative genes only, excluding those genes with low mean abundance, it was not necessary to account for the gene-wise mean-variance relationship (see above). Asymptotic tests for the equality of CVs were run using the R package cvequality20. All resulting p-values, used to further categorise genes, were adjusted for multiple testing using the Benjamini-Hochberg method18.
Functional annotation
Functional enrichment analysis was performed on gene sets, to uncover their biological functions, using The Database for Annotation, Visualization and Integrated Discovery (DAVID)21,22 version 6.8. Only genes which were annotated in the draft genome for M. agrestis could be included in gene sets and analysed in this way. The population background against which gene sets were analysed was composed of all annotated genes in the draft genome (n = 13426 M. agrestis genes mapped to n = 14161 Mus musculus genes). Benjamini-Hochberg corrected p-values and gene counts are reported alongside ontology terms of interest, including Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways23,24,25, to indicate their level of enrichment (significance cut-off = 0.1). A small number of genes (e.g. predicted genes) were unmapped in DAVID and were excluded from the enrichment analysis.
Age-specific analysis
In order to begin to investigate the relative importance of genetic variation versus prior stimulation for shaping patterns of variation in response, the same analysis was performed separately on immature and mature voles. As we had more samples from mature voles (n = 43) than immature voles (n = 19), we randomly sampled the mature population (with replacement) 1000 times and averaged across these samples. The number (immatures) or mean number (matures) of genes in a response category for each of these age classes is presented in the text.
Results
Stimulation with an immune agonist causes a widespread response
1150 genes (5% of all genes in the field vole genome and 85% of informative genes, those genes which were more strongly expressed in the spleen) fell into one or more of the response categories set out in Fig. 1. As expected, given that anti-CD3 and anti-CD28 antibodies are known to stimulate T-cell proliferation12, they were enriched with genes (hereafter ‘markers’) associated with the T-cell receptor (TCR) signalling pathway (n = 26; p < 0.001) and other T cell-related terms: positive regulation of T-cell proliferation (n = 12; p = 0.08), TCR complex (n = 7; p < 0.01), positive thymic T-cell selection (n = 7; p = 0.01), negative thymic T-cell selection (n = 6; p = 0.05) and alpha-beta TCR complex (n = 5; p < 0.01). For the majority of these markers, a significant positive linear relationship was found between baseline and stimulated abundance (n = 844). Only a single marker, Fam193b, demonstrated a significant negative linear relationship between baseline and stimulated abundance.
There are three main categories of response to stimulation
Three main categories of response to stimulation were identified based on the dependence of an individual’s response on its baseline abundance. Each of these categories demonstrates a unique pattern (Fig. 1):
Uncorrelated response
Markers for which individuals taken from the wild differ in their baseline abundance, but the responses of different individuals are variable and independent of their baseline, such that the slope of the relationship between baseline and stimulated abundance is not significantly different from zero.
Constant response
Markers for which individuals taken from the wild also differ in their baseline abundance, but the responses of different individuals are (on average) constant and independent of their baseline, such that the slope of the relationship between baseline and stimulated abundance is not significantly different from one and the intercept (indicating the level of response) is significantly greater than zero.
Baseline-dependent response
Markers for which individuals taken from the wild again differ in their baseline abundance, but the responses of different individuals vary as a function of their baseline, either as a linear function of their baseline (slope significantly different from one), or as a quadratic function of their baseline, where stimulated abundance either increases exponentially as a function of baseline abundance or becomes saturated at some upper limit.
We also identified markers, across these three categories, which displayed high individual variability in response. These highly variable markers can be further divided into two categories (Fig. 1):
Convergent response
Markers for which individual variability in baseline abundance is significantly greater than individual variability in stimulated abundance.
Divergent response
Markers for which individual variability in stimulated abundance is significantly greater than individual variability in baseline abundance.
The baseline-dependent response category is most common and is significantly enriched in components of conventional immune pathways
The baseline-dependent response category was the most common (Table 1), and included a majority of markers for which stimulated abundance was a linear function of baseline abundance (n = 539), and a remainder for which it was a quadratic function (n = 160). The majority of quadratic response markers showed evidence for saturation (n = 138), indicating some upper limit on stimulated abundance. Functional enrichment analysis revealed significant enrichment for a number of immune ontology terms in the linear response category, including T helper cell surface molecules (n = 6; p = 0.08), leukocyte transendothelial migration (n = 13; p = 0.03) and HTLV-1 infection (n = 20; p = 0.05). Only the TCR signaling pathway was enriched in the quadratic response category (n = 7; p = 0.03; Fig. 2).
The uncorrelated response category is least common and lacks enrichment in components of conventional immune pathways
A number of markers showed no evidence for a relationship between baseline and stimulated abundance (n = 47; Table 1). For the majority of these, mean abundance was significantly greater for stimulated than for baseline samples (n = 39; p < 0.05), suggesting that these markers were (on average) responding to stimulation, but to an individually variable degree, independent of baseline abundance. These markers lacked any enrichment for immune-related terms (Fig. 2).
A number of markers, including Zap70 , show particularly high individual variability in response
For a number of markers, variability in baseline and stimulated abundance was significantly different, leading to high individual variability in response (n = 244). The vast majority of these markers showed a divergent (n = 237), rather than a convergent (n = 7) response. Within the (stimulated) TCR signalling pathway, the highest level of variability in individual response, and the highest level of divergence, was demonstrated by Zap70 (Fig. 3). All convergent markers fell into one of the three main categories of response. However, over a third of divergent markers (n = 98), did not fall into any of these categories (Table 1), appearing instead as markers which (on average) did not respond to stimulation. There was also no significant difference in the mean abundance of these markers in stimulated and baseline samples (p > 0.05).
Immature voles show more individual variability in response than mature voles
An age-specific analysis showed that high individual variability in response to stimulation (whether divergent or convergent) was more common among immature voles (numbers of markers; divergent = 108; convergent = 6) than mature voles (mean numbers of markers and empirical 95% intervals; divergent = 50 [0–338.2]; convergent = 0.11 [0–1]).
Response to stimulation is not limited to components of conventional immune pathways
Functional enrichment analysis revealed that non-immune related ontology terms were also enriched in the response categories, including: toxin transport (n = 6; p = 0.02 in divergent response category; n = 6; p = 0.05 in constant response category) and dioxygenase (n = 7; p = 0.002 in baseline-dependent response category). The top convergent response marker, Pdk1, is also a component of the insulin signalling pathway (Fig. 2).
Discussion
The need to better understand variation in immune response in natural populations is now widely accepted26,27,28,29. Our understanding of immune responses in laboratory settings comes from animals that vary little either genetically or in prior experience. By contrast, animals in natural populations vary (perhaps extensively) in both of these. In this study, we describe the patterns of natural variation in response to a standardised immune challenge in a wild population of rodents. We identify three main categories of response: uncorrelated response, constant response and baseline-dependent response. We also identify markers, across these categories, which show particularly high individual variability in response.
The baseline-dependent response category is the largest. Markers in this category show a relationship between baseline and stimulated abundance across individuals, and their response to stimulation is (to a lesser or a greater extent) dependent on their baseline abundance. In some cases, individuals already expressing the greatest abundance of a marker in their natural setting went on to exhibit the greatest response to stimulation by an agonist. In others, the opposite was true, and these individuals exhibited the smallest response to stimulation. Similarly, previous work on humans has identified baseline (transcriptional) predictors of influenza vaccination response30,31. These differences in baseline abundance could be driven by either genetic variation or individual differences in past experience. In humans, genetic determinants of baseline immune cell population frequencies have been identified32. Even though the stimulation we describe here was not antigen specific, previous challenge by a pathogen might also lead to changes in the baseline T-cell population within an individual’s spleen, affecting its response to any subsequent challenge. In fact, we find that voles infected with Babesia microti (a blood parasite, common in our population33) have larger spleens than uninfected voles34. This prior experience may prime an individual, enabling a greater response to subsequent challenge (e.g. slope greater than one; Fig. 1). However, individuals may also have an upper limit on the number of immune cells they have available35,36. An individual that is already mounting an immune response to a pathogen, and has a large number of activated T cells, could therefore respond less to a similar challenge than an ‘immunologically naïve’ individual (slope less than one; Fig. 1). Membership of the baseline-dependent response category recapitulates the known biology of the immune response (being highly enriched for immune ontology terms). In doing so, it validates the approach we use here, as a way of identifying markers of immune significance.
For both the uncorrelated response category and the constant response categories, individuals varied in their natural abundance of a marker but their response was unrelated to this. For the uncorrelated response category, the majority of markers occurred at a significantly higher mean abundance in stimulated samples than in baseline samples suggesting that they did nevertheless respond to stimulation. This category, which contains a moderate number of markers, also lacks any enrichment for immune-related ontology terms. This suggests that markers in this category are not conventional immune markers but could be of immune significance. We warn against omitting such markers from studies of immune response in the laboratory. They could play an important part in our understanding of the immune response, indicating for example, genetic variation in response among individuals, which is independent of baseline abundance.
The constant response category was the second largest and was moderately enriched for immune ontology terms. Markers in this category showed (on average) a constant response to stimulation across individuals. This pattern may be characteristic of markers that are of critical importance in mounting a successful immune response to a pathogen in vivo, and are therefore under tight regulation. For example, one of the top markers in the constant response category, Fyn, plays a central role in initiating T-cell activation and expansion37,38. Fyn-deficient mice have a severely compromised response to staphylococcal enterotoxin B, demonstrating the importance of this gene for mounting a successful response39. However, this is not the case for all markers in this category. Some markers showed individual variation around the average (constant) response and, for almost a third, responses of different individuals were highly variable (see below).
Cutting across this categorisation, a large number of markers displayed a pattern in which variation between individuals was particularly strong. We describe two types of such markers, both of which could be used in future studies as indicators of natural variability in immune response. Markers in the less common, convergent, response category showed much greater variation naturally than following stimulation. This pattern may be characteristic of markers showing variable levels of prior activation and some maximum or optimum abundance, resulting in a stabilisation of marker levels across the population following stimulation. We found that convergent patterns were more common among immature voles.
Assuming that individuals converge immunologically over time as they become colonised by the most common pathogens in the environment, we would expect young voles to be more variable in their levels of prior activation, as a result of their shorter and more variable exposure histories40. This could make them more likely to converge. Equally though, this could suggest that they are more constrained in the energy they have available (as a result of the competing energetic demands of growth and development) or the number of immune cells they have available (as a result of a developing immune system). Either resource constraint could result in a maximum abundance, also making them more inclined to converge. Due to the costly nature of the immune response, individuals often trade-off their investment in different arms of the immune system41,42. Different types of immune response are therefore likely to be associated with different optimum abundances (or regions) and an individual already mounting an immune response, but to a different type of challenge (associated with different cell types), may respond by down-regulating expression.
Divergent markers, which were more common, showed much greater variation following stimulation than there was naturally. This pattern may be characteristic of (but not limited to) markers showing genetic variation in response to the agonist, independent of baseline abundance e.g. subsets of animals that appear similar but respond more strongly to stimulation than others. Our own recent work, where we found an association between polymorphism in a single gene and a marker of a more tolerant immune response43, is an example of such genetic variation in response. Further supporting this hypothesis, here, we found more divergent markers among immature voles than mature voles. Assuming that individuals diverge, rather than converge, immunologically over time as a result of the cumulative influence of environmental exposure4, we would expect genetic effects to be more easily detectable in younger voles, with shorter exposure histories. Equally, though, divergent patterns could be the result of differences in early life experiences. One would also expect these to be more easily detectable in immature voles.
The divergent category (predominantly) included markers for which individuals made (on average) the same response to stimulation, and markers that did not respond (on average) to stimulation. Standard differential expression analysis would miss the individual variation present in the former group, and would fail to pick up the latter group of markers altogether. Both warn against looking at average (population-level) response, and point instead, to the value of looking at individual-level differences in response. This is particularly important because divergent markers may act as critical regulators of pathways. For example, Zap70, which demonstrates particularly high levels of individual variability in response and is centrally located in the TCR signalling pathway, interacts with many other markers (Fig. 3). We suggest that Zap70 expression could be used as a marker of response in larger studies. Indeed, it is already linked to major seasonal immune variation in wild fish44 and is being used as a prognostic marker for B-cell chronic lymphocytic leukemia in humans, with potential implications for determining a patient’s treatment path (recently reviewed in45). This example supports our assertion, in this study, that markers for which response is most variable among individuals are likely to be most useful for understanding why individuals vary in susceptibility. Other potential prognostic (or diagnostic) markers, which may have been missed using standard differential expression analysis, may be present in this category and warrant further investigation. Individual markers do not work in isolation though, and the immune system is a ‘multi-faceted defence system’46. For example, one study showed that 24 immunological parameters explained 60% of the variation in resistance to streptococcus group B in laboratory mice when combined in a multivariate analysis, providing insight which would have been missed by analysing each immunological parameter separately47. Similarly, groups of markers of prognostic and/or diagnostic value may be identified within this category by applying multivariate statistical methods.
The response categories we describe here are based on splenocytes (i) stimulated with anti-CD3 and anti-CD28 antibodies, (ii) sampled at 24 hours and (iii) assayed by RNASeq. Firstly, we chose to use a broad spectrum stimulus, as opposed to a pathogen-specific challenge, as we were looking for broadly applicable response categories. One caveat on our results is that our choice of agonist specifically activates and expands T cells. Our results may therefore exclude or underrepresent host responses that are independent of, or less strongly driven by, T cells. Furthermore, unlike an immune agonist, a pathogen is a dynamic, living organism which can multiply, evolve and interact with the host’s immune response by e.g. evasion or manipulation. In order to fully understand why individuals vary in their susceptibility to infection then, we must consider the pathogen as well as the host’s response46. We hope that future work will verify the categories we set out here using other immune agonists or, ideally, a range of pathogen-specific challenges.
Secondly, we sampled our spleen cells at 24 hours post stimulation, but the choice of time point may affect the relative frequency of the response categories reported here. Markers are known to follow different response trajectories, with some immediately responding and reaching peak activation, and others taking longer to reach this point48. Sampling at a later time point, then, when the ‘slower’ markers have reached peak activation, may lead to more convergence than reported here. In order to fully account for this temporal variation, multiple time points need to be averaged across. We argue that both time-specific and averaged responses are of functional significance, but hope others will extend our work.
Thirdly, we chose to assay expression by RNASeq, again, in order to give a broad view of the response to stimulation. Previous work has shown that transcript levels generally correlate with protein levels across genes49. However, more work is needed to confirm these response categories at this functional level50. In future, quantitative PCR (Q-PCR) or protein-level data could be used in order to include weakly expressed markers, which were excluded here as a result of the heteroscedasticity inherent in RNASeq data. Single-cell RNASeq could also be used to quantify differences in individual response resulting solely from differences in cell-specific activity.
Markers that responded to stimulation were not limited to immune pathways as conventionally defined. They included markers involved in energy generation, including the insulin signalling pathway. This is in line with previous studies, which suggest that insulin plays a key role in coordinating an organism’s response to infection, influencing, in particular, the allocation of resources51,52. The marker in question, Pdk1, was among the top convergent markers. This could be representative of the high levels of variability in the (baseline) nutritional status of individuals in a natural population, coupled with an upper limit on the processes involved in glucose metabolism.
The response categories we presented here, therefore, highlight markers not traditionally associated with immune functions, and offer a promising avenue for identifying potential prognostic (or diagnostic) factors for disease, like Zap70. They also point to both genetics and prior experience as important drivers of immune variability, providing the next challenge for understanding why individuals vary in their response to infection. The field population we study has been the subject of extensive previous study on population ecology and pathogen dynamics33,34,53,54. Our future work will, therefore, place the patterns of natural variation described here into this context, in order to decompose them into those driven by genetics and prior experience.
Data availability
RNASeq data are available from the European Nucleotide Archive (study accession number PRJEB23617).
References
Sarre, C. et al. Comparative immune responses against Psoroptes ovis in two cattle breeds with different susceptibility to mange. Vet. Res. 46, (2015).
Mahler, M., Janke, C., Wagner, S. & Hedrich, H. J. Differential susceptibility of inbred mouse strains to Helicobacter pylori Infection. Scand. J. Gastroenterol. 37, 267–278 (2002).
Eugui, E. M. & Allison, A. C. Malaria infections in different strains of mice and their correlation with natural killer activity. Bull. World Health Organ. 57, 231–238 (1979).
Brodin, P. et al. Variation in the human immune system is largely driven by non-heritable influences. Cell 160, 37–47 (2015).
Hooper, L. V., Littman, D. R. & Macpherson, A. J. Interactions between the microbiota and the immune system. Science 336, 1268–1273 (2012).
Lochmiller, R., Vestey, M. & Nash, D. Gut associated lymphoid tissue in the cotton rat (Sigmodon hispidus) and its response to protein restriction. J. Wildl. Dis. 28, 1–9 (1992).
Sellers, R. S., Clifford, C. B., Treuting, P. M. & Brayton, C. Immunological variation between inbred laboratory mouse strains: points to consider in phenotyping genetically immunomodified mice. Vet. Pathol. 49, 32–43 (2012).
Abolins, S. R., Pocock, M. J., Hafalla, J. C., Riley, E. M. & Viney, M. E. Measures of immune function of wild mice, Mus musculus. Mol. Ecol. 20, 881–892 (2011).
Beura, L. K. et al. Normalizing the environment recapitulates adult human immune traits in laboratory mice. Nature 532, 512–6 (2016).
Wanelik, K. et al. IgE receptor polymorphism predicts divergent, sex-specific inflammatory modes and fitness costs in a wild rodent Preprint at https://www.biorxiv.org/content/10.1101/841825v1 (2019).
Jackson, J. A. et al. The analysis of immunological profiles in wild animals: a case study on immunodynamics in the field vole, Microtus agrestis. Mol. Ecol. 20, 893–909 (2011).
Frauwirth, K. A. & Thompson, C. B. Lymphocyte signal transduction Activation and inhibition of lymphocytes by costimulation. Cancer Res. 109, 295–299 (2002).
Trapnell, C., Pachter, L. & Salzberg, S. L. TopHat: discovering splice junctions with RNA-Seq. Bioinformatics 25, 1105–1111 (2009).
Hoff, K. J., Lange, S., Lomsadze, A., Borodovsky, M. & Stanke, M. BRAKER1: unsupervised RNA-Seq-based genome annotation with GeneMark-ET and AUGUSTUS. Bioinformatics 32, 767–769 (2015).
Liao, Y., Smyth, G. K. & Shi, W. featureCounts: an efficient general purpose program for assigning sequence reads to genomic features. Bioinformatics 30, 923–930 (2014).
R Core Team. R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. https://www.R-project.org/. (2018).
Law, C. W., Chen, Y., Shi, W. & Smyth, G. K. Voom: precision weights unlock linear model analysis tools for RNA-seq read counts. Genome Biol. 15, 1–17 (2014).
Benjamini, Y. & Hochberg, Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J. R. Stat. Soc. Ser. B 57, 289–300 (1995).
Kilmer, J. T. & Rodríguez, R. L. Ordinary least squares regression is indicated for studies of allometry. J. Evol. Biol. 30, 4–12 (2017).
Marwick, B. & Krishnamoorthy, K. cvequality: tests for the equality of coefficients of variation from multiple groups. R package version 0.2.0. https://CRAN.R-project.org/package=cvequality (2019).
Huang, D. W., Sherman, B. T. & Lempicki, R. A. Bioinformatics enrichment tools: paths toward the comprehensive functional analysis of large gene lists. Nucleic Acids Res. 37, 1–13 (2009).
Huang, D. W., Sherman, B. T. & Lempicki, R. A. Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat. Protoc. 4, 44–57 (2009).
Kanehisa, M. & Goto, S. KEGG; Kyoto Encyclopedia of Genes and Genomes. Nucleic Acids Res. 28, 27–30 (2000).
Kanehisa, M., Sato, Y., Kawashima, M., Furumichi, M. & Tanabe, M. KEGG as a reference resource for gene and protein annotation. Nucleic Acids Res. 44, D457–D462 (2016).
Kanehisa, M., Furumichi, M., Tanabe, M., Sato, Y. & Morishima, K. KEGG: new perspectives on genomes, pathways, diseases and drugs. Nucleic Acids Res. 45, D353–D361 (2017).
Jackson, J. A. et al. Immunomodulatory parasites and toll-like receptor-mediated tumour necrosis factor alpha responsiveness in wild mammals. BMC Biol. 7, 16 (2009).
Pedersen, A. B. & Babayan, S. A. Wild immunology. Mol. Ecol. 20, 872–880 (2011).
Turner, A. K. & Paterson, S. Wild rodents as a model to discover genes and pathways underlying natural variation in infectious disease susceptibility. Parasite Immunol. 35, 386–395 (2013).
Viney, M. & Riley, E. M. The immunology of wild rodents: current status and future prospects. Front. Immunol. 8, 1–9 (2017).
HIPC-CHI Signatures Project Team & HIPC-I Consortium. Multicohort analysis reveals baseline transcriptional predictors of influenza vaccination responses. Sci. Immunol. 2, eaal4656 (2017).
Tsang, J., Schwartzberg, P., Kotliarov, Y. & Al., E. Global analyses of human immune variation reveal baseline predictors of post-vaccination respones. Cell 157, 499–513 (2014).
Orrù, V. et al. Genetic variants regulating immune cell levels in health and disease. Cell 155, 242–256 (2013).
Telfer, S. et al. Species interactions in a parasite community drive infection risk in a wildlife population. Science 330, 243–246 (2010).
Taylor, C. H. et al. Physiological, but not fitness, effects of two interacting haemoparasitic infections in a wild rodent. Int. J. Parasitol. 48, 463–471 (2018).
Lochmiller, R. L. & Deerenberg, C. Trade-offs in evolutionary immunology: just what is the cost of immunity? Oikos 88, 87–98 (2000).
Sheldon, B. C. & Verhulst, S. Ecological immunology: costly parasite defences and trade-offs in evolutionary ecology. Trends Ecol. Evol. 11, 317–321 (1996).
Palacios, E. H. & Weiss, A. Function of the Src-family kinases, Lck and Fyn, in T-cell development and activation. Oncogene 23, 7990–8000 (2004).
Salmond, R. J., Filby, A., Qureshi, I., Caserta, S. & Zamoyska, R. T-cell receptor proximal signaling via the Src-family kinases, Lck and Fyn, influences T-cell activation, differentiation, and tolerance. Immunol. Rev. 228, 9–22 (2009).
Appleby, M. W. et al. Defective T cell receptor signaling in mice lacking the thymic isoform of p59fyn. Cell 70, 751–763 (1992).
Sol, D., Jovani, R. & Torres, J. Parasite mediated mortality and host immune response explain age-related differences in blood parasitism in birds. Oecologia 135, 542–547 (2003).
Reiner, S. L. Development in motion: helper T cells at work. Cell 129, 33–36 (2007).
Abbas, A. K., Murphy, K. M. & Sher, A. Functional diversity of helper T lymphocytes. Nature 383, 787–793 (1996).
Wanelik, K. M. et al. A candidate tolerance gene identified in a natural population of field voles (Microtus agrestis). Mol. Ecol. 27, 1044–1052 (2018).
Brown, M. et al. Seasonal immunoregulation in a naturally-occurring vertebrate. BMC Genomics 17, (2016).
Liu, Y., Wang, Y., Yang, J., Bi, Y. & Wang, H. ZAP-70 in chronic lymphocytic leukemia: a meta-analysis. Clin. Chim. Acta 483, 82–88 (2018).
Sadd, B. M. & Schmid-Hempel, P. Principles of ecological immunology. Evol. Appl. 2, 113–121 (2009).
Keil, D., Luebke, R. W. & Pruett, S. B. Quantifying the relationship between multiple immunological parameters and host resistance: probing the limits of reductionism. J. Immunol. 167, 4543–4552 (2001).
Tate, A. T. & Graham, A. L. Dissecting the contributions of time and microbe density to variation in immune gene expression. Proc. R. Soc. B Biol. Sci. 284, (2017).
Tian, Q. et al. Integrated genomic and proteomic analyses of gene expression in mammalian cells. Mol. Cell. Proteomics 3, 960–969 (2004).
Viney, M. E., Riley, E. M. & Buchanan, K. L. Optimal immune responses: immunocompetence revisited. Trends Ecol. Evol. 20, 665–669 (2005).
Garcia, N. W. et al. Exogenous insulin enhances humoural immune responses in short-day, but not longday, Siberian hamsters (Phodopus sungorus). Proc. R. Soc. B Biol. Sci. 277, 2211–2218 (2010).
Lee, S. et al. Modulatory upregulation of an insulin peptide gene by different pathogens in C. elegans. Virulence 9, 648–658 (2018).
Burthe, S. et al. Cowpox virus infection in natural field vole Microtus agrestis populations: delayed density dependence and individual risk. J. Anim. Ecol. 75, 1416–1425 (2006).
Smith, A., Telfer, S., Burthe, S., Bennett, M. & Begon, M. Trypanosomes, fleas and field voles: ecological dynamics of a host-vector–parasite interaction. Parasitology 131, 355–365 (2005).
Acknowledgements
The authors wish to thank those involved in obtaining and processing samples from the field: Rebecca Turner, Lukasz Lukomski, Stephen Price, Sarah Gore, Ed Parker, William Foster, Ann Lowe and Anna Thomason. They also wish to thank the Forestry Commission for access to the study sites and the Centre for Genomic Research at the University of Liverpool for sequencing samples. This research was funded by the Natural Environment Research Council (NERC) award NE/L013452/1 to SP, MB, JAJ and JEB.
Author information
Authors and Affiliations
Contributions
M.B., J.E.B., J.A.J. and S.P. designed the study. E.A. undertook the stimulatory assays. K.M.W. analysed the data. All authors wrote the manuscript.
Corresponding author
Ethics declarations
Competing interests
The authors declare no competing interests.
Additional information
Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Supplementary information
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.
About this article
Cite this article
Wanelik, K.M., Begon, M., Arriero, E. et al. Transcriptome-wide analysis reveals different categories of response to a standardised immune challenge in a wild rodent. Sci Rep 10, 7444 (2020). https://doi.org/10.1038/s41598-020-64307-7
Received:
Accepted:
Published:
DOI: https://doi.org/10.1038/s41598-020-64307-7
This article is cited by
-
Managing host-parasite interactions in humans and wildlife in times of global change
Parasitology Research (2022)
Comments
By submitting a comment you agree to abide by our Terms and Community Guidelines. If you find something abusive or that does not comply with our terms or guidelines please flag it as inappropriate.