IgE receptor polymorphism predicts divergent, sex-specific inflammatory modes and fitness costs in a wild rodent

Sexual dimorphism is widespread in the animal kingdom, with males and females differing in their physiology, morphology and behaviour. However, relatively little is known about the importance of sex for the expression of genotype in natural populations, particularly in the context of health and disease. We combine data on genotype, immune gene expression, infection incidence and pedigree of individuals in a wild population of field voles (Microtus agrestis). We identify a polymorphism in the Immunoglobulin E (IgE) receptor with sexually dimorphic effects on immune phenotype, health and, ultimately, fitness. While males with the polymorphism upregulate their pro-inflammatory response with a detrimental effect on their reproductive success, females upregulate their anti-inflammatory response, increasing their risk of infection, but with no apparent reproductive cost. Our work demonstrates the potential for the same genetic difference to have profoundly different consequences for the health and fitness of males and females in natural populations.


Introduction
Mutation, and subsequent selection, are thought to play an important role in shaping the immune system, with some of the fastest evolving genes being immune related (Nielsen et al., 2005). However, the effects of any mutation will depend on the context in which it occurs. Sexual dimorphism is widespread in the animal kingdom; with males and females representing different contexts or 'environments' within which genes (or mutations within genes) act. Sex can therefore have an important effect on the phenotypic expression of genotype (genotype by sex interactions). However, the majority of evidence for genotype by sex interactions comes from laboratory experiments conducted on model organisms (e.g. Mackay and Anholt 2006). Relatively little is known about the importance of sex for the expression of genotype in natural populations, particularly in the context of health and disease.
Studies of natural populations have explored the effects of genotype on immune phenotype, with consequences for susceptibility to infection. Most notably, variability in the genes of the Major Histocompatibility Complex has been associated with resistance to intestinal nematodes in domestic sheep (Paterson, Wilson, & Pemberton, 1998) and with resistance to malaria, hepatitis, and AIDS in humans (Carrington et al., 1999;Hill et al., 1991;Thursz, Thomas, Greenwood, & Hill, 1997). The role of polymorphism elsewhere in the genome, for shaping immune phenotype, has also been studied (Turner, Begon, Jackson, & Paterson, 2012;Wanelik et al., 2018). However, in the vast majority of these studies, genotype by sex interactions were not considered. In some studies, particularly those of natural populations, this is due to low sample size (and hence limited power), making it difficult to test for datahungry genotype by sex interactions (Ober, Loisel, & Gilad, 2008). In some studies, sex is considered a confounding variable which is controlled for, in order not to hide any subtle genetic associations (Paterson et al., 1998). Other studies focus on a single sex for the sake of simplicity (e.g. Jackson et al. 2014). Where sex by genotype interactions are documented, the mechanisms are often unknown (e.g. ACE genotype-sex interaction in humans).
However, the sexes differ not only in their morphology and behaviour, but less conspicuously, in their physiology, biochemistry and in their genetic architecture. Sex-biases in the expression of thousands of autosomal genes (present in both males and females) have been documented, suggesting that males and females differ in their regulation of expression (Ellegren & Parsch, 2007;Rinn & Snyder, 2005; X. Yang et al., 2006). Sex-biases in the co-expression of genes have also been documented, suggesting that males and females differ not only in their regulation of single genes, but whole networks of interacting genes (Sutherland, Prokkola, Audet, & Bernatchez, 2018;Van Nas et al., 2009). The same polymorphism then, could be expressed as different immune phenotypes, and subsequently have very different consequences for the health of males and females in the wild. Furthermore, in the real world, males and females differ in the environments (physical and social) they themselves inhabit and the selective pressures they face, all of which may also interact with genotype. One might therefore expect genotype by sex interactions to be even more important. We argue then that by excluding sex, in studies of natural populations especially, one is in danger of underestimating the effect of genetic variation not only on the health, but on the fitness of individuals.
In humans, immunoglobulin E (IgE) is involved in defence against helminths and parasitic worms (Gounni et al., 1994). However, it has also been implicated in inflammatory disease, particularly allergy (Tomassini et al., 1991). In order to effect a response, antigen (or allergen)-bound IgE must bind to its receptor, Fcer1, which is found on the surface of mast cells, basophils, eosinophils, and Langerhans cells. On activation, these immune cells degranulate, releasing histamine and other inflammatory mediators. IgE-mediated responses are therefore controlled by the Fcer1 receptor (Daeron & Nimmerjahn, 2014). Naturally occurring polymorphisms in this receptor are known to affect an individual's serum IgE levels with likely consequences for their susceptibility to infection (Granada et al., 2012;Weidinger et al., 2008) and their risk of developing inflammatory disease (Hasegawa et al., 2003;Niwa et al., 2010;Potaczek et al., 2006;Zhou et al., 2012). Furthermore, serum IgE levels (Weiss, Pan, Abney, & Ober, 2006) and the incidence of IgE-mediated inflammatory disease (e.g. Chen et al. 2008) are both sexually-dimorphic traits, suggesting that any polymorphism in this pathway is likely to experience very different contexts in males and females. Indeed, one study found evidence for a sex-specific effect of a single nucleotide polymorphism (in the alpha chain of the Fcer1 receptor) on susceptibility to systematic lupus erythematosus (a chronic inflammatory disease) in a Chinese Han population (Yang et al. 2013).
In a previous study of a natural population of field voles (Microtus agrestis), we found an association between a polymorphism in the alpha chain of the Fcer1 receptor and an immunological marker of tolerance to infection in males . Here, we show that this polymorphism also has important consequences for resistance to infection in the same population, via its effects on inflammation. We also explore the effects of this polymorphism in females. In doing so, we provide a rare example of polymorphism, in the wild, with sexually dimorphic effects on inflammatory phenotype, resulting in asymmetrical costs in males and females.

