Early-life immune expression profiles predict later life health and fitness in a wild rodent

Individuals differ in the nature of the immune responses they produce, affecting disease susceptibility and ultimately health and fitness. These differences have been hypothesised to have an origin in events experienced early in life that then affect trajectories of immune development and responsiveness. Here we investigate early life influences on immune expression profiles using a natural population of field voles, Microtus agrestis, in which we are able to monitor variation between and within individuals though time by repeat (longitudinal) sampling of individually marked animals. We analysed the co-expression of 20 immune genes in early life to create a correlational network consisting of three main clusters, one of which (containing Gata3, Il10 and Il17) was associated with later life reproductive success and susceptibility to chronic bacterial (Bartonella) infection. More detailed analyses supported associations between early life expression of Il17 and reproductive success later in life, and of increased Il10 expression early in life and later infection with Bartonella. We also found significant association between an Il17 genotype and the early life expression of Il10. Our results demonstrate that immune expression profiles can be manifested during early life with effects that persist through adulthood and that shape the variability among individuals in susceptibility to infection and fitness widely seen in natural populations.

Introduction 51 52 Differences at birth, or experienced during neonatal or juvenile development, may 53 impact an individual's ability as an adult to respond to environmental allergens or 54 pathogens, or to maintain health and reproduction (measures of which we define 55 below). In medical sciences, the developmental origin of health and disease (DOMD) 56 theory -which proposes that adult diseases can be traced back to childhood -has 57 been highly influential since its inception and has provided insight into the aetiology 58 of chronic inflammatory and allergic disorders in humans, such as atopy, asthma and 59 autoimmune diseases 1,2 . Similarly, the development of the immune system during 60 early life is an important influence on subsequent health and fitness 3 . The immune 61 system is shaped by environmental exposure to allergens and pathogens, which, in 62 turn, affect the immune responses made following subsequent exposure. Coupled 63 with genetic polymorphisms affecting immune function 4 , early life effects associated 64 with the immune system can give rise to substantial variation among individuals in 65 their immune responses. Such environmental effects can act very early in life, even 66 at birth. For example, children born by caesarean section are exposed to different 67 microbes to vaginal delivery. This has been shown to modify their gut microbiota and 68 the production of pro-inflammatory factors, and ultimately to increase their risk of 69 opportunistic infections 5 and chronic inflammatory disorders in later life e.g. asthma, 70 type 1 diabetes and celiac disease 6,7 . 71 72 While effective immune responses are essential to health by protecting individuals 73 against pathogens, immunopathology (such as tissue damage from inflammatory or 74 over-reactive responses) contributes to significant ill-health in chronic inflammatory 75 8 our sample size) we regressed age against eye lens weight for those voles definitely 174 first sampled prior to breeding (for which we had a birth date) that were culled in later 175 life as part of our cross-sectional study (n = 99). The relationship between age and 176 eye lens weight has been previously described 20 . A quassipoisson GLM with 177 quadratic and cubic terms for eye lens weight, and log link, provided a good model fit 178 for our data (r 2 = 0.56; Supplementary Fig. S1). We used this regression to predict 179 the ages (and hence birth dates) for all other culled voles. Of these, 59 voles were 180 inferred to have first been caught prior to breeding, since their birth date was no 181 more than 6 weeks before their date of first capture, of which 28 had a blood sample 182 taken at first capture. Therefore, the final number of voles first sampled prior to 183 breeding (known and inferred) was 223. 184 185

Immune gene expression assays 186
We used SYBR green based Q-PCR to measure the expression levels of a panel of 187 20 immune-associated genes (Table 1)   T-lymphocyte activation antigen CD86 Cd99l2 CD99 antigen-like protein 2 Foxp3 Forkhead box protein P3 Gata3 GATA binding protein 3 Ifng Interferon gamma Irf2 Interferon regulatory factor 2 Irf9 Interferon regulatory factor 9 Tollip Toll-interacting protein Retnlg Resistin-like gamma Tgfb1 Transforming growth factor beta 1 Il4 Interleukin-4 Il10 Interleukin-10 Il17 Interleukin-17 Il1b Interleukin-1 beta Il1rap Interleukin-1 receptor accessory protein Il1rn Interleukin-1 receptor antagonist protein Apobr Apolipoprotein B receptor Orai1 Calcium release-activated calcium channel protein 1 ccn2 (ctgf) Cellular communication network factor 2 200 All primer sets were designed de novo in-house and validated (to confirm specific 201 amplification and 100±10% PCR efficiency under assay conditions). Ywhaz and Actb 202 were employed as endogenous control genes. We extracted RNA from blood 203 conserved in RNAlater using the Mouse RiboPure Blood RNA Isolation Kit 204 (ThermoFisher), according to manufacturer's instructions. RNA extracts were DNAse 205 treated and converted to cDNA using the High-Capacity RNA-to-cDNA™ Kit 206 (ThermoFisher), according to manufacturer's instructions, including reverse 207 transcription negative (RT-) controls for a subsample. SYBR green-based assays 208 were pipetted onto 384 well plates by a robot (Pipetmax, Gilson) using a custom 209 programme and run on a QuantStudio 6-flex Real-Time PCR System 210 (ThermoFisher) at the machine manufacturers default real-time PCR cycling 211 conditions. Reaction size was 10 µl, incorporating 1 µl of template and 212 PrecisionFAST qPCR Master Mix with low ROX and SYBR green (PrimerDesign) 213 and primers at the machine manufacturer's recommended concentrations. We used 214 three standard plate layouts for assaying, each of which contained a fixed set of 215 target gene expression assays and the two endogenous control gene assays (the 216 same sets of animals being assayed on matched triplets of the standard plate 217 layouts). Unknown samples were assayed in duplicate wells and calibrator samples 218 in triplicate wells and no template controls for each gene were included on each 219 plate. Template cDNA (see above) was diluted 1/20 prior to assay. A main calibrator 220 sample (identical on each plate) was created by pooling cDNA from blood samples 221 taken from many different voles from the study site. As Tollip, Il1rap and Irf2 were 222 relatively poorly represented in this main calibrator sample, a synthesised 478 bp 223 gene fragment containing the amplification target for each of these genes was used 224 as an additional calibrator sample (at 10 × 10 5 copies µl -1 ) in these cases. Samples Half of individuals present in our pedigrees (n = 325) were found to have no 281 offspring. We expect the majority of these to be true zeros (representing actual 282 reproductive failure) as we sampled a large proportion of the total population within 283 clear-cuts. We minimised the chance of false zeros by excluding from the pedigree 284 those individuals (e.g., at the periphery of a study grid) for which we recorded no 285 relatives (including offspring) likely because we had not sampled in the right place. 286

Correlation networks 289
We constructed an expression correlation matrix for all immune genes measured 290 prior to breeding, in early life. We used Spearman Rank correlation coefficients (1) in 291 case of any non-linear relationship between genes, and (2) because our expression 292 data were not normally distributed. There was some missing data in our expression 293 dataset (640/4460 values = 14%). In order to maximise the sample size used to 294 generate each correlation coefficient, we used pairwise complete observations. This 295 resulted in 72% of correlations based on at least 75% of the total data. For each 296 correlation coefficient, we randomly permuted the data 1000 times to calculate a p-297 value. We also adjusted all p-values for multiple testing using the Benjamini-298 Hochberg method. We thresholded our correlation network using these corrected p-299 values -only keeping those edges with a significant corrected p-value (p 0.05). 300 Genes were clustered on edge betweenness, using the edge betweenness algorithm 301 in igraph 27 . The edge betweenness score of an edge is a measure of the number of 302 shortest paths that go through it. The algorithm identifies densely connected 303 modules by gradually removing edges with highest edge betweenness scores 28 . 304

305
We then ran an exploratory analysis, repeating the steps above to construct a To test whether these associations are strictly early-life or also present in adult 370 (breeding) samples, we re-ran the same best models for reproductive success and 371 proportion of later life infected with Bartonella on adult blood samples. We 372 maximised the sample size for these analyses by including all individuals sampled as 373 adults, whether or not they were first sampled prior to breeding and therefore had an 374 approximate birth date. For this reason, we were unable to account for birth month 375 adult sample, we selected one adult sample at random from which to measure 381 expression. To ensure that our results were robust to different random draws, we 382 repeated this process 1000 times, re-running the best models on 1000 random 383 draws. We report the number of random draws for which each association was 384 significant (p 0.05) compared to that expected by chance. We also tested for an 385 association between whether or not an individual expressed each immune parameter 386 of interest early in life, and whether or not they expressed it in later life using a 387 Fisher's exact test. We constructed a correlation network for 20 immune genes whose expression was 404 measured in our voles prior to breeding. Following thresholding, one immune gene, 405 Il4, dropping out of the network. Using clustering analysis we identified three main 406 immune clusters in early life: (1) immune cluster 1 (Il1rap, Il1rn, Il1b, Orai1, Retnlg,  407 Apobr, Cd8a, Ifng, Tgfb1, Cd86), (2) immune cluster 2 (Irf2, Irf9, Foxp3, Tollip) and 408 (3) immune cluster 3 (Gata3, Il10, Il17; Fig. 1a).  spp.) were associated with immune cluster 3 (cluster highlighted in green in Fig. 1b). which an individual was born was the only other fixed effect present in the best 465 model for, and was also significantly associated with, Bartonella infection (Fig. 2b). 466 The association between Il10 expression in later life and probability of infection was 467 also significant (n = 993/1000 random samples for which p 0.05; many more than 468 the 50/1000 expected by chance) indicating that this association was not specific to 469 early-life expression. We found no association between whether or not an individual We have used a wild rodent system to provide evidence that, (1) the reproductive 508 success (a measure of fitness), and (2) the response to infection of mature 509 individuals (a measure of health) is foreshadowed by immune expression in the early 510 stages of their lives. We have also found evidence to suggest that these patterns 511 vary between different parts of the immune network. 512 We were able to derive three clusters of genes related to immune function based on 514 co-expression data and it was unsurprising to see genes involved in different 515 pathways, and with mixed pro-and anti-inflammatory functions, included in the same 516 clusters. This is expected given the complexity of immune regulatory networks. For 517 example, some responses might increase at the same time as their negative 518 feedbacks, or responses within one pathway might lead to host states that favour the 519 expression of another pathway. Nonetheless, we saw very clear functional biases, 520 with the largest cluster containing predominantly pro-inflammatory markers such as 521 Il1 signalling cascade members 33 , Ifng 34 , Tgfb1 35 and Retnlg 36 , but also genes with 522 more nuanced functions in the promotion and suppression of inflammation, including 523 Tgfb1 37 itself and Orai1 38 . This "pro-inflammatory" cluster was distinct from a second 524 cluster containing predominantly immune down-regulatory markers including 525 Foxp3 39 , Tollip 40 and Irf2 41 , but also Irf9 42 a driver of interferon effector responses. A 526 third cluster containing Gata3, Il10 and Il17 may give new insights into immune 527 function, given that it contained a counterintuitive set of genes and was, at the same 528 time, the cluster best linked to host health and fitness. We have previously found 529 GATA3 to be associated with immune tolerance phenotypes, and the association 530 with the anti-inflammatory cytokine IL10 provides further support of this finding 43,44 . 531 Il10 expression is, in turn, positively associated with Il17 expression. IL17, which we 532 found to be an early-life indicator of reproductive success, is typically associated with 533 a pro-inflammatory immunopathology but may also be protective against the 534 numerous primary microbial insults that young voles, lacking acquired immunity, may 535 encounter 45 . The correlation between IL10 and IL17 may thus be one whereby IL10 536 is produced in response to IL17 in order to safely regulate its inflammatory effects 537 and prevent immunopathology 46-49 . We argue that wild rodent studies can help us to 538 define which of the possible immunological interactions known from laboratory 539 studies predominate and shape the structure of immune phenotypes in individuals 540 exposed to pathogens, nutritional challenge or other stressors within the natural 541 environment. As in the case of IL10 and IL17, wild rodents may also give insights 542 into immune regulation within a natural context. 543

544
Having identified the correlational structure of immunological traits in a wild 545 population, it is then a natural extension to include other phenotypic traits. We used 546 longitudinal data to address our aim to determine how immune phenotypes may be 547 laid down in early life to affect later life health and fitness. We found two pathways 548 through which immune phenotypes in early life act on later life traits (Fig. 3). In the 549 first, increased expression of Il10 in blood samples taken during early life (i.e., prior 550 to breeding) was subsequently associated with both a ten-fold increased 551 susceptibility to Bartonella infection later in life and a three-fold reduction in 552 survivorship. The association of a genotypic effect on Il10 expression, a 553 polymorphism in Il17, suggests that early-life effects may be, at least in part, laid 554 down at birth in this case. A challenge for future work will be to tease out causation 555 from correlation in this case. An adverse Il17 genotype may make an individual more 556 likely to become infected with Bartonella, with infection stimulating Il10 expression. 557 Alternatively, that genotype may predispose an individual to express Il10 early in life 558 The second pattern that we found is for increased Il17 expression early in life to be 575 associated with a subsequent five-fold increased likelihood of reproduction. By 576 contrast, Il17 expression in breeding adults showed no such association with 577 reproduction, and indeed, no association either between Il17 expression early and 578 later in life. In contrast to Il10 expression, no genetic effects on Il17 expression were 579 detected; at least from the deliberately limited set of cytokine polymorphisms 580 considered, which included Il17 polymorphism itself. This suggests that early life 581 expression of Il17 may be the consequence of an environmental effect leading to an 582 immune phenotype that acts on juveniles and influences their development as they 583 transition into reproductively mature adults. This could be, for example, a bacterial 584 infection or composition of the gut microbiome early in life. least over a much more extended period. Il10 differences confirm a genetic basis for 598 immune expression that, as previously shown in this system and others, act through 599 the network of interacting cytokines such that polymorphism at one cytokine gene 600 can affect expression of other cytokines and immune phenotypes 23,51-56 , with these 601 effects further percolating to affect health and fitness more broadly in later life. By 602 contrast, the Il17 differences appear to be more reminiscent of the Jesuit motto "Give 603 me the boy until he is seven and I will give you the man": that independently of any 604