Living with chronic infection: Persistent immunomodulation during avirulent haemoparasitic infection in a wild rodent

Apicomplexans are a protozoan phylum of obligate parasites which may be highly virulent during acute infections, but may also persist as chronic infections which appear to have little fitness cost. Babesia microti is an apicomplexan haemoparasite that, in immunocompromised individuals, can cause severe, potentially fatal disease. However, in its natural host, wild field voles (Microtus agrestis), it exhibits chronic infections that have no detectable impact on survival or female fecundity. How is damage minimized, and what is the impact on the host's immune state and health? We examine the differences in immune state (here represented by expression of immune‐related genes in multiple tissues) associated with several common chronic infections in a population of wild field voles. While some infections show little impact on immune state, we find strong associations between immune state and B. microti. These include indications of clearance of infected erythrocytes (increased macrophage activity in the spleen) and activity likely associated with minimizing damage from the infection (anti‐inflammatory and antioxidant activity in the blood). By analysing gene expression from the same individuals at multiple time points, we show that the observed changes are a response to infection, rather than a risk factor. Our results point towards continual investment to minimize the damage caused by the infection. Thus, we shed light on how wild animals can tolerate some chronic infections, but emphasize that this tolerance does not come without a cost.

must invest in resistance and/or tolerance of the infection to mitigate the impacts of its continued presence (Råberg, 2014).
We can describe an immune response in terms of an animal's "immune state"-that is, the numbers and activity of immune cells and associated molecules within different tissues at any given time. An "immune response" is then the changes in immune state that are associated with a particular challenge (Abolins et al., 2018). Like the dynamics of the infection itself, this response can vary during the time course of an infection, with some short-term changes to immune state occurring rapidly during acute infection and others persisting for longer (e.g., Chen et al., 2000;Lazarevic et al., 2005). Surprisingly little is known about the dynamics of response to chronic infection in wild hosts, and in particular the long-term impact on immune state. Any such impacts can have great importance as they in turn can affect, for example, energy expenditure (Martin et al., 2003), reproductive investment (Festa-Bianchet, 1989), inflammation (Finney et al., 2009), oxidative damage (Costantini, 2022) and response to coinfections (Lo et al., 2018). Thus, changes in immune state following infection in a wild individual have the potential for substantial impacts on longterm health and fitness (Budischak et al., 2018;Nussey et al., 2014).
Much of our knowledge of immune response to infection comes from laboratory experiments that use inbred strains low in genetic diversity, under near constant environmental conditions. However, to understand the true impact on wild individuals, infections must be studied in their natural setting, where immune responses can differ markedly from those recorded in the laboratory (Abolins et al., 2017). Abolins et al. (2018) examined a wide range of immune markers and found a remarkably weak link between serological measures of infection and immune state. However, they acknowledged that use of serological infection data (reflecting both historical and recent infections) may well have obscured important relationships.
In more targeted studies, there have been various demonstrations of links between particular immune markers and chronic helminth infections in the wild Jackson et al., 2009;Sparks et al., 2018). However, these, like most ecoimmunology studies, involve observational data, making it difficult to separate cause from effect given the wide array of factors which may simultaneously be acting on and responding to an individual's immune state. Using an experimental approach, removal of the helminth Heligmosomoides polygyrus did not result in major changes to gene expression in wild wood mice Apodemus sylvaticus (Babayan et al., 2018), again leaving the causal relationship between infection and immune response uncertain in that instance. In addition to experiments, longitudinal studies are a useful but under-used means to provide insights into causality since it is possible to draw inferences from the sequence of events without having to manipulate the host's conditions. In general terms, it is still unclear to what extent chronic infections impact a wild animal's immune state.
The apicomplexan parasite Babesia microti is widespread among various rodent species (Bajer et al., 2014;Barnard et al., 1993;Bown et al., 2008). It invades erythrocytes, sometimes causing a malarialike fever and anaemia (Homer et al., 2000;Lykins et al., 1975). If the initial acute infection is not fatal, it can then persist as a chronic infection for weeks or months with few symptoms. For example, field voles, Microtus agrestis, are commonly found to be infected with B. microti and, at least during the chronic phase, the infection shows low virulence: it causes splenomegaly but has no known negative effects on survival or fecundity . The parasite is transmitted by ticks, which can also pass B. microti to humans, causing life-threatening disease for older and immunocompromised individuals (Krause et al., 2008;Vannier et al., 2004).
Here, we describe the relationship between immune state and infection in wild M. agrestis, with our analysis divided into cross-sectional and longitudinal components. In the cross-sectional component, we obtained detailed physiological and immunological measurements from culled individuals and examined their association with several common chronic infections, while accounting for other influences such as seasonal variation. In the longitudinal component, we took advantage of an extensive data set of marked individuals with high recapture rates to examine changes within individuals over time while remaining in situ in their natural population. With this powerful data set, we can infer likely causal directions by examining sequences of events within individuals. We find that chronic infection with B. microti is associated with significant changes in expression of a wide range of measured genes. Other infections, including Bartonella spp. and helminths, demographic factors and condition show remarkably little association with our immune markers. We further show that changes to immune state follow B. microti infection and are reversed upon clearance. This indicates that, despite having low virulence in this population , B. microti is an important hidden driver of immune state variation.