Field methods
We studied M. agrestis in Kielder Forest, Northumberland (55°130N, 2°330W), using livetrapping of individual animals from natural populations. Trapping was performed from 2015-2017 across seven different sites, each a forest clear-cut. At each site, 150-197 Ugglan small mammal traps (Grahnab, Sweden) were laid out in a grid spaced approximately 3-5 m apart.
Our study was divided into longitudinal and cross-sectional components. Every other week, traps were checked every day, once in the morning and once in the evening. Newly trapped field voles were injected with a Passive Integrated Transponder (PIT) tag (AVID, UK) for unique identification. This approach allowed us to build up a longitudinal record for voles that were caught on multiple occasions. A total of 2, 881 voles were individually marked in this way. We also took a small tissue sample from the tail, for genotyping of individuals and a drop of blood from the tail which we put into 500 ml of RNAlater (Fisher Scientific, UK), for use in parasite detection (see below). On initial and subsequent captures, we administered general-use anti-helminthic (Profender, Bayer AG) and anti-ectoparasitic treatment (Frontline, Boehringer Ingelheim) to half of the trapped population (those voles with PIT tags ending in an even number, except for pregnant females). All animal procedures were performed with approval from the University of Liverpool Animal Welfare Committee and under a UK Home Office licence (PPL 70/8210 to SP).

Cross-sectional data
For the cross-sectional component of the study, captured animals (n = 544) were returned to the laboratory where they were killed by a rising concentration of CO2, followed by exsanguination. This component of the study allowed us to take a more comprehensive set of gene expression measurements, parasite measurements and biomarker measurements (see below). A small tissue sample was also taken from the ear of cross-sectional animals for genotyping.

Splenocyte cultures
Spleens of cross-sectional animals were removed, disaggregated, and splenocytes cultured under cell culture conditions equivalent to those used in Jackson et al. (2011). Unstimulated splenocytes, taken from 62 cross-sectional animals collected between July and October 2015, were initially used to assay expression by RNASeq (see below). Given the known link between polymorphisms in the Fcer1 receptor and risk of developing inflammatory disease, we then exposed splenocytes from the same animals to inflammatory stimulation in vitro, namely the TLR2 agonist heat-killed, Listeria monocytogenes (HKLM; concentration 5 × 107 cells ml-1; Invivogen, Paisley, UK), and the TLR7 agonist, imiquimod (IMI; concentration 10 μg ml-1; Invivogen). Stimulated splenocytes were used to assay expression by Q-PCR.

RNASeq
Full details of the methods used for RNA preparation and sequencing can be found in Wanelik et al. (2018). Briefly, samples were sequenced on an Illumina HiSeq4000 platform.
High-quality reads were mapped against a draft genome for M. agrestis (GenBank Accession no. LIQJ00000000) and counted using featureCounts (Liao, Smyth, & Shi, 2014). Genes with average log counts per million, across all samples, of less than one were considered unexpressed and removed from the data (n = 8, 410). Following filtering, library sizes were recalculated, data were normalized and MDS plots were generated to check for any unusual patterns in the data.

Q-PCR
We used two-step reverse transcription quantitative real-time PCR (Q-PCR) to measure the expression levels of a panel of 20 genes in HKLM-and imiquimod-stimulated splenocytes from our cross-sectional animals. We did this, in part, to validate our RNASeq results in an independent dataset. We used the observed expression profile as a measure of the potential responsiveness of the immune system to an inflammatory stimulation in vivo. The choice of our panel of genes was determined by an important known immune-associated function in mice, combined with significant sensitivity to environmental or intrinsic host variables in our previous studies (Jackson et al., 2011(Jackson et al., , 2014 or in a recent differential expression analysis (DEA) of RNASeq data (not reported here).
Primers (23 sets, including 2 endogenous control genes) were designed de novo and supplied by PrimerDesign (Chandler's Ford, UK) (15 sets) or designed de novo in-house (8 sets) and validated (to confirm specific amplification and 100 ± 10% PCR efficiency under assay conditions). All PrimerDesign primer sets were validated under our assay conditions before use. The endogenous control genes (Ywhaz and Sdha) were selected as a stable pairing from our previous stability analysis of candidate control genes in M. agrestis splenocytes (Jackson et al. 2011). We extracted RNA from splenocytes conserved in RNAlater using the RNAqueous Micro Total RNA Isolation Kit (ThermoFisher), following the manufacturer's instructions. RNA extracts were DNAse treated and converted to cDNA using the High-Capacity RNA-to-cDNA™ Kit (ThermoFisher), according to manufacturer's instructions, including reverse transcription negative (RT-) controls for a subsample. SYBR green-based assays were pipetted onto 384 well plates by a robot (Pipetmax, Gilson) using a custom programme and run on a QuantStudio 6-flex Real-Time PCR System (ThermoFisher) at the machine manufacturers default 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 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). 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. The calibrator sample (identical on each plate) was created by pooling cDNA derived from across all splenocyte samples. Samples from different sampling groups were dispersed across plate pairs, 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.

Parasite detection
After culling, the fur of cross-sectional animals was examined thoroughly under a binocular dissecting microscope to check for ectoparasite infections, which were recorded as direct counts of ticks and fleas. Guts of cross-sectional animals were transferred to 70% ethanol, dissected and examined under a microscope for gastrointestinal parasites. Direct counts of cestodes were recorded.
In the longitudinal component of our study, we quantified infections by microparasites (Babesia microti and Bartonella spp.) in blood using a strategy analogous to the SYBR green host gene expression assays above and targeting pathogen ribosomal RNA gene expression normalized to the expression of host endogenous control genes. Blood samples were derived from tail bleeds. RNA was extracted from blood samples stored in RNAlater at -70°C using the Mouse RiboPure Blood RNA Isolation Kit (ThermoFisher), according to manufacturer's instructions, and DNAse treated. It was then converted to cDNA using the High-Capacity RNA-to-cDNA™ Kit (ThermoFisher), according to manufacturer's instructions. For B.
microti we used the forward primer TBCbabF and reverse primer TBCbabR targeting the 18S ribosomal RNA gene and for Bartonella spp. we used the forward primer TBCbartF and reverse primer TBCbartR targeting the 16S ribosomal RNA gene. Assays were pipetted onto 384 well plates by a robot (Pipetmax, Gilson) using a custom programme and run on a QuantStudio 6-flex Real-Time PCR System (ThermoFisher) at the machine manufacturers default real-time PCR cycling conditions. Reaction size was 10 μl, incorporating 1 μl of template (diluted 1/20) and PrecisionFAST qPCR Master Mix (PrimerDesign) with low ROX and SYBR green and primers at the machine manufacturer's recommended concentrations.
Alongside the pathogen assays we also ran assays for two host genes (Ywhaz and Actb) that acted as endogenous controls. We used as a calibrator sample a pool of DNA extracted from indexed to the calibrator sample. Melting curves and amplification plots were individually inspected for each well replicate to confirm specific amplification. We validated our diagnostic results by comparing our PCR RQ values to independent data for a subset of crosssectional voles with mapped genus-level pathogen reads from RNAseq analysis of blood samples (n = 44), finding the two data sets strongly corroborated each other (Bartonella, r2 = 0.81, p < 0.001; Babesia, r2 = 0.81, p < 0.001).

SOD1 measurement
Given the well-established link between inflammation and oxidative stress, we measured antioxidant superoxide dismutase 1 (SOD1) enzymatic activity in blood samples taken from cross-sectional animals. Assays were carried out using the Cayman Superoxide Dismutase kit and following manufacturer's instructions except where otherwise indicated below. Blood pellets from centrifuged cardiac bleeds were stored at -70°C and thawed on ice prior to assay.
A 20 μl aliquot from each pellet was lysed in 80 μl of ultrapure water and centrifuged (10000 G at 4°C for 15 minutes) and 40 μl of the supernatant added to a 1.6 × volume of ice cold chloroform/ethanol (37.5/62.5 (v/v)) (inactivating superoxide dismutase 2). This mixture was then centrifuged (2500 G at 4°C for 10 minutes) and the supernatant removed and immediately diluted 1:200 in kit sample buffer. A seven-point SOD activity standard was prepared (0, 0.005, 0.010, 0.020, 0.030, 0.040, 0.050 U/ml) and the assay conducted on a 96 well microplate with 10 μl of standard or sample, 200 μl of kit radical detector solution and 20 μl of xanthine oxidase (diluted as per manufacturer's instructions) per well. Plates were then covered with a protective film, incubated at room temperature for 30 minutes on an orbital shaker and read at 450 nm on a VERSAmax tuneable absorbance plate reader (Molecular Devices), subtracting background and fitting a linear relationship to the standard curve in order to estimate SOD activity in unknown samples.

Genotyping
We identified 346 single nucleotide polymorphisms (SNPs) in 127 genes. See Wanelik et al. (2018) for details of the approach used to select these SNPs. Our list included two synonymous SNPs in the immune gene Fcer1a (alpha chain of the high-affinity receptor for IgE) which we had previously identified as a candidate tolerance gene in a natural population of M. agrestis. We identified three haplotypes present in the population GC, AC and AT, at frequencies of 0.12, 0.76 and 0.07 respectively. We also identified the GC haplotype as being of particular interest, given its significantly lower expression level of the transcription factor Gata3 (a biomarker of tolerance to infection in our population) compared to the other haplotypes . DNA was extracted from an ear sample (cross-sectional component) or tail sample (longitudinal component) taken from the vole using DNeasy Blood and Tissue Kit (Qiagen). Genotyping was then performed by LGC Biosearch Technologies (Hoddesdon, UK; http:// www.biosearchtech.com) using the KASP SNP genotyping system. This included negative controls (water) and duplicate samples for validation purposes.

Pedigree reconstruction
A pedigree was reconstructed using the sequoia package (Huisman, 2017) and a subset of our SNP dataset. SNPs which violated assumptions of Hardy Weinberg equilibrium were removed from the dataset. For pairs of SNPs in high linkage disequilibrium (most commonly within the same gene), the SNP with the highest minor allele frequency (MAF) was chosen.
A minimum MAF cut-off of 0.1 and call rate of > 0.7 was then applied, and any samples for which more than 50% of SNPs were missing were removed. This resulted in a final dataset including 113 SNPs for 2, 674 samples.
Life history information, namely sex and month of birth, for samples was also inputted into sequoia where possible. Juvenile voles weighing less than 12 g on first capture were assigned a birth date two weeks prior to capture. Juvenile voles weighing between 12 and 15 g on first capture were assigned a birth date four weeks prior to capture. Finally, adult voles breeding on first capture, were assigned a birth date six weeks prior to capture (minimum age at first breeding; Begon et al., 2009;Burthe et al., 2010). Adult voles not breeding on first capture could not be assigned a birth date, as it was not known whether they had previously bred or not. Virtually all samples (99%) were assigned a sex, and approximately half (54%) were assigned a birth month. As we sampled individuals from across seven different clear-cut areas of the forest, each several kilometres apart, these were assumed to be independent, closed populations with negligible dispersal. Site-specific pedigrees were therefore generated.
The resulting pedigrees included 1, 211 known individuals with some familial relationship.
As expected from the biology of M. agrestis, the majority of predicted parent-offspring pairs were born in the same year (87%) and were found in similar areas (i.e. along the same transects; 90%). There was no difference in the probability of appearing in a pedigree for individuals with versus without our haplotype of interest (Chi-squared test: 2 = 0.09, df = 1 and p = 0.76). For each individual present in a pedigree, the number of offspring was counted to provide a measure of their reproductive success. Very few individuals were first trapped as juveniles (n = 167), with the majority trapped as adults which had already recruited into the population. Our measure of reproductive success then, more closely resembles the number of recruited (rather than new-born) offspring per individual. Although we sampled a large proportion of the total population within clear-cuts, our sampling was not exhaustive, resulting in incomplete pedigrees. For example, over half of individuals present in our pedigrees (n = 671) were found to have no offspring. Although some of these are true zeros, others may be false i.e. individuals that produced offspring which were not detected. We therefore consider our measure of reproductive success to be relative (rather than absolute).

Statistical analyses
Fcer1a haplotypes were inferred from SNP data, and tested for associations with (i) gene expression, (ii) parasite infection, (iii) reproductive success and (iv) SOD1 concentration. Not all individuals appeared in all datasets, therefore sample sizes vary between analyses.
Throughout, residuals from regression models were checked for approximate normality and homoscedasticity and all non-genetic covariates were tested for independence using variance inflation factors (VIFs). All analyses were performed in R statistical software version 3.5.2 (R Core Team, 2018).

Differential expression analysis
Differential expression analysis was performed on filtered count data using the edgeR package (Robinson, McCarthy, & Smyth, 2010), the aim being to identify individual genes which are differentially expressed between those individuals with and without a copy of the GC haplotype (hGC). Samples from different sexes were analysed separately. As this was an exploratory analysis, used in conjunction with more targeted measurement of immune gene expression (see below), model specification was kept as simple as possible. Only site and season were included as additional covariates, in order to account for any spatial and/or temporal autocorrelation.
Given the known link between polymorphism in the Fcer1 receptor and risk of developing inflammatory disease, we tested for enrichment of pro-and anti-inflammatory genes in our results; more specifically, the gene ontology terms, 'positive regulation of inflammatory response' (GO:0050729; n = 143) and 'negative regulation of inflammatory response' (GO:0050728; n = 149). The aim here was to answer the question are pro-or antiinflammatory genes more highly ranked relative to other genes, when we compare individuals with and without the haplotype? This enrichment analysis was performed using the limma package (Ritchie et al., 2015) and genes were ranked on log fold change.

Association between Fcer1a haplotype and immune gene expression
Initially, we ran multiple models to look for associations between Fcer1a haplotype and immune gene expression using the hapassoc package (Burkett, Graham, & McNeney, 2006).
The hapassoc package allows likelihood inference of trait associations with SNP haplotypes and other attributes, adopts a generalized linear model framework and estimates parameters using an expectation-maximization algorithm. If the haplotype combination of an individual cannot, with certainty, be inferred from its genotyping data (i) because it is heterozygous at two or more markers or (ii) because it has missing data for a single marker, the approach implemented in hapassoc is to consider all possible haplotype combinations for that individual. Standard errors accounting for this added uncertainty are calculated using the Louis' method (Louis, 1982). Hapassoc models (here and below) accounted for other covariates, assumed an additive genetic model and were run separately for males and females.
Other covariates accounted for in hapassoc models and considered potential drivers of immune gene expression, were informed by our previous work (Jackson et al., 2011). These included snout-vent length (SVL), age estimated by eye lens weight, whether or not an individual was given anti-parasite treatment, reproductive status (males: scrotal or non-scrotal testes; females: pregnant or not), body condition (estimated by regressing body weight against SVL and its quadratic term; males only), site, year and season. As we ran a total of 20 models per immune agonist (one model per gene), the Benjamini and Hochberg method of correction was applied to all p-values (Benjamini & Hochberg, 1995). Resulting q-values (FDR-corrected p-values) are reported in the text, alongside original p-values.
We were unable to include any random variables in the initial haplotype association analysis due to the framework adopted by hapassoc (Burkett et al., 2006). So, in order to confirm these results, linear mixed effects models (LMMs), including a random term for assay plate number, were run for those immune genes for which expression appeared to be associated with haplotype (q ≤ 0.1). Immune gene expression variables were also Yeo-Johnson transformed (Yeo & Johnson, 2000) to achieve more normal and homoscedastic residuals (as in Wanelik et al. 2018). Genotype was coded as the presence or absence of the GC haplotype, previously identified as being of interest . Therefore, only those individuals for which haplotype could be inferred with certainty could be included.

Association between haplotype and parasite infection
The three macroparasite measurement taken from cross-sectional animals (counts of ticks, fleas and cestodes) were summarised as a single principal component (explaining 37% of total variation). This combined measure of macroparasite burden was modelled using the hapassoc package. See Jackson et al. (2014) and Wanelik et al. (2018) for full details of this approach. As microparasite infection status (infected or not) was assessed multiple times for the majority of individuals in the longitudinal component of the study (mean = 2.8; range = 1-11), a generalised linear mixed-effects model (GLMM) with a binomial error distribution and a random term for individual, was used to test for an association between Fcer1a haplotype and the probability of an individual being infected with a microparasite. Other covariates, considered potential drivers of parasite infection were, again, informed by our previous work (Taylor et al., 2018). These included year and site.

Association between Fcer1a haplotype and SOD1 concentration
SOD1 concentration was modelled using the hapassoc package. Other covariates, considered potential drivers of SOD1 concentration, were the same as those listed under 'Association between Fcer1a haplotype and immune gene expression'.

Association between Fcer1a haplotype and reproductive success
Initially, we ran a model in hapassoc to look for an association between Fcer1a haplotype and our measure of reproductive success. We included birth month as a covariate in this model, given that autumn-born voles have been shown to have a lower chance to reproduce than spring-born voles (Wang, Ebeling, & Hahne, 2019). Other covariates included in this analysis were, whether or not an individual was culled for the cross-sectional component of this study (again, reducing the opportunity to reproduce), whether or not an individual was given antiparasite treatment, site and year. However, reproductive success was found to be zeroinflated and, due to the framework adopted by hapassoc (Burkett et al., 2006), we were unable to account for this. So, to confirm these results, negative binomial, zero-inflated Poisson (ZIP) and zero-inflated negative binomial (ZINB) generalised linear models (GLMs) were run instead. These models were compared on AIC and the normality and homoscedasticity of residuals. The negative binomial GLM was found to be the best fit.

Males with GC haplotype have a stronger pro-inflammatory response, lower reproductive success and higher concentration of SOD1
Differential expression analysis followed by gene set enrichment performed on the RNASeq data showed that pro-inflammatory genes tended to be upregulated in males with hGC compared to males without (p = 0.02; Fig 1A). More specifically males with hGC expressed the pro-inflammatory cytokine, IL33, at a level nearly seven times higher than males without hGC (rank = 3; logFC = 2.8; q = 0.06; Fig 1B). Further confirmation of the role of hGC in modulating the pro-inflammatory phenotype of males, was provided by the Q-PCR data which indicated that, on exposure to an inflammatory stimulus (IMI), Gata3 expression levels were significantly associated with the GC haplotype (p < 0.01; q = 0.1; n = 231). More specifically, males with hGC expressed Gata3 at a lower level than males without hGC (Fig.   1C).
We then tested the consequences of the GC haplotype for reproductive success and probability of infection. Strikingly, males with at least one copy of hGC had 0.62 times lower reproductive success than males without (p = 0.03; n = 233). However, there was no effect of hGC on the probability of a male being infected by any parasite (e.g. Babesia microti; p = 0.86; n = 470; Fig 1D). The GC haplotype was also associated with 1.3 times higher levels of the oxidative stress biomarker SOD1 (3.27 U/ml blood; SE = 1.30) than the AT haplotype (2.51 U/ml blood; SE = 1.21; p = 0.03; n = 235) in males. represent 95% CI; significance codes: n.s. non-significant . ≤ 0.1 * < 0.05 ** < 0.01).

Females with GC haplotype have a stronger anti-inflammatory response, higher risk of infection but the same reproductive success
Anti-inflammatory genes tended to be upregulated in females with the GC haplotype compared to females without (p < 0.01; Fig 2A). The top-ranking gene for females was the anti-inflammatory suppressor of cytokine signaling SOCS3 (logFC = 1.2; q = 0.16; Fig 2B).
On exposure to an inflammatory stimulus (HKLM), IL10 and ORAI1 expression levels were significantly associated with Fcer1a haplotype (p < 0.01; q = 0.04; n = 100). More specifically, females with hGC expressed IL10 and ORAI1 at a lower level than females without hGC ( Fig 2C).
As for males, we tested the consequences of the GC haplotype for reproductive success and probability of infection. There was an effect of hGC on the probability of a female being infected -females with at least one copy of hGC were 1.5 times more likely to be infected with Babesia microti than females without (p < 0.01; n = 421; Fig 2D). There was no effect of GC haplotype on the reproductive success of females (p = 0.38; n = 417; Fig 2E). Association of hGC with reproductive success (error bars represent 95% CI; significance codes: n.s. non-significant . ≤ 0.1 * < 0.05 ** < 0.01).

Discussion
In this study, we have found a rare example of a polymorphism, in the wild, with sexually dimorphic effects on immune phenotype, health and fitness. While males with the GC haplotype upregulate their pro-inflammatory response with a detrimental effect on their reproductive success, females upregulate their anti-inflammatory response, increasing their risk of infection, but with no apparent reproductive cost.
Our work also provides a rare example of a single locus conferring both resistance and tolerance in the wild. In our previous study, we showed that the GC haplotype was associated with lower levels of a biomarker of tolerance (Gata3) in males , Here, we confirm this and show the same haplotype is also associated with (i) resistance to Babesia infection, and (ii) inflammatory mode. Although our analyses were not able to identify causal relationships, the most likely scenario is that females with the GC haplotype have a higher risk of infection with Babesia (i.e. they are less resistant) as a result of their more antiinflammatory phenotype (it is difficult to envisage a scenario in which a higher risk of Babesia infection would result in a more anti-, rather than pro-, inflammatory response in these females). Although we did not find any association between the GC haplotype in males and infection risk among the panel of parasites that we measured, our panel was far from exhaustive. It is possible though, that males with the GC haplotype benefit from being more resistant to some unmeasured parasite, as a result of their more pro-inflammatory phenotype.
If so, males with the haplotype would be less tolerant and more resistant, consistent with evidence for a negative genetic correlation (or trade-off) between these two defence strategies (Best, White, & Boots, 2008;Miller, White, & Boots, 2006;Raberg, Sim, & Read, 2007;Restif & Koella, 2003) The parasite in question could be currently circulating in the population, or a recurrent parasite infection (currently absent from the population, or present at very low prevalence). Microtine rodents are known to exhibit cyclic population dynamics, and evidence suggests that parasites, like this, may play an important role in causing such host population cycles (Hudson, Dobson, and Newborn 1998; though see Lambin et al. 1999).
More interestingly perhaps, our results suggest opposing effects of the same polymorphism on the defence strategy of males and females. In order to fully understand this sex by genotype interaction, it is important to consider the fundamental differences in the immune systems of males and females. While males tend to be biased towards the Th1 immune response, females are Th2-biased (Girón-González et al., 2000), particularly during pregnancy, when females suppress their Th1 immune response further, inducing maternal tolerance and protecting the foetus from attack by the maternal immune response (Wegmann, Lin, Guilbert, & Mosmann, 1993). The minimum age of first breeding in M. agrestis is six weeks (Begon et al., 2009;Burthe et al., 2010). The gestation period is around three weeks (Corbet, 1966) and the inter-litter interval is just 20 days (Alibhai & Gipps, 1991;Corbet & Harris, 1991;Morgan Ernest, 2003). Therefore, females in our study population are likely to be pregnant for a large proportion of their lives. Given their depressed Th1 immune response, their defence (and resulting susceptibility) against microparasites may be more sensitive to any shifts in immune architecture. This could explain why we only see impacts of this polymorphism on microparasite infection risk in females. Previous studies have also highlighted the important role of species interactions in the parasite community in driving infection risk in this study population, warning against considering parasites in isolation . It follows then, that although we observed that females with the GC haplotype were more likely to be infected with the parasite, Babesia microti, this particular parasite likely represents a community of parasites, to which this haplotype alters susceptibility in females.
Although the IgE-mediated immune response is known to be important in defence against parasitic worms (Gounni et al., 1994), we found no effect of this IgE receptor polymorphism on risk of worm infection in males or in females. However, IL33, a key pro-inflammatory cytokine involved in the Th2 immune response against macroparasites, was significantly upregulated in males with the GC haplotype. IL33 is known to induce the release of proinflammatory mediators by mast cells, which themselves are activated by Fcer1-bound IgE (Daeron & Nimmerjahn, 2014). Consistent with our previous work, this suggests that in males, the polymorphism alters the strength of the Th2 immune response . However, as is the case for IgE, IL33 has also been implicated in allergy (Cayrol & Girard, 2014). It is possible then, that rather than causing a functional change in the male immune response, this polymorphism is causing an allergic immunopathology.
Our results suggest the same polymorphism is having the opposite effect on the inflammatory response of females. Indeed, suppressor of cytokine signaling-3 (SOCS3) was found to be the most upregulated gene in females with the GC haplotype. SOCS3, as the name suggests, plays an important role in the inhibition of cytokine signal transduction (Carow & Rottenberg, 2014). Given the naturally more immunosuppressed state of females, this suggests that the polymorphism in question may be exaggerating existing sex-differences in immune response. If the polymorphism is (more generally) a risk factor for allergic immunopathology, it is possible that females are better able to suppress this pathology as a result of their existing immunosuppressive architecture, making them more resilient to change (at the cost of an increased susceptibility to Babesia). These findings could therefore shed light on the evolutionary significance of allergic immunity, and could be used to inform animal models of allergic disease. In fact, the Fcer1-binding site of IgE has already been the target for therapies to prevent hypersensitivity (Djukanović et al., 2004). If individual vary widely in their susceptibility to allergy, as a result of their genotype and sex, a personalised medicine approach may be most appropriate.
Effects of the polymorphism on reproductive success were visible in males only. Both reproducing and mounting a pro-inflammatory response are energetically costly activities (Sheldon & Verhulst, 1996). One possibility then, is that because males with the GC haplotype invest more energy in inflammation, they have less available to invest in reproduction. Inflammation and oxidative stress are also closely interrelated. Persistent oxidative stress can lead to chronic inflammation and, consequently, to chronic diseases such as cancer and diabetes (Reuter, Gupta, Chaturvedi, & Aggarwal, 2011). Inflammatory cells including eosinophils, neutrophils, and macrophages are also known to generate reactive oxygen species at sites of inflammation, in turn, leading to oxidative stress (Collins, 1999).
The increased concentration of SOD1 (a biomarker for oxidative stress) in males with the GC haplotype could therefore be the cause, or effect, of a more pro-inflammatory phenotype.
Either way, vertebrate spermatozoa are particularly susceptible to damage by oxidative stress.
As a result, oxidative stress can lead to reduced sperm motility and sperm velocity (e.g. Losdat et al. 2011). Sperm competition is an important feature in the reproduction of Microtine voles, making it particularly important for a male not only to produce many sperm but sperm of high quality, in order to successfully outcompete the ejaculate of other males.
Another possibility then, is that the higher levels of oxidative stress experiences by males with the GC haplotype is reducing the quality of their sperm and, consequently, their fertility. Despite the immune system playing an important role in determining the health and fitness of an individual, natural selection has not converged on a single immune optimum. Instead, individuals in natural populations vary widely in their response to infection. Graham (2013) suggests that there are three different 'scales of enquiry' for understanding this immunoheterogeneity: (1) how, mechanistically individuals vary in their immune response, (2) what, proximately, determines this heterogeneity, and (3) why this heterogeneity is expected to be maintained in natural populations. In this study we have so far addressed the first two scales of enquiry. However, understanding why we see this immunoheterogeneity is a key question in eco-immunology, with potentially important implications for biomedicine (Graham 2013). Sexually antagonistic selection, where a mutation is beneficial to one sex and harmful to the other, is thought to be one mechanism by which balancing selection is generated, and genetic variation in immunity is maintained within natural populations (Graham et al. 2010). Based on samples collected between 2008-2010, we previously identified the three haplotypes at this locus, GC, AC and AT, at frequencies of 0.12, 0.76 and 0.07 respectively . In the present study, conducted between 2015-2017, these frequencies have remained relatively unchanged (0.08, 0.81 and 0.10). The fact that the GC haplotype remains in the population, despite its detrimental effects on male reproductive success, suggests that it may be under balancing selection. Our future work will investigate this possibility further.
In summary, we demonstrate the potential for small genetic differences to have large knockon effects for both the health and fitness of individuals in natural populations. We also demonstrate the importance of considering sex by genotype interactions in this context, with the potential for the same genetic difference to have profoundly different effects for the health and fitness of males and females in natural populations.