| MATERIAL S AND ME THODS
All animal procedures carried out as part of this 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. Voles were marked with unique passive integrated transponders (PIT tags; Avid) on first capture, inserted subcutaneously on the back. On the first capture of each session for a given individual, we recorded sex and reproductive status (ascended or descended testes for males; pregnant, lactating, perforate or imperforate vagina for females) from visual inspection. We measured body length from snout to vent, and body mass. A small blood sample (up to 50 μl) was taken from the tip of the tail and stored in RNAlater at −20°C for a maximum of 5 days before being transferred to −80°C for longer term storage. This blood sample was used in gene expression and microparasite assays (see below). Individuals were then released unless selected for the cross-sectional component of the study.

| Cross-sectional sampling
On the last day of a session, a subset of individuals (up to a maximum of 20 per session) were retained in individual cages and transported to the University of Nottingham. They were then housed individually for one or two nights with access to bird seed mix, chopped carrot and water ad libitum. After this, voles were culled by a rising concentration of CO 2 , with death confirmed by exsanguination. These individuals formed the cross-sectional component of the study, allowing us to take more invasive measurements and samples (further details below): blood and spleen samples for measurement of gene expression, eye lens mass, packed cell volume (PCV), and counts of macroparasites in the fur and GI tract. In addition, we recorded sex, reproductive status, snout-to-vent length and body mass as described above for the longitudinal component of the study.

| Eye lenses
Eye lens mass can be used as a proxy for age in rodents, because lenses grow relatively continuously and predictably over time (Augusteyn, 2014). Eyes were removed after culling and stored in formalin at room temperature. Later, we removed the lenses by dissection and dried them at 60°C for 48 h. We then measured mass of both lenses together to 0.1 mg accuracy with an electronic balance.

| Packed cell volume
PCV measures the concentration of erythrocytes in circulation, with lower values indicative of anaemia (Wallerstein Jr., 1987).
A small volume of blood (<50 μl) was drawn up in a heparinised capillary tube from the thoracic cavity immediately after culling.
The sample was spun for 3 min in a microhaematocrit centrifuge (Hawksley) to separate cells from plasma and the ratio of the two (PCV) was calculated.

| Gene expression in blood
Blood samples from cross-sectional animals were taken from the thoracic cavity immediately after culling. Serum was separated by centrifugation at 13,000 rpm for 5 min and the cell pellets were placed in RNAlater and stored below −70°C. These cell pellets, plus the tail bleed samples from the longitudinal animals, were analysed by two-step reverse transcription SYBR Green-based quantitative real-time PCR to quantify levels of expression of a panel of 28 genes.
For our measurement panel, we curated a set of immuneassociated genes, or genes with other physiological functions of interest that might impinge on immunity based upon variation in our own data and known comparative functions in the immunology literature. Genes within the set were selected on the basis of: (a) significant tissue-specific sensitivity of expression to environmental (including infection) or intrinsic host variables in our previous studies (Jackson et al., 2011(Jackson et al., , 2014; (b) significant tissue-specific sensitivity of expression to environmental or intrinsic host variables in a recent differential expression analysis of RNASeq data based upon a small subset of our present cross-sectional samples (n < 70) (not reported here, and not including any of the infections considered here as explanatory terms as the relevant infection data was not available at the time of selection); and (c) representation of a distinctive immunological pathway, or pathway that impinges upon immunity, based upon the immunological literature. Compared to an alternative transcriptome-wide approach, targeted curation of a panel of candidate genes in the above way, based on baseline knowledge and pilot study data, may mitigate type 2 errors due to multiplicity correction in the case of small effects. Primers (30 sets, including 2 endogenous control genes) were designed de novo in-house and validated (to confirm specific amplification and 100 ± 10% PCR efficiency under our assay conditions). The endogenous control genes (ywhaz and actb) were selected as a stable pairing from our previous stability analysis of candidate control genes in M. agrestis splenocytes (Jackson et al., 2011).
RNA was extracted from blood samples stored in RNAlater at −70°C using the Mouse RiboPure Blood RNA Isolation Kit (ThermoFisher), according to the manufacturer's instructions. RNA extracts were DNAse treated and converted to cDNA using the High-Capacity RNA-to-cDNA kit (ThermoFisher), according to the manufacturer's instructions, including reverse transcription negative (RT-) controls for a subsample. Assays were pipetted onto 384 well plates by a robot (Pipetmax, Gilson) using a custom program and run on a QuantStudio 6-flex real-time PCR System (ThermoFisher) at the machine manufacturers default real-time PCR cycling conditions. Reaction size was 10 μl, incorporating 1 μl of template and PrecisionFAST qPCR Master Mix with low ROX and SYBR green (PrimerDesign) and primers at the machine manufacturer's recommended concentrations. We used three standard plate layouts for assaying, each of which contained a fixed set of target gene expression assays and the two endogenous control gene assays (the same sets of animals being assayed on matched triplets of the standard plate layouts). Unknown samples were assayed in duplicate wells and calibrator samples in triplicate wells and no template controls for each gene were included on each plate. Template cDNA (see above) was diluted 1/20 prior to assay. We used as a main calibrator sample a single pool of cDNA derived from the splenocytes of several hundred Kielder M. agrestis samples collected in (Jackson et al., 2011. As sesn3, tollip, il1rap, irf2 and pgrmc1 were relatively poorly represented in this main calibrator sample, a synthesized 478 bp gene fragment containing the amplification target for each of these genes was used as an additional calibrator sample (at 10 × 10 5 copies/μl) in these cases. Samples from different sampling groups were dispersed across plate triplets, avoiding confounding of plate with the sampling structure. Gene relative expression values used in analyses are RQ values calculated by the QuantStudio 6-flex machine software according to the ∆∆Ct method, indexed to the calibrator sample. Melting curves and amplification plots were individually inspected for each well replicate to confirm specific amplification.

| Microparasite assays
We quantified infections by microparasites (B. microti and Bartonella spp.) in blood using a strategy analogous to the host gene expression assays above, targeting pathogen ribosomal RNA gene expression and normalizing to the host endogenous control genes. We included two extra sets of primers in the blood qPCR assays described above (in one of the standard plate layouts for each plate triplet above). For

| Gene expression in splenocytes
Spleens were removed by dissection immediately after culling and stored temporarily in RPMI 1640. Spleen cells were disaggregated using a 70 μm cell strainer and treated with erythrocyte lysis buffer.
The remaining leucocytes were washed three times in RPMI 1640 and plated in duplicate at 2 × 10 6 cells/ml. The culture medium was supplemented with 10% heat-inactivated foetal calf serum, 2 mM lglutamine, 100 μg/ml penicillin, 100 μg/ml streptomycin and 60 μm monothioglycerol. Cells were cultured at 37°C and 5% CO 2 for 48 h, and then stored in RNAlater at −80°C. Note that the culture period was important as these cells acted as the controls for further samples that received stimulation during culture, but those data are not analysed for this study; see Wanelik et al. (2020) for further details.
Levels of gene expression were analysed in the cultured cells for a panel of 21 genes, again with the selection of genes based on known immune function and sensitivity to environmental or intrinsic host variables in our previous work. This led to the selection of a distinct, but overlapping, panel when compared to the blood assays (10 genes were common to both panels), due to differences in gene expression patterns in the two tissues. Primers (23 sets, including the endogenous control genes ywhaz and sdha, which were selected as a stable pairing from our previous stability analysis of candidate control genes in M. agrestis splenocytes (38)) were designed de novo and supplied by PrimerDesign (16 sets) or designed de novo in-house (7 sets). All primers (including those from PrimerDesign) were validated to confirm specific amplification and 100 ± 10% PCR efficiency under our assay conditions. We extracted RNA from the splenocytes using the RNAqueous micro total RNA isolation kit (ThermoFisher), following the manufacturer's instructions. We then followed the same qPCR protocol as described above for the blood samples, to obtain relative expression values for each gene. We used two standard plate layouts for assaying, each of which contained a fixed set of target gene expression assays and the two endogenous control gene assays (the same sets of animals being assayed on matched pairs of the standard plate layouts). In this case, the calibrator sample was the main calibrator sample used for gene expression assays in the blood (above) but without the addition of a synthesized gene fragment.

| Macroparasites
Small and large intestines were removed by dissection immediately after culling and stored in ethanol at room temperature. Subsequently, the intestines were dissected and examined for the presence of macroparasites, which were counted and identified from visual features. Gut parasites identified were nematodes Heligmosomoides laevis, Trichuris arvicolae and Syphacia nigeriana, and tapeworms (Cestoda: Hymenolepidae and Anoplocephalidae, not identified to species level). Of these, H. laevis and T. arvicolae had very low prevalence (<10%) and therefore were not included in statistical analyses.
After dissection, we also searched for ectoparasites by a thorough visual inspection of the fur under a dissecting microscope. We recorded counts of ticks and fleas, being potentially important vectors with relevance to B. microti and Bartonella spp. Several species of small mites were also recorded but not included in statistical analysis here.

| Derived variables
Outputs from qPCR were converted into fold changes in gene expression relative to the calibrator sample by the ΔΔCT method described by Vandesompele et al. (2002). These relative quantification (RQ) values were log-transformed to improve normality of residuals from linear models (see below). Zero values were treated separately during modelling (see below).
Parasite variables were converted to binary measures of presence/absence. For B. microti, and to a lesser extent Bartonella spp., distributions of nonzero RQ values were bimodal with a secondary peak very close to zero. Therefore, for these two parasites the cutoff was drawn at a RQ value of 1, which we judged to narrow "presence" down to biologically meaningful examples.
Seasonal variation was modelled using two sinusoid curves (sine and cosine) with a period of one year. Including these two variables in a model allows for a curve to be fitted that varies smoothly between a maximum at any point in the year and a corresponding minimum 183 days apart. The two seasonal variables were considered as one for the purposes of model simplification. More detailed parameterisation of seasonal variation would have been possible but was not considered constructive here given (a) the lack of data from winter months and (b) that seasonality was not the primary focus of the study.
Scaled mass index (SMI; Peig & Green, 2009) was calculated from body mass and length measures as an index of body condition.
The continuous variables interval length, SMI, body length and lens weight were each scaled to a mean of 0 and standard deviation of 1 when used as predictors in models.
We simplified our categories of reproductive status to "active" or "inactive" for both males and females. Males were classed as active if they had descended testes, and females if they were pregnant, lactating or had a perforate vagina.

| Statistical modelling
Analyses were carried out in R version 3.6.1 (R Core Team, 2019), making use of the tidyverse packages for data processing and graphical presentation (Wickham, 2016). Recognizing that expression levels will be correlated across multiple genes, we considered the option of adopting multivariate methods for analysing gene expression variation. Unfortunately, missing values scattered across genes meant that there were very few complete cases when all genes were considered together (only 41 of 544). Techniques such as multiple imputation (Josse et al., 2011) can work around the missing data but at the cost of reducing the signal to noise ratio. Exploring the complete cases using principal components analysis (PCA) revealed only moderate gains from dimension reduction when compared against randomized examples (n = 1000) in which any correlational structure between variables was removed (Björklund, 2019). The first three components explained 32% of the variation, as compared to a median of 25% in the randomized null models.
Instead, we adopted a univariate approach, with each gene from the blood or spleen panels analysed with a separate generalized linear mixed model (GLMM) implemented in package lme4 (Bates et al., 2015). The use of GLMMs enabled the use of some extra features, such as random predictors and hurdle models, which are difficult or impossible to apply to multivariate analyses. The main drawback of our approach is that after correction for multiple testing, our results will be conservative since repeated tests on correlated variables will contain overlapping information. However, as demonstrated in our exploration of PCA results above, the degree of overlap is not great.
We used QQ plots to assess normality of residuals, scale-residual plots to assess homoscedasticity, and observed-residual plots to check for nonlinear relationships. Using the transformations described under "derived variables" above, we did not find any major deviation from these assumptions.
We initially constructed full models with all fixed and random effects included. First for random and then for fixed terms, we tested all possible combinations of terms, and discarded any combinations with ΔAIC >2 from the model with the lowest AIC score. We then averaged the remaining set of "best" models using package MuMIn (Barton, 2016) to produce our final simplified model.
Within a set of models, we corrected for multiple tests of the same predictor across multiple gene expression variables using the Benjamini-Hochberg method (Benjamini & Hochberg, 1995) to produce "q values". We used a conservative threshold of q < 0.05 for significance, which indicates a false discovery rate of 5%.

| Cross-sectional gene expression models
We fitted hurdle GLMMs with gene expression variables as the response. For genes with >15% zero values, we fitted a binomial model (the "hurdle") for whether or not the gene was expressed at all. Then we fitted a separate linear mixed model (i.e., Gaussian model) to the non-zero values only, provided N > 100. Fixed effects were sex, lens weight (proxy for age), reproductive state (active or inactive), SMI (proxy for body condition), infection with parasites B. microti, Bartonella spp., cestodes and S. nigeriana, year (categorical), season (continuous; see "derived variables" above) and number of days captive (categorical: 1 or 2). Random effects were site and sampling session.

| Cross-sectional PCV models
We selected all individuals recorded as positive for B. microti and fitted a linear mixed model with PCV squared as the response (this transformation yielded residuals closer to a normal distribution). To keep numbers of predictors low, we selected a single representative gene expression variable in light of the general gene expression models above. We selected irf2 as having the strongest positive association with B. microti in the cross-sectional blood data (Figure 2), as well as strong correlations with other important variables tollip and irf9 (Pearson's r = .80 and 0.55, respectively, both p < .0001).
Other fixed effects were intensity of B. microti infection (qPCR value), sex, condition (SMI) and age (lens weight). We included an interaction between B. microti intensity and irf2 expression to test whether individuals that express irf2 more strongly suffer less anaemia for equivalent infection intensities, which would support a tolerance function (Råberg, 2014). Random effects were site and sampling session.

| Longitudinal models
We transformed longitudinal data to refer to time intervals between two successive records for a single individual. We then discarded any intervals greater than 20 days in length (i.e., that did not correspond to an adjacent pair of sessions). This was important for standardization as well as to minimize our chances of missing an important change in status (e.g., new infection, pregnancy). Due to high recapture rates, 82% of intervals were 20 days or less and were retained.
For models of B. microti acquisition, we then selected only intervals when the individual was free from infection at the start, and categorized them according to whether or not an infection had been acquired by the end of the interval. This binary category was then used as the response variable in models. Fixed effects were Bartonella spp. infection, sex, reproductive state, body length (proxy for age as lens mass not available), SMI, tick and flea presence, interval length, season and year. Random effects were site and sampling session. All were measured from the start of the interval.
We initially attempted to run models on all such intervals using individual ID as a random effect, however these models did not converge effectively. Therefore, we instead subsampled the data randomly so that only one interval was retained per individual, thus avoiding problems of pseudoreplication. To mitigate the extra randomness introduced by this approach, we repeated models on 100 different randomly drawn subsamples, and report the proportion of those repeats in which a predictor had q < 0.05.
An initial model of B. microti acquisition was run without any immune-related predictors. After model simplification, which left only year as a significant effect, each gene expression variable was added separately to this base model. This strategy was to avoid problems of overfitting and correlations among predictors if all gene expression variables were added at once.
For models of change in gene expression, we calculated log fold change in expression of a given gene from the start to the end of the interval, for use as our response variable. It was necessary to exclude intervals that included zero values at either the beginning or end, since it is not possible to calculate log fold change in those instances. This led us to reject seven genes which had too high a proportion of intervals containing a zero value, and accounted for 38% of intervals in the remaining genes. Fixed effects were B. microti and Bartonella infection status, reproductive status, sex, body length, interval length, season and year. Random effects were site, session and individual ID. To confirm that the exclusion of zero values was not strongly affecting our results, we also conducted a simpler parallel analysis in which changes were categorized as increasing, decreasing or stable (<10% change), which allowed the inclusion of zero data. Using this ordinal response variable, we fitted cumulative link mixed models equivalent to the above LMMs, using package ordinal (Christensen, 2019). Results were closely comparable with those presented in the main results section (Table S7).

| Expression of immune related genes is associated with B. microti infection and temporal variation
From 544 culled field voles across 3 years and seven sites, we measured expression of panels of genes in the spleen (21) and blood (28) as a representation of systemic immune state. These genes were selected based on known importance to immune function in mice and in other data sets from these vole populations (see Section 2), including various cytokines, cell markers and transcriptions factors (Table 1). For each separate gene, we tested the effects of various predictors including a range of parasite (presence/ absence of B. microti, Bartonella spp., cestodes and S. nigeriana), host (sex, age, reproductive state and body condition) and environmental (site, season and year) variables, as well as a technical variable of length of time in captivity, since voles were kept for either 1 or 2 days prior to culling. Among these predictors, the broadest effects on immune state in both the blood and the spleen were from infection with B. microti (significant effects on expression of 8 of 28 genes in the blood, Table S1, and 9 of 21 in the spleen, Table S2) and temporal variation (season and year between them affecting 22 of 28 in the blood, Table S1, and 18 of 21 in the spleen, Table S2). Host variables and infection with Bartonella spp. were associated with expression of only small numbers of genes in the blood (≤4 each, Table S1) and none in the spleen (Table S2), and no effects of gut parasites on systemic immune state were detected (Tables S1 and S2, Figure 1). For 10 genes which were included in both the spleen and the blood panels, only two (apobr and foxp3) showed significant correlations between expression levels in the two tissues (Table S3).
Babesia microti was among the most prevalent of the observed parasites (Table 2) and was generally associated with increases in expression of anti-inflammatory and antioxidant genes in the blood (Figure 2). In the spleen there was both up-and downregulation of genes of a range of functions, including inflammatory and macrophage-related genes ( Figure 2, Table 1, Table S2).

| Individuals infected with B. microti exhibit anaemia
Looking for possible negative effects of B. microti on host physiology, we used PCV as a measure of erythrocyte numbers. Individuals had lower PCV, indicating anaemia, when infected with B. microti (p = .012). PCV correlated negatively with intensity of B. microti infection (Figure 3, Table S4), indicating anaemia. We then asked whether this harm associated with the infection showed any correlation with immune state. PCV was negatively associated with expression of irf2 (interferon regulatory factor 2; picked as a representative of B. microti-associated immune state as it was the top TA B L E 1 Genes comprising the blood and spleen panels assayed by quantitative PCR (qPCR) Note: B and/or S indicate inclusion in the blood and/or spleen panels, respectively."Type" and "Immune relevance" columns refer to the products coded for by the listed gene.
fold-difference gene in the blood), even after accounting for infection levels ( Figure 3, Table S4). However, the interaction term between infection intensity and irf2 expression was not retained after model simplification (ΔAIC = 2.14), providing no evidence that individuals with higher irf2 expression were able to tolerate greater infection intensities than those for which expression was low.

| Gene expression changes associated with B. microti infection occur within longitudinal samples from individual voles
To gain some insight into the possible causal relationships between B. microti infection and immune state, we then used longitudinal data F I G U R E 1 Heatmap showing the modelled effects (log fold change from linear mixed models) of the listed predictors on expression of the spleen and blood gene panels. Nonsignificant effects (q > 0.05) are left blank. Symbols indicate significant effects in the equivalent zero model (logistic generalized linear mixed model), comparing "expressed" with "not expressed", with "+" indicating a positive effect on expression (i.e., fewer zeros) and "-" indicating a negative effect. Grid lines are for ease of interpretation only.  TXNRD2  TOLLIP  TGFB  SESN3  RETNLG  PGRMC1  ORAI1  NOS2  IRF9  IRF2  INSR  IL4  IL1RN  IL1RAP  IL1B  IL17  IL10  IFNG  GPX1  GATA3  FOXP3  COL1A1  CD99L2  CD8A  CD86  CCN2  AR  APOBR   TNFA   TGFB   TBET  In light of the cross-sectional analysis, we tested whether the associations between B. microti and gene expression could be explained by individuals with certain immune profiles having an increased chance of acquiring an infection. Each gene expression variable was used as a predictor for B. microti acquisition in separate generalized linear models alongside host, parasite and temporal F I G U R E 2 Differences in gene expression associated with Babesia microti infection in cross-sectional samples. Points show mean log fold difference in infected individuals against those uninfected, and ranges show 95% confidence intervals from 1000 bootstrap samples. Genes are ordered by effect size. Red text and triangle point symbols indicate genes that appear in both tissue panels to help with cross-references between the two panels. Greyed columns indicate changes in expression that were not significant, with q > 0.05, in models that included additional predictors beyond infection status.
We then looked for variables that correlated with changes in gene expression within an individual. Nine genes changed in expression when individuals either acquired a new infection of B. microti or cleared an existing one (Table S6). For all but one of these genes (retnlg), directions of estimated effects were reversed when the infection was cleared compared to when it was acquired (Figure 4). Relative expression of IRF2 F I G U R E 4 Changes in gene expression when Babesia microti infection status changes. Points show the mean log fold change in expression between two trapping sessions (<20 days) within an individual; vertical lines indicate standard error. "Acute" infection status indicates an individual that was free from infection at the start of the interval but infected at the end, and "cleared" indicates the opposite. Genes are ordered by "acute" mean values. Greyed columns indicate changes in expression that were not significant, with q > 0.05 for both acute and cleared intervals, in models that included additional predictors beyond infection status. See Table S6 for further details. Other predictors, including season, were not associated with changes in gene expression within an individual (Table S6). For confirmation, we also ran a parallel analysis in which, rather than being treated as continuous, changes were categorized as increasing, decreasing or stable (see Section 2). The results from this simpler analysis were broadly comparable with those presented above (Table S7).

| DISCUSS ION
Our results show that wild field voles infected with B. microti un- microti is a cause, rather than a consequence, of the observed differences in immune state. We also demonstrate that these changes persist throughout the potentially long period of infection, only reverting to base levels on the relatively rare occasions when the infection was cleared. We therefore establish this prevalent chronic infection as an important driver of immune state in our study population.
Babesia microti can have substantial physiological effects on its host, most notably causing major splenomegaly (Lykins et al., 1975;Watkins et al., 1991), and in particular the proliferation of macrophages in the spleen (Cullen & Levine, 1987). These responses are probably part of a coping strategy, as asplenic hosts typically suffer worse disease (Krause et al., 2008;Lykins et al., 1975). It is therefore unsurprising that we see shifts in gene expression in the spleen, likely related to changes in proportions of cell types and their states of activation. For example, the marked increase in expression of marco, which codes for a protein macrophage receptor, probably reflects an increase in the proportion and activity of macrophages, which may be involved in clearance of B. microti (Djokic et al., 2018).
Several other genes undergo changes that may be related to antiinflammatory functions (upregulation of il10 and downregulation of il1b and il6).
In the blood, immune state is also impacted by B. microti, although correlation between genes expressed in the two tissues is weak (Table S3), indicating that we are not simply seeing a "spillover" effect of cells and cytokines from the spleen. Genes with antiinflammatory (tollip) and antioxidant (sesn3, txnrd2, gpx1) roles are among those with the strongest positive association with B. microti infection. We note that while expression of ifng, known to play a role in clearing acute infections of B. microti (Igarashi et al., 1999), itself shows no significant change, we see a major increase in expression of irf2 and irf9 in the blood. Both of these are transcription factors involved in the regulation of interferons and interferon-induced gene expression (Antonczyk et al., 2019). Although irf transcription factors primarily directly regulate type I, rather than type II interferons, such as ifng, they may have indirect effects within ifng-driven transcriptional programmes (Waddell et al., 2010) and irf2 has been implicated in helminth modulation of immunity to other apicomplexan blood parasites (malaria; Metenou et al., 2012). Here, any response to the acute infection involving ifng itself was either absent or short-lived (and hence missed by our sampling intervals).
A further role for irf2 is stimulating erythropoiesis (Mizutani et al., 2008). We found that individuals infected with B. microti exhibit mild anaemia, which is also associated with higher irf2 expression.
This may demonstrate an attempt by the host to compensate for the loss of erythrocytes that have been invaded by B. microti. We did not find any evidence that irf2 increases levels of tolerance towards B. microti, as there was no interaction between infection intensity and irf2 expression in their effect on PCV. Causal relationships among infection, anaemia and irf2 are not clear from our data as we were unable to measure PCV from longitudinal data. Elsewhere, increased erythropoiesis has been observed in laboratory mice infected with B. microti (Hussein, 1982), in which instance the response was sufficient to compensate for the erythrocyte loss. In light of the associations in our data and the potential mechanistic link between irf2 and erythropoiesis, we believe this possible tolerance mechanism merits investigation in future studies.
Taken together, the observed changes in the blood suggest that the response is aimed, at least in part, at tolerance of infection and minimizing tissue damage in the blood. A similar signal forms a component of the response to malaria in birds, where genes with functions related to antioxidation were upregulated during infection (Videvall et al., 2015). Our previous study indicates that, at least in field voles, B. microti has remarkably low virulence, and despite the substantial effects on physiology and immunology, there is no clear fitness cost . Anti-inflammatory activity is considered a likely contributor to disease tolerance (Sears et al., 2011), and there is some evidence for both Th2  and anti-inflammatory (Manzoli et al., 2018) pathways contributing to tolerance in wild animals. The idea that the inflammatory response is potentially being actively suppressed during B. microti infection to reduce harmful side-effects is consistent with our results.
Equally, this does not rule out a role for "resistance"-type responses to B. microti, which could act in parallel. Elevated marco expression in the spleen and the increase in irf9 expression in the blood point towards macrophage and interferon activity respectively, both of which could be targeted at lowering parasite numbers.
While the observed immune response is not accompanied by obvious fitness impacts , it may nonetheless be important for the infected individual. Due to their long-term nature, impacts can sometimes be difficult to detect. For example, it was long thought that impact from the chronic phase of avian malaria infection was minimal (Asghar et al., 2012), but more recent work showed longterm fitness impacts associated with telomere shortening (Asghar et al., 2015). Immune responses lead to immunopathology (Bertrand et al., 2006) and energetic costs (Lochmiller & Deerenberg, 2003), but the impact of these costs on fitness might only become apparent when resources are scarce. Alternatively, long-term changes in immune state could lead to occasional but substantial fitness costs through altered susceptibility to other, more virulent infections.
Host variables included in our models (age, sex, condition and reproductive activity) showed very little association with any of our measured markers. In some cases, this is surprising when compared to other similar studies. For example, while some other wild mammal populations show variation in immune state with age (Abolins et al., 2018;Clerc et al., 2019;Watson et al., 2016), here we find very little association with age. A likely explanation is that our data are not well suited to detecting an effect of immunosenescence given that they do not include many individuals at the extremes of the possible age range where effects might be expected to be strongest (Clerc et al., 2019;Nussey et al., 2012). Likewise, in contrast to our limited association between body condition and immune gene expression, Abolins et al. (2018) found strong associations between immune cell counts and condition. Animals with greater relative body mass have greater raw numbers of immune cells (Abolins et al., 2018) but our results suggest that the activity of those cells is similar regardless of condition. Meanwhile, our results agree with numerous other studies in finding significant seasonal variation in our immune markers (Beldomenico et al., 2008;Jackson et al., 2011;Lochmiller et al., 1994;Xu et al., 2018), although as we did not sample during winter months we are limited as to detailed interpretation.
Through a powerful longitudinal data set, our study demonstrates the major impact that chronic infection with B. microti can have on the immune state of a wild rodent. Babesia microti stands out among several other prevalent infections that by contrast appear to have limited influence on the host's immune state. This influence of a chronic haemoparasite is especially noteworthy given the lowvirulence nature of the infection; effects such as these might easily be overlooked unless specifically targeted for study. More broadly, we improve our understanding of the factors that determine an individual's immune state, and ultimately why individuals may differ in

ACK N O WLE D G E M ENTS
We would like to thank Maria Capstick, Noelia Dominguez Alvarez, William Foster, Sarah Gore, Ann Lowe, Lukasz Lukomski, Ed Parker, Benoit Poulin, Stephen Price, Anna Thomason, Rebecca Turner and Susan Withenshaw for their roles in collecting and processing samples from the field. We thank the Forestry Commission for access to the study sites in Kielder Forest. Our work was funded by the Natural Environment Research Council, award NE/L013452/1 to JEB, MB, JJ and SP.

DATA AVA I L A B I L I T Y S TAT E M E N T
Data are available from the NERC EDS Environmental Information Data Centre: https://doi.org/10.5285/e5854 431-6fa4-4ff0-aa02-3de68 763c952 Paterson et al. (2022). Benefit sharing: Benefits from this research accrue from the sharing of our data and results on public databases as described above.