The magenta module has a modest, but strongly significant correlation to the fraction of cells that are eosinophils (*r* = 0.28, *p* = 0.005). A recent review highlighted the importance of eosinophils in host-helminth interactions. Specifically, at the site of host tissue damage, eosinophils are primed by fibronectin, and produce a variety of proteins which are toxic to helminths. Furthermore, the diversity of eosinophil cell surface receptors makes them central to mediate the inflammation response. Finally, helminths appear to have a number of anti-eosinophil proteins which both evade and dampen host response to the presence of helminths (58).

Constructing separate signed coexpression networks for each population reveals dramatic differences in network structure.

Using the same construction parameters for each population, the Rob network resembles the combined population network, showing strong module-trait correlations for worm size, cell population phenotypes, ROS production, and sex across many different modules. In stark contrast, the Gos network is much more static, with only a subset of traits which were previously significant correlated to a single module. Overall, Gos fish have many fewer total modules and correlations between modules and traits are much weaker. To estimate relationships between networks, we calculated module similarity between every pair of Rob and Gos modules by shared gene membership. Broadly, large modules in Gos are split into many smaller modules in Rob. Thus, we conclude that there are two possible outcomes: (1) the Rob fish have evolved a more modular and dynamic repertoire of expression with which to fend off cestodes, or (2) the Gos fish evolved to maintain constitutive gene expression in response to cestodes, instead adopting a tolerance strategy (See Figures S2–10 in Supplementary Material for additional details on shared and population specific signed networks, and module similarity heat map).

Our WGCNA results confirm observed population differences in immune function and patterns of gene expression in the linear modeling. Furthermore, they support our hypothesis that there are two traits involved in stickleback resistance to cestodes: (1) innate immune response to prevent cestode establishment and (2) limiting worm growth once cestodes become established. These two traits separate into distinct modules of gene expression, each enriched for genes with immunological function matching *a priori* expectations. This two-trait perspective refines the question of variation in cestode prevalence among stickleback populations by focusing attention to both early and late stage infection. Future studies will be needed to describe the full timeseries of gene expression as exposure and infection proceeds in these study populations, and to establish directionality of the interaction between cestodes and stickleback. While others have argued that cestodes are the primary drivers of coevolution (47), only by sequencing both host and parasite mRNAs can we hope to detail this interaction at the molecular level. Host genotype by parasite genotype interactions offers a promising opportunity for further study. Such GxG interaction was recently described in the stickleback-cestode system, but only documented growth phenotypes, and no attempt to describe the genetic basis for such traits has been made (59). We maintain that this type of study provides a means to identify evolutionarily labile genes, which underlie beneficial shifts in immune function across the geographic mosaic of host-parasite coevolution.

### Summary

Using a large scale controlled laboratory infection experiment, we find changes in gene expression between two host populations, and as a function of infection status. For a smaller portion of genes, the expression response to infection differed between the two host populations. These findings are consistent with observations of host immune function in the same infection experiment (29). ROS production and recycling, B cell activation and targeting, and fibrosis appear to play important roles in stickleback defense against cestodes. Our analysis also suggests that host resistance involved two components; response to challenge by cestodes, and control over cestode growth. Furthermore, differences in coexpression networks between populations suggest that either Rob fish have evolved a more elaborate expression profile or that Gos fish are shutting down expression to tolerate cestodes. Our results not only suggest a mechanistic link between host immune phenotypes and candidate genes, but also provide the foundation for studying the direct effects of host alleles on parasite fitness.

### MATERIALS AND METHODS

We obtained mRNA from the head kidneys of stickleback from three experimental groups: unexposed controls (*N* = 16 and 19, Rob and Gos fish, respectively), exposed but ultimately uninfected stickleback (*N* = 21 and 16), or exposed and infected (*N* = 17 and 9) fish. Tissues from the latter two groups were harvested 42 days post exposure. A brief summary is presented below, but see Weber et al. (29) for full experimental methods. We focus on expression in head kidneys as it is the major site of immune cell development in fish (47, 60–62), and head kidney cell cultures were used to measure immune function independently of gene expression (stickleback have two head kidneys).

### Summary of Breeding and Infection As Described in Weber et al. (29)

Weber et al. (29) generated 17 pure Rob and 16 pure Gos families by *in vitro* fertilization of wild caught parents. F1 fish were raised to adulthood. Members of each family were either sham exposed or exposed to cestodes *via* uninfected/infected copepods, respectively. Infection and screening of copepods was done by hand in a laboratory setting. After exposure, fish were grown for 42 days before being euthanized. For the present study, we subsampled from these fish, attempting to balance both population and infection status. We incorporated at least one individual from each treatment per family. Our final sample sizes (TagSeq libraries post outlier removal) are: Gos: Control: 19, Exposed: 21, Infected: 17, across 11 families. Rob: Control: 16, Exposed: 16, Infected: 9 across, 8 families (detailed in Table S1 in Supplementary Material).

### Sample Collection, Sequence Library Construction, and Analysis of Flow Cytometry Data

Briefly, head kidneys were dissected from stickleback; one head kidney was used for flow cytometry (measuring ROS and describing cell populations) and the other was preserved in RNAlater at −20°C. Detailed protocols can be found in Weber et al. (29). From the stored head kidney, RNAseq libraries were constructed according to Lohman et al. (31). This method focuses sequencing effort on short "tags" adjacent to the poly-A tail of mRNAs rather than distributing sampling effort across the entire mRNA as in traditional RNAseq. Samples were sequenced on the Illumina HiSeq 2500 at the Genome Sequence and Analysis Facility at the University of Texas at Austin, producing ~6.7M raw reads per sample. See Weber et al. (29) for full ROS and flow cytometry methods. Flow cytometry data were analyzed using FlowJo software (Treestar). Granulocyte and lymphocyte populations were defined based linear forward scatter (FSC) and side scatter (SSC) gating described in Weber et al. (29). Precursor, myeloid, lymphoid, and eosinophil populations were defined using linear FSC and logarithmic SSC gating as described in Wittamer et al. (63).

### Bioinformatics

TagSeq reads were processed (removal of PCR duplicates, adapter contamination, and average base quality score –*n* 20) according to the iRNAseq pipeline (64) using version 79 of the stickleback transcriptome from Ensemble. Resulting genes were filtered for mean counts greater than 1 among all samples, producing 9077 genes among all samples. Transcriptome annotations were based on the UniProtKB database (http://www.uniprot.org/ help/uniprotkb) and followed previously described procedures (64). Code for the iRNAseq pipeline can be found here: https:// github.com/z0on/tag-based\_RNAseq. Code for the annotation pipeline can be found here: https://github.com/z0on/annotating Transcriptomes.

## Statistical Analysis with DESeq2

We scanned for outliers using arrayQualityMetrics (65) and removed one sample because of insufficient read depth (final *N* = 98). To test for differential gene expression, we used the following model in DESeq2:

```
Yij ∼ βBatch P +β +β +β opulation InfectionStatus Population*InfectionStatus + εij (1)
```
where *Yij* is the count of gene *i* in individual *j*, βPopulation is a fixed effect with two levels: Rob and Gos, βInfectionStatus is a fixed effect with three levels: control, exposed (but not infected), and infected, and full-sibling families are nested within populations. βBatch is the lane on which samples were sequenced. An additional predictor βSex was included for genes when appropriate (lower AIC score) and improved the model fit of 839 genes total. We fit the full model (including sex) to all genes and then extracted only the 839 that were improved by the addition of sex and looked for significant *p* values for main effects and interactions. With the full model, 67 of these "sex improved" genes were significantly different between populations. No genes were significant for either exposure or infection, and one gene was significant for the interaction of population and infection status (myosin 5ab, ENSGACG00000006025: log2foldchange = 2.98, Wald *p* = 0.07). All *p*-values were multiple test corrected using 10% FDR (Benjamini–Hochberg). Although fish from the controlled infection experiment were exposed to three different parasite genotypes (each family exposed to only one parasite genotype), we are only interested in the host response to any parasite, and therefore average across parasite genotypes by simply not including this as a term in our linear model.

### GO with GO\_MWU

We used the Mann–Whitney *U* test for GO analysis. This approach has been described (66) and the code for analysis can be found here: https://github.com/z0on/GO\_MWU.

### Weighted Coexpression Gene Network Analysis

Raw read counts were normalized using limma (67) for input into WGCNA (57). All genes included in the post filter data were included in the WGCNA analysis. Following the walkthrough in Langfelder and Horvath (57), we built a signed network with a soft thresholding power of 7, and a minimum module size of 30 genes. Following dynamic tree cut, we merged modules with greater than 80% similarity, producing 14 modules. We separated Rob and Gos samples and repeated this process with the same parameters.

### Caveats and Limitations

Our study flips the traditional search for immune candidate genes from inbred lab strains to wild populations, using historical natural selection as a tool to screen for changes in gene expression associated with parasite infection. While our host-parasite model system is powerful, it does have some limitations. The reference genome is of generally good quality but annotation is lacking (approximately 22.5% of the entries in the stickleback genome are either unnamed or labeled as novel genes). Thus, GO analysis is performed after assigning GO accession terms by BLAST homology, rather than functional verification, a common solution for non-model systems. The features of the stickleback genome may be missing potentially interesting immunological genes which are sufficiently diverged from human or mouse genes and therefore may be unannotated. In particular, the number and location of MHC II paralogs remains uncertain, illustrating need for genome sequence improvement. With respect to MHC II in particular, the unknown number of gene copies in the populations of interest may potentially complicate our estimates of gene expression.

Our linear modeling with DESeq2 employs appropriate FDR correction, but we choose to accept higher than "standard" *p*-values associated with LFC because of the direct connection between candidate genes and independently observed immune phenotypes. If, for example, we had not measured immune phenotypes, we would not accept log2 fold changes in expression with associated *p*-values greater than 0.05 but less than 0.1. Furthermore, we chose a very low base min mean filter because we have high confidence in detecting lowly expressed genes (31).

### REFERENCES

We also wished to include more genes in our enrichment analysis and this also detracts from our power due to multiple test correction. We accordingly accept slightly larger than normally allowed *p*-values. Our TagSeq based approach has been shown to be at least as good as total RNAseq methods (having an equal or higher correlation between observed and known values of a spike in control) but does not account for splice variants or copy number variation, which may be potentially important in the evolution of immune responses.

Our study used tissue from a single organ (head kidneys) for both gene expression and immune phenotype measures. Head kidneys are a crucial hematopoetic organ in fish, but analysis of other tissues may produce different results. Moreover, head kidneys contain multiple immune cell populations that we are unable to sort effectively for cell-type-specific expression studies. We do use cell population counts (proportion granulocytes vs. lymphocytes) as a covariate, which as noted for MHC weakly contributes to expression variation of a few genes. But, lacking monoclonal antibodies to many immune cell receptors in stickleback, we cannot readily distinguish among finer subdivisions of cell types. This resource limitation, typical of most non-model organisms, limits our ability to statically detect effects of cell population composition on expression.

### DATA ARCHIVAL LOCATION

Instructions for accessing the code for analysis and plotting, and raw sequence data from UT Austin's permanent data archive are located in supplementary materials.

### AUTHOR CONTRIBUTIONS

BL built, sequenced, and analyzed TagSeq libraries. NS and JW provided samples, flow cytometry data, and comments on analysis and interpretation of results. BL and DB wrote the manuscript with comments from NS and JW. All authors approved the final version.

### ACKNOWLEDGMENTS

This work was supported by the Howard Hughes Medical Institute (DB). We thank John Lovell and Marie Strader for helpful comments during data analysis and writing.

### SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at http://journal.frontiersin.org/article/10.3389/fimmu.2017.01071/ full#supplementary-material.

immunity at large. *Annu Rev Immunol* (2006) 24:353–89. doi:10.1146/ annurev.immunol.24.021605.090552

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

*Copyright © 2017 Lohman, Steinel, Weber and Bolnick. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*

# The Kinetics of the humoral and interferon-gamma immune responses to experimental *Mycobacterium bovis* infection in the White rhinoceros (*Ceratotherium simum*)

### *Edited by:*

*Andrew Steven Flies, University of Tasmania, Australia*

#### *Reviewed by:*

*Francisco Javier Salguero, University of Surry, United Kingdom David E. MacHugh, University College Dublin, Ireland Maria A. Risalde, Instituto Maimonides de Investigación Biomédica de Cordoba (IMIBIC), Spain*

> *\*Correspondence: Sven D. C. Parsons<EMAIL_ADDRESS>
#### *†Present address:*

*Ross McFadyen, Centre for Proteomic and Genomic Research (CPGR), Cape Town, South Africa*

#### *Specialty section:*

*This article was submitted to Comparative Immunology, a section of the journal Frontiers in Immunology*

> *Received: 18 June 2017 Accepted: 04 December 2017 Published: 22 December 2017*

#### *Citation:*

*Parsons SDC, Morar-Leather D, Buss P, Hofmeyr J, McFadyen R, Rutten VPMG, van Helden PD, Miller MA and Michel AL (2017) The Kinetics of the Humoral and Interferon-Gamma Immune Responses to Experimental Mycobacterium bovis Infection in the White Rhinoceros (Ceratotherium simum). Front. Immunol. 8:1831. doi: 10.3389/fimmu.2017.01831*

*Sven D. C. Parsons1,2,3\*, Darshana Morar-Leather4 , Peter Buss5 , Jennifer Hofmeyr5 , Ross McFadyen1,2,3,6†, Victor P. M. G. Rutten4,7, Paul D. van Helden1,2,3, Michele A. Miller1,2,3 and Anita Luise Michel4,8*

*1DST/NRF Centre of Excellence for Biomedical Tuberculosis Research, Cape Town, South Africa, 2SAMRC Centre for TB Research, Cape Town, South Africa, 3Division of Molecular Biology and Human Genetics, Faculty of Medicine and Health Sciences, Stellenbosch University, Cape Town, South Africa, 4Bovine Tuberculosis and Brucellosis Research Programme, Faculty of Veterinary Science, Department of Veterinary Tropical Diseases, University of Pretoria, Onderstepoort, South Africa, 5Veterinary Wildlife Services, South African National Parks, Kruger National Park, Skukuza, South Africa, 6Centre for Proteomic and Genomic Research (CPGR), Cape Town, South Africa, 7 Faculty of Veterinary Medicine, Department of Infectious Diseases and Immunology, Utrecht University, Utrecht, Netherlands, 8National Zoological Gardens of South Africa, Pretoria, South Africa*

*Mycobacterium bovis* is the cause of tuberculosis (TB) in a wide range of species, including white rhinoceroses (*Ceratotherium simum*). Control of the disease relies on the indirect detection of infection by measuring pathogen-specific responses of the host. These are poorly described in the white rhinoceros and this study aimed to characterize the kinetics of immune responses to *M. bovis* infection in this species. Three white rhinoceroses were infected with *M. bovis* and their immune sensitization to this pathogen was measured monthly for 20 months. Cell-mediated immunity was characterized in whole blood samples as the differential release of interferon-gamma in response to bovine purified protein derivative (PPDb) and avian PPD (PPDa) as well as the release of this cytokine in response to the *M. bovis* proteins 6 kDa early secretory antigenic target (ESAT-6)/10 kDa culture filtrate protein (CFP-10). Humoral immunity was quantified as the occurrence or the magnitude of antibody responses to the proteins ESAT-6/CFP-10, MPB83, MPB83/MPB70, and PPDb. The magnitude and duration of immune reactivity varied between individuals; however, peak responses to these antigens were detected in all animals circa 5–9 months postinfection. Hereafter, they gradually declined to low or undetectable levels. This pattern was associated with limited TB-like pathology at postmortem examination and appeared to reflect the control of *M. bovis* infection following the development of the adaptive immune response. Measurement of these markers could prove useful for assessing the disease status or treatment of naturally infected animals. Moreover, immune responses identified in this study might be used to detect infection; however, further studies are required to confirm their diagnostic utility.

Keywords: immune response, interferon-gamma, serology, *Mycobacterium bovis*, tuberculosis, white rhinoceros

## INTRODUCTION

White rhinoceroses (*Ceratotherium simum*) are classified as "Near Threatened" by the International Union for Conservation of Nature, with the majority of animals occurring in South Africa (1). Of these, a substantial number occur in the greater Kruger National Park (KNP) and the Hluhluwe-iMfolozi Park (HiP). However, because of their economic value and threatened conservation status, animals from these populations are regularly translocated to other reserves and privately owned collections. Importantly, movement of animals from these areas could present a risk of translocating *Mycobacterium bovis*, a major cause of tuberculosis (TB). This pathogen can infect a wide variety of domestic and wildlife hosts, including rhinoceros species (2) and has become established in both the KNP and HiP (3). The recent detection, in the KNP, of a case of severe pulmonary TB in a black rhinoceros (4) and confirmed *M. bovis* granulomas in lymph nodes of four white rhinoceroses (unpublished data) highlights the potential risk of movement of these species.

Tuberculosis is slowly progressive and the causative organisms may initially be contained within well circumscribed granulomas (5). For this reason, detection of the pathogen can be challenging and, as in other species, TB in rhinoceroses might only be diagnosed postmortem or once animals have developed advanced disease (2). Infection is therefore often diagnosed indirectly by measuring the host's adaptive immune response toward *M. bovis* antigens. This is commonly done by quantifying either the *in vivo* or the *in vitro* immune response to purified protein derivative (PPD), a preparation containing a broad range of *M. bovis* antigens (6, 7). Alternatively, recombinant proteins that are more specific to *M. bovis* can be utilized as test antigens and these include 6 kDa early secretory antigenic target (ESAT-6), 10 kDa culture filtrate protein (CFP-10), MPB70, and MPB83 (8). Assessment of immune responses to these and other antigens might also be used to distinguish between latent infection and progressive disease or to monitor treatment in both humans and animals (8, 9).

The present study forms part of a broader project that characterized the clinical features and associated gross and histopathology of the experimental infection of three white rhinoceroses with *M. bovis* (10). In all cases, 20 months after infection, animals had shown no clinical signs of TB disease and had limited TB-like pathology (10). We hypothesized that quantifying the adaptive immune responses to *M. bovis* in these animals would provide an indirect measure of their infection and disease status. As such, we aimed to characterize (i) their immune sensitization to selected antigens following infection and (ii) the kinetics of their humoral and cell-mediated immune responses.

### MATERIALS AND METHODS

### Animals

The capture, maintenance, chemical immobilization, infection, and sampling of animals in the KNP, as well as biohazard containment, have previously been described in detail (10). Briefly, three subadult male white rhinoceroses, identified as PB1, PB2, and PB4, were infected by endoscopic endobronchial instillation of the *M. bovis* strain SB0121, a genotype commonly isolated from wildlife in the KNP. The inocula for the three animals contained approximately 2.1 × 103 colony forming units (cfu), 1.8 × 102 and 1.4 × 103 cfu, respectively. Each month, from 3 months prior to infection until 20 months postinfection (PI), animals were chemically immobilized and blood was collected from the radial vein into lithium heparin and serum vacutainer tubes (Fisher Scientific, Suwanee, GA, USA) and used in the immunological assays described below (**Table 1**). On these occasions, endoscopic bronchoalveolar lavages were performed for mycobacterial culture. Twenty months after infection, animals were euthanized, postmortem examinations performed, and systematic tissue sampling conducted to determine the presence or absence of *M. bovis* by histopathology, mycobacterial culture, and polymerase chain reaction (PCR) and findings have been previously reported in detail (10). Approval for the study was obtained from the Animal Ethics Committees of the South African National Parks and Stellenbosch University (proposal SU-ACUM12-00012) as well the South African National Department of Agriculture, Forestry and Fisheries in terms of Section 20 of the Animal Diseases Act (Permit 12/11/1/1/6/1).

## Rhinoceros-Specific Interferon-Gamma (IFN-**γ**) Release Assay

Blood collected in heparin tubes was processed within 12 h after collection as previously described (11). Briefly, whole blood samples were incubated at 37°C in 5% CO2 for 24 h with *M. bovis* PPD (PPDb, 20 µg/ml) and *M. avium* PPD (PPDa, 20 µg/ml), PMA/CaI (0.1/2 μg/ml—positive control), culture medium (Nil—negative control), respectively. Plasma was harvested following centrifugation at 1,088 × *g* for 5 min and stored at −80°C until tested, in duplicate, in the rhinoceros-specific IFN-γ capture enzyme-linked immunoassay (ELISA) (11). The results of the ELISA were determined at 490 nm using an ELISA plate reader (BioTek, Powerwave XS2, Gen5 software). An initial reference/ blank reading was performed on the plate prior to the blocking step, at the same wavelength, and final optical density (OD) values were determined by subtracting the mean reference value from the mean test value for each sample well. Antigen- and mitogenspecific release of IFN-γ was calculated as the OD value derived from the PPD- and PMA/CaI-stimulated samples minus that derived from the Nil sample. *M. bovis*-specific release of IFN-γ was calculated as the PPDb value minus the PPDa value (ΔPPDb-a).

### Modified QuantiFERON TB Gold (In-Tube) (QFT) Assay

One milliliter of heparinized whole blood was transferred to a "Nil" tube (containing saline) and a "TB Antigen" tube (coated

**Abbreviations:** CFP-10, 10 kDa culture filtrate protein; ELISA, enzyme-linked immunoassay; ESAT-6, 6 kDa early secretory antigenic target; HiP, HluhluweiMfolozi Park; IFN-γ, interferon-gamma; KNP, Kruger National Park; OD, optical density; PBS, phosphate-buffered saline; PI, postinfection; PPDa, avian purified protein derivative; PPDb, bovine purified protein derivative; QFT, modified QuantiFERON TB Gold (In-Tube) assay; TMB, 3, 3′, 5, 5′- tetramethylbenzidine.

TABLE 1 | Immunoassays utilized to measure immune sensitization to selected antigens in *Mycobacterium bovis*-infected white rhinoceroses.

*PPD, purified protein derivative; IGRA, interferon-gamma release assay; QFT, QuantiFERON TB Gold (In Tube); N/A, not applicable; ELISA, enzyme-linked immunosorbent assay; LFD, lateral flow device.*

with peptides simulating ESAT-6, CFP-10, and TB 7.7) of the QFT system (Qiagen, Hilden, Germany). In addition, as a measure of cell viability, 1 ml blood was incubated with phytohemagglutinin (PHA) (Sigma-Aldrich, St. Louis, MO, USA) in phosphatebuffered saline (PBS) at a final concentration of 10 µg/ml. The tubes were shaken according to the manufacturer's instructions and incubated for 20–24 h at 37°C. Hereafter, the tubes were centrifuged at 1,600 × *g* for 10 min and plasma was harvested and stored at −80°C. Plasma samples were assayed in duplicate using a commercial bovine IFN-γ ELISA cross-reactive with IFN-γ of sheep and horses (Kit 3115-1H-20; Mabtech AB, Nacka Strand, Sweden) that has previously been shown to detect recombinant rhinoceros IFN-γ (data not shown). Reactions were visualized using 3,3′,5,5′-tetramethylbenzidine (TMB) (BD Biosciences, NJ, USA) as a color substrate. The IFN-γ concentration in each sample was measured as the OD of each well, at a wavelength of 450 nm, using a Labtech LT-4000 microplate reader (Lasec, Cape Town, South Africa). The *M. bovis*-specific release of IFN-γ was calculated as the mean OD derived for plasma harvested from the TB Ag tube minus the mean OD derived for plasma harvested from the Nil tube.

### PPD ELISA

Flat-bottomed 96-well microtiter plates (Nunc, New York, NY, USA) were coated with 100 µl of a 10 µg/ml PPDb solution (Prionics, Schlieren-Zurich, Switzerland) in 0.05 M carbonatebicarbonate buffer (pH 9.6) and incubated overnight at 4°C. For each well coated with antigen, a corresponding control well was coated with 100 µl blocking buffer (BB) consisting of 5% milk powder (Clover, Roodepoort, South Africa) in PBS + 0.05% Tween-20 (Sigma-Aldrich, St. Louis, MO, USA). After incubation, plates were decanted and washed five times with PBS containing 0.05% Tween-20, then blocked with 200 μl/well of BB for 1 h at room temperature (RT). The plates were washed as above and serum samples diluted 1:200 in BB were added to duplicate wells (100 μl/well). Plates were incubated at RT for 1 h and then washed five times as above. Plates were then incubated with 100 μl/well of peroxidase-conjugated recombinant protein A/G (Thermo Scientific, MA, USA) diluted 1:100,000 in PBS, for 1 h at RT. Plates were washed as above before addition of 100 μl/well of TMB and subsequently incubated in the dark for 15 min. The reaction was stopped using 50 μl/well of 2 M H2SO4 and the OD of each well was measured at 450 nm using a LT-4000 Microplate Reader (Lasec). For each animal, an assay result was calculated as the mean OD value derived from duplicate PPD-coated wells minus that of duplicate BB-coated wells.

### ElephantTB STAT-PAK**®** Assay

Rhinoceros sera were tested using the ElephantTB STAT-PAK® assay (Chembio Diagnostic Systems, Inc., Medford, NY, USA). The assay has been optimized for the detection of *M. tuberculosis* infection in elephants but is not specific for a particular host species and has previously been found useful in measuring antibody responses in a black rhinoceros infected with *M. tuberculosis* (12). Briefly, 30 µl of serum was added to the sample well followed by three drops of diluent/antibody detection conjugate. If the sample contained antibodies to *M. bovis* antigens (ESAT-6/ CFP-10/MPB83), a positive line appeared as a blue band. Any visible band observed in the test line area by two independent observers and read at 20 min was considered as an antibody positive result.

### Dual Path Platform (DPP)**®** VetTB and Bovid DPP Assays

The DPP® VetTB and Bovid DPP assays (Chembio Diagnostic Systems, Inc.) were performed according to the manufacturer's instructions as previously described (12). After 15 min, the test results were read. A DPP optical reader device (Chembio Diagnostic Systems, Inc.) was used to measure the reflectance of test strips and a result was quantified as a numerical score, represented as reflectance units (RU). The DPP® VetTB assay includes two test lines containing the antigens MPB83 and ESAT-6/CFP-10, along with a positive control. In the Bovid DPP assay, a single test line contained MPB83/MPB70, along with a positive control.

### RESULTS

## Rhinoceros-Specific IFN-**γ** Assay

The release of IFN-γ in unstimulated blood was negligible, resulting in median ELISA OD values of 0.08, 0.05, and 0.07 for PB1, PB2, and PB4, respectively. All animals displayed strong IFN-γ responses to PMA/CaI stimulation, with mean mitogen-specific OD values of 0.86, 0.83, and 0.86, respectively. The kinetics of PPDa and PPDb-specific IFN-γ release in whole blood from these animals are illustrated in **Figure 1A**. Immune sensitization to these antigens was negligible prior to infection and first observed at 2–3 months PI. The release of IFN-γ in response to both antigens was greatest between 4 and 10 months PI and the *M. bovis*-specific

FIGURE 1 | The kinetics of the cell-mediated immune response to *Mycobacterium bovis* infection in the white rhinoceros. Following *M. bovis* infection, whole blood from three animals (PB1, PB2, and PB4) was incubated overnight at 37°C without stimulation, or stimulation with bovine PPD purified protein derivative (PPDb), avian PPD (PPDa), and 6 kDa early secretory antigenic target (ESAT-6)/10 kDa culture filtrate protein (CFP-10)/TB7.7 peptides, respectively. Plasma interferon-gamma (IFN-γ) was measured by enzyme-linked immunoassay, in duplicate, as the optical density (OD), and antigen-specific release calculated as the OD derived for the unstimulated sample subtracted from that derived for each antigen (mean ΔOD ± SD). (A) Antigen-specific IFN-γ release in response to PPDb and PPDa, and the differential response to these antigens (ΔPPDb-a); (B) Antigen-specific IFN-γ release in response to ESAT-6/CFP-10/TB7.7 peptides.

PPD responses (ΔPPDb-a) peaked at circa 5–6 months PI in PB1 and PB2 and at 11 months PI in PB4. Subsequently, IFN-γ release in responses to both PPDa and PPDb were low until 20 months PI when all animals again showed substantial responses to these antigens. At this time, IFN-γ release was notably greater in response to PPDa than PPDb in PB2 and PB4, while in PB1, the PPDb response was higher than that to PPDa.

### Modified QFT Assay

The release of IFN-γ in unstimulated blood was negligible, resulting in median ELISA OD values of 0.09, 0.08, and 0.09 for PB1, PB2, and PB3, respectively. As a result of experimental error, measurements of PHA-induced IFN-γ responses were available for only 14/21 sampling occasions for each animal. On these occasions, median mitogen-specific OD values were 0.81, 0.35, and 0.89, respectively. The kinetics of ESAT-6/CFP-10/TB7.7 specific IFN-γ release in whole blood from PB1, PB2, and PB4 are illustrated in **Figure 1B**. These antigen-specific responses were first observed at 1 month PI in PB2 and at 2 months PI in PB1 and PB4. Following infection, peak responses occurred at circa 6–9 months PI, after which IFN-γ release decreased gradually over time. At 20 months PI, PB2 displayed a moderate but distinct increase in IFN-γ release in response to ESAT-6/ CFP-10/TB7.7.

## PPD ELISA, ElephantTB STAT-PAK**®**, and DPP Assays

The kinetics of the humoral responses to *M. bovis* antigens in PB1, PB2, and PB4 are illustrated in **Figure 2**. Using the ElephantTB STAT-PAK® assay, PB1 displayed humoral sensitization to the combined antigens ESAT-6/CFP-10/MPB83 at all time points from 1 to 10 months PI (**Figure 2A**). In contrast, antibodies to this antigen pool were not detected in PB2 for the duration of the study and only on a single occasion in PB4, i.e., 2 months PI (**Figure 2A**). Using the DPP assays, only PB1 displayed a strong and sustained humoral response to the *M. bovis* antigens ESAT-6/CFP-10, MPB83, and MPB83/MPB70, with sensitization to all three antigens detected at 5 and 6 months PI. Hereafter, the humoral response to MBP83 alone was not again detected in this animal; however, responses to MBP83/MPB70 were sustained until 10 months PI and responses to ESAT-6/CFP-10 peaked again at 14 months PI. Both PB2 and PB4 displayed very low antibody quantities specific for these antigens, with PB4 showing only a moderate response to MPB83/MPB70 at 17 months PI.

Prior to infection, levels of circulating antibodies to PPDb (**Figure 2B**) were similar in all animals, with levels rising from 2 to 4 months PI. A peak in the PPDb-specific humoral response occurred circa 6 months PI, after which this gradually decreased in all animals.

FIGURE 2 | The kinetics of the humoral response to *Mycobacterium bovis* infection in the white rhinoceros. Following *M. bovis* infection, serum antibodies of three animals (PB1, PB2, and PB4) were measured using selected commercial assays and an in-house bovine purified protein derivative (PPDb) enzyme-linked immunoassay (ELISA). (A) Serological responses to MPB83 and 6 kDa early secretory antigenic target (ESAT-6)/10 kDa culture filtrate protein (CFP-10) (Dual Path Platform (DPP)® VetTB assay) and MPB83/70 (Bovid DPP assay); seroreactivity to ESAT-6/CFP-10/MPB83 (ElephantTB STAT-PAK® assay) is indicated (\*). (B) Serological responses to PPDb, measured in duplicate by indirect ELISA, and calculated as the difference in optical density (OD) between control wells and PPDb-coated wells (mean ΔOD ± SD).

### DISCUSSION

Following experimental endobronchial instillation of *M. bovis* in 3 white rhinoceroses, these developed time-dependent humoral and cell-mediated immune responses toward the antigens PPDa, PPDb, ESAT-6/CFP-10, MPB83, MPB83/MPB70, and ESAT-6/ CFP-10/MPB83. Notably, cell-mediated responses toward ESAT-6/CFP-10 and PPD showed greater consistency between individuals than did humoral responses toward MPB83, MPB83/ MPB70, and ESAT-6/CFP-10/MPB83. While the specific patterns of immune reactivity varied between individuals, peak responses were detected in all three animals circa 5–9 months PI, after which they gradually declined to low or undetectable levels.

In mice, the onset of adaptive immunity toward *M. tuberculosis* requires the active replication of the pathogen in lymph nodes draining the site of infection (13). Moreover, the progression and magnitude of this response is associated with the amount of antigen produced in these lymph nodes (13, 14). Similarly, *M. bovis* infection in cattle and deer invariably results in the development of granulomatous disease in lymph nodes draining the site of infection, suggesting that these tissues are the primary site of antigen presentation and immune activation (15). In the present study, during the first 5–9 months PI, the rhinoceroses displayed a consistent rise in immune reactivity toward *M. bovis* antigens, strongly suggesting the establishment of this infection in all three animals.

The magnitude and duration of *M. bovis*-specific immune reactivity appeared to reflect the disease status of each animal. PB1, which displayed the most robust and sustained responses to all test antigens was the only animal from which viable *M. bovis* was isolated, i.e., from a bronchoalveolar lavage sample collected 5 months following infection (10). Moreover, at postmortem examination of PB1, acid fast bacterial rods and *M. bovis* DNA were detected in tracheobronchial and lung lesions, respectively (10). In contrast, PB4 displayed substantially lower responses of both humoral and cell-mediated immunity (CMI) and no definitive confirmation of *M. bovis* infection was made in this animal by culture, histopathology, or PCR (10).

Following an initial peak circa 5–9 months PI, all rhinoceroses displayed a gradual decline in the magnitude of their immune responses, and antigen-induced IFN-γ release reached low levels within 12–16 months after infection. In cattle, the severity of TB disease correlates with both the magnitude of ESAT-6-specific CMI (5) and the humoral response to *M. bovis* (6). Similarly, in rabbits experimentally infected with *M. tuberculosis*, progressive disease is associated with strongly rising titers of anti-PPD antibodies (16), while control of the infection is characterized by low seroreactivity to PPD that declines over time (17). A similar pattern of waning sensitization to ESAT-6/CFP-10 has been documented in naturally infected black rhinoceroses during treatment for *M. tuberculosis* infection (8). Such patterns are associated with a decrease in pathogen burden and antigenic load (17). Notably, for PB1, the peak in CMI responses at 5 months PI occurred at the only time point that viable *M. bovis* was isolated from the respiratory tract. Hereafter, the gradual decline in *M. bovis*-specific immune responses, as for the other two experimentally infected animals, was associated with limited TB-like pathology at postmortem examination and a failure to isolate viable organisms from multiple tissues. Together, these immune response patterns appear to reflect the control of *M. bovis* infection in these animals following the development of the adaptive immune response.

Humoral responses broadly mirrored those of CMI, with PB1 showing strong seroreactivity to all test antigens. However, PB2 and PB4 showed limited detectable sensitization to the *M. bovis*specific antigens MPB83, MPB83/MPB70, and ESAT-6/CFP-10. Similar differences in seroreactivity of individual animals are seen following experimental infection in other species, including cattle and badgers (18, 19). As for these hosts, this phenomenon in rhinoceroses could be related to the infection dose administered to each animal or reflect a natural animal-to-animal diversity in antigen recognition (18, 19). Nonetheless, following infection, PB2 and PB4 did display a transient elevation in seroreactivity to PPDb. In part, this may reflect the far greater range of antigenic epitopes present in PPDb compared to the more specific antigens such as MPB83 and ESAT-6. Alternatively, because PPD includes antigens that are shared between various mycobacterial species (20), the detection of immunological sensitization to PPDb might in part reflect the boosting of immune responses to cross-reactive antigens. Evidence of such cross-reactivity is indicated by the distinct CMI response of PB1 to PPDa following *M. bovis* infection.

At 20 months PI, all rhinoceroses displayed a rapid increase in IFN-γ release in response to the antigens PPDb, PPDa, and ESAT-6/CFP-10/TB7.7. This boosting of the immune response might have reflected acute progression of TB disease in these animals (5); however, this is unlikely given the paucity of *M. bovis*associated pathology at postmortem examination shortly after this sampling time point (10). Alternatively, this observation may be related to viral pneumonia observed in PB2 and PB4 (10). The immune response associated with acute viral infections can induce non-specific "bystander" activation and proliferation of lymphocytes (21), giving rise to IFN-γ production (22). Notably, at this time, PB2 and PB4 displayed a greater response to PPDa than to PPDb. This may be further evidence that this phenomenon resulted from a non-specific mycobacterial or viral infection rather than an *M. bovis*-specific event.

Results from the present study provide insight into the temporal patterns of immunological responses to *M. bovis* infection in white rhinoceroses. However, limitations in the study design must be considered when extrapolating findings to animals under natural conditions. First, the nature of the study precluded statistical analysis of data, partly because of the limited sample size and partly because of the substantial variation in immune responses between individuals. Unlike studies using inbred species such as experimental mice (13), this limitation is commonly experienced in studies investigating outbred animals, especially wildlife (19). Second, differences in the infective doses administered to each animal may have accounted for the differences in pathological outcomes and immune response profiles and these may not reflect natural infection. Nonetheless, it is notable that while the magnitude of these responses differed between individuals, their kinetics were comparable. Lastly, the three experimentally infected individuals were all relatively young, had been treated for ectoparasites, were generally free from severe stressors and were well fed (10). As such, the immunological profiles described in this study may reflect an optimal outcome of infection and not necessarily a scenario in which these and other factors might affect either disease progression or resolution.

In summary, this study has characterized antigen-specific immune response patterns of white rhinoceroses that appear to reflect the effective control of *M. bovis* infection by this host. Importantly, the measurement of these responses might prove useful for assessing the disease status or response to treatment of naturally infected animals. Moreover, the immunological markers identified in this study might be used for the detection of infection; however, further studies investigating these markers in *M. bovis*-uninfected animals are required to confirm their diagnostic utility.

### ETHICS STATEMENT

This study was carried out with the approval of the Animal Ethics Committees of the South African National Parks (SANParks) and Stellenbosch University. Animal care was in accordance with the SANParks Standard Operating Procedures for the Capture, Transportation and Maintenance in Holding Facilities of Wildlife.

### AUTHOR CONTRIBUTIONS

AM, MM, VR, and PH designed the study; MM, PB, and JH provided animal care and performed chemical immobilization and blood collection; SP, DM-L, JH, RM, and MM performed laboratory assays; data analysis was done by SP, DM-L, VR, PH, MM, and AM; the first draft of the manuscript was written by SP and all authors contributed to editorial revision.

## ACKNOWLEDGMENTS

The authors are grateful for the expert support of Guy Hausler of Stellenbosch University and Nomkosi Mathebula of SANParks Veterinary Wildlife Services (VWS), for their contribution to the processing of blood. We are also thank Dr. Markus Hofmeyr, Marius and Melandie Kruger and the VWS boma team for the husbandry and care of the rhinoceroses.

### FUNDING

Research reported in this publication was supported by the South African Medical Research Council (SAMRC) and the South African National Research Foundation (grants 41744, 86949 as well as the Incentive Funding for Rated Researchers). The content is solely the responsibility of the authors and does not necessarily represent the official views of the SAMRC or the NRF.

### REFERENCES

**Conflict of Interest Statement:** The research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

*Copyright © 2017 Parsons, Morar-Leather, Buss, Hofmeyr, McFadyen, Rutten, van Helden, Miller and Michel. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*

# Feeding immunity: Physiological and Behavioral responses to infection and resource limitation

*Sarah A. Budischak1 \*, Christina B. Hansen1 , Quentin Caudron1 , Romain Garnier1,2, Tyler R. Kartzinel3 , István Pelczer <sup>4</sup> , Clayton E. Cressler <sup>5</sup> , Anieke van Leeuwen1,6 and Andrea L. Graham1*

*1Department of Ecology and Evolutionary Biology, Princeton University, Princeton, NJ, United States, 2Department of Veterinary Medicine, University of Cambridge, Cambridge, United Kingdom, 3Department of Ecology and Evolutionary Biology, Brown University, Providence, RI, United States, 4Department of Chemistry, Princeton University, Princeton, NJ, United States, 5School of Biological Sciences, University of Nebraska, Lincoln, NE, United States, 6NIOZ Royal Netherlands Institute for Sea Research, Department of Coastal Systems, and Utrecht University, Texel, Netherlands*

#### *Edited by:*

*Andrew Steven Flies, University of Tasmania, Australia*

#### *Reviewed by:*

*Claudio Silva, Federal University of Uberlandia, Brazil Jamie Kopper, Michigan State University, United States William Parker, Duke University, United States*

> *\*Correspondence: Sarah Budischak<EMAIL_ADDRESS>
#### *Specialty section:*

*This article was submitted to Comparative Immunology, a section of the journal Frontiers in Immunology*

*Received: 29 September 2017 Accepted: 14 December 2017 Published: 08 January 2018*

#### *Citation:*

*Budischak SA, Hansen CB, Caudron Q, Garnier R, Kartzinel TR, Pelczer I, Cressler CE, van Leeuwen A and Graham AL (2018) Feeding Immunity: Physiological and Behavioral Responses to Infection and Resource Limitation. Front. Immunol. 8:1914. doi: 10.3389/fimmu.2017.01914*

Resources are a core currency of species interactions and ecology in general (e.g., think of food webs or competition). Within parasite-infected hosts, resources are divided among the competing demands of host immunity and growth as well as parasite reproduction and growth. Effects of resources on immune responses are increasingly understood at the cellular level (e.g., metabolic predictors of effector function), but there has been limited consideration of how these effects scale up to affect individual energetic regimes (e.g., allocation trade-offs), susceptibility to infection, and feeding behavior (e.g., responses to local resource quality and quantity). We experimentally rewilded laboratory mice (strain C57BL/6) in semi-natural enclosures to investigate the effects of dietary protein and gastrointestinal nematode (*Trichuris muris*) infection on individual-level immunity, activity, and behavior. The scale and realism of this field experiment, as well as the multiple physiological assays developed for laboratory mice, enabled us to detect costs, trade-offs, and potential compensatory mechanisms that mice employ to battle infection under different resource conditions. We found that mice on a low-protein diet spent more time feeding, which led to higher body fat stores (i.e., concentration of a satiety hormone, leptin) and altered metabolite profiles, but which did not fully compensate for the effects of poor nutrition on albumin or immune defenses. Specifically, immune defenses measured as interleukin 13 (IL13) (a primary cytokine coordinating defense against *T. muris*) and as *T. muris*-specific IgG1 titers were lower in mice on the low-protein diet. However, these reduced defenses did not result in higher worm counts in mice with poorer diets. The lab mice, living outside for the first time in thousands of generations, also consumed at least 26 wild plant species occurring in the enclosures, and DNA metabarcoding revealed that the consumption of different wild foods may be associated with differences in leptin concentrations. When individual foraging behavior was accounted for, worm infection significantly reduced rates of host weight gain. Housing laboratory mice in outdoor enclosures provided new insights into the resource costs of immune defense to helminth infection and how hosts modify their behavior to compensate for those costs.

Keywords: *Trichuris muris*, resource–immune trade-offs, compensatory feeding, DNA metabarcoding, nuclear magnetic resonance spectroscopy metabolite profiling, rewilding mice

## INTRODUCTION

Ecologists have long-studied energy and nutrient flows in ecosystems to understand their function and how they may respond to environmental changes (1–4). In ecological communities, these flows are often a core currency of species interactions in food webs. Parasite–host interactions represent a trophic interaction in their own right, and elucidating resource flows within infected hosts can reveal crucial processes that determine immunity and infection outcomes. Resources ingested by hosts must first be metabolized and used for maintenance (i.e., baseline metabolism), and can subsequently be used for host growth (including growth of immune cells), or be diverted to parasite growth and reproduction (**Figure 1**).

Resource flow models also capture the total cost of infection, including nutrients going to immunity, parasites, and host tissue repair. By modeling the priority of resource allocation to the host's immune system vs. resources captured by parasites, Cressler et al. (5) illustrated how increasing resource acquisition can have qualitatively different effects on host immunity

and parasite loads: the immune priority scenario (i.e., when additional resources go first into immune system components) quickly clears infections, while the parasite priority scenario (i.e., additional resources are first nabbed by the parasite), results in chronic infections. At molecular and cellular levels, immunologists are increasingly describing how nutrients and metabolites affect particular immune pathways (6–8). For example, receptors for glucose and leptin, a signal of body fat (9), are found on T- and B-lymphocytes, macrophages, neutrophils, and natural killer cells and can stimulate inflammatory responses [reviewed in Ref. (6, 10)]. While great progress has been achieved in understanding how nutrition affects immunity (7), scaling-up our understanding of resource–immune interactions from molecules and cells to entire organisms and populations remains challenging (11). Here, we use semi-natural enclosures to investigate resource flow through parasitized hosts by examining trade-offs among immunity, host condition (protein and fat stores), and parasite growth and survival under two levels of resource availability.

Because of their tractability and the plethora of tools available for studying immune pathways, laboratory (lab) mice have been integral in building our basic understanding of immunology, as well as how parasite-infected hosts respond to energy, macronutrient, and micronutrient limitation (12–17). However, when scaling to organism-level questions of nutrition and resource flows during infection, it is becoming apparent that lab mouse experiments do not recapitulate some critical biological features (**Table 1**). The *ad lib*, energy-rich, readily accessible resource conditions of lab mice differ from those of most human and wildlife populations. Additionally, ties between immunity, metabolism, and the gut microbiome are increasingly recognized (18, 19), and the low diversity microbiomes and low activation state of T cells of lab mice most closely resemble neonates (20) and diverge widely from wild mice (21). Furthermore, gastrointestinal (GI) helminths also play a role in training and modulating immunity (22, 23), with increasing risks of allergic and autoimmune conditions in human populations devoid of their coevolved worms (24). Thus, this and other studies of helminth infection in

Table 1 | Studying laboratory mice in semi-natural enclosures provides a tractable experimental system that recapitulates more aspects of wild systems than traditional laboratory experiments while avoiding confounding complexities such as unknown exposure histories and coinfections.

*a Unless breeding pairs are purposely put together.*

*bNone if housed in single-sex enclosures, but possibly reproduce if both sexes are cohoused.*

rewilded mice1 provide the opportunity to study how these key components of immunological regulation interact to affect the health of humans and other animals.

Fortunately, recent studies offer a promising compromise between realism and tractability for studying immune–resource interactions in lab mice. For example, immune traits of lab mice can be made more representative by generating more diverse microbiomes or exposure histories in the laboratory (20, 25). An even greater degree of realism can be achieved by putting lab mice in outdoor enclosures (**Table 1**) that arguably mimic aspects of their evolutionary history as commensals of humans engaged in agriculture (26). Outdoors, lab mice develop more diverse microbiomes, elevated T cell responses and higher parasite loads compared with indoor lab mouse controls (See text footnote 1). Additionally, in outdoor enclosures, mice experience more natural variation in activity and thermo-regimes that may make trade-offs between immunity and other physiological processes more apparent than under lab conditions (**Table 1**).

Host behavior is central to scaling-up immune–resource interactions from the cellular to organismal level. Mice on low-quality (i.e., protein) diets may consume a greater quantity of chow, which can alter their body fat composition and immune profiles (27). From livestock studies, we know that GI nematodes, or host immune responses to them, can reduce host appetite, decreasing the energy budget the host has to allocate to defense, repair, and other physiological functions (28). Reciprocally, hosts can alter foraging during infection by consuming medicinal plants (29) or by increasing foraging to mitigate costs of infection and immunity (30). In free-ranging populations, increased foraging may come with increased energetic costs, parasite exposure, social stress, and predation risks (11, 28, 31, 32). Thus, assessing feeding behavior is critical for understanding how organisms respond to resource limitation and the resultant fitness consequences (7).

Using a semi-natural field system (See text footnote 1), we manipulate resource availability to examine the resource costs of infection and immunity. Our goal was to investigate the costs of infection (including resources diverted to parasites and immunity) and to learn how hosts may use foraging behavior to mitigate those costs in lab mice (C57BL/6) infected with the GI nematode *Trichuris muris*. *T. muris* is a whipworm that lives in the cecum, and is a congener of the parasite *Trichuris trichiura* that infects over 450 million people worldwide (33). To assess how hosts respond to infection under resource limitation, we manipulate levels of dietary protein, which is known to have strong effects on host immune defenses (7, 17, 34–37). We might expect mice fed the high-protein diet to have stronger immune responses and lower parasite loads than those on the low-protein diet. Alternatively, mice on the high-protein diet could tolerate infection while mice on the low-protein diet resist (17), which would turn the expected observation around and show stronger immune responses and lower parasite loads in the low-protein treatment. At the individual level-scale, we track whether mice compensate for potential joint costs of infection or a low-protein diet by altering foraging behavior. We predict that there will be trade-offs between food intake, investment in immunity, and parasite load (**Figure 1**). Infected mice may either feed less due to infection-induced inappetence (28, 38, 39), or feed more to compensate for costs of infection and immunity. Our findings do indeed suggest complex repercussions of behavioral changes for the flow of resources into infection and immunity.

### MATERIALS AND METHODS

### Experimental Design

Our experiment to examine trade-offs between infection, immunity, and within-host resources (i.e., diet quality) included four treatment groups: high-protein infected, high-protein uninfected, low-protein infected, and low-protein uninfected. Eightyeight female C57Bl/6 mice aged 5–6 weeks were obtained from Jackson Laboratories and individually identified with both ear tags and RFID tags (see below). All animal care was conducted in accordance with protocols approved by the Princeton University Animal Care and Use Committee (Protocol no. 1982-14). Mice were housed in groups of five in the laboratory and randomly assigned to the two diet treatments as well as two cohorts that were staggered by 2 days for logistical purposes (**Figure 2**). The high-protein (HP; 20%; Envigo Teklad Custom Diet TD.91352) and low-protein (LP; 6%; Envigo Teklad Custom Diet TD.90016) diets had the same energy density (3.8 kcal/g) and micronutrient composition. The typical chow fed to lab mice (e.g., PicoLab® Rodent Diet 20) has very similar composition to the HP diet. For 10 days mice were fed the assigned diets in the lab while temperature and light cycles were gradually altered to mimic outdoor conditions (June/July in New Jersey, USA: 26 ± 1°C with a 15-h light–9-h dark cycle; **Figure 2**).

Next, the 22 mice in each diet–cohort combination were transported to four outdoor enclosures (**Figure 2**), two of which contained the HP chow and the other two the LP chow (**Figure 2** inset). The enclosures are replicate pens of approximately 180 m2 , fenced in by zinced iron walls extending 1.5-m high and 80-cm deep, and topped with electric fencing and reflective aluminum pans to deter ground and aerial predators, respectively (See text footnote 1). Diets were provided *ad libitum* at feeding stations monitored by two RFID readers to track each individual's time spent feeding. The natural environment could serve as an additional source of food (e.g., berries, seeds, insects). Each enclosure provided two watering stations inside a small (180 cm × 140 cm × 70 cm) straw-filled shed (See text footnote 1). Mice were weighed at the start of the experiment, the day of release, and approximately weekly thereafter as they were trapped overnight in the outdoor enclosures using chow-baited Longworth traps. At each weekly trapping, fecal samples were collected and blood samples were taken *via* shallow cutaneous tail snips into heparinized capillary tubes.

Mice were acclimated to the enclosures for 2 weeks (14–17 days) prior to *T. muris* infection (**Figure 2**). A dose of 200 embryonated *T. muris* eggs (strain E) was then given *via* oral gavage to the first 16 mice trapped per enclosure. If more than 16 mice

<sup>1</sup>Leung JM, Budischak SA, Chung HT, Hansen C, Bowcutt R, Neill R, et al. The shock of the new: rapid environmental effects on gut nematode susceptibility in re-wilded mice. Under review.

were trapped on a given day, individuals were randomly assigned to infection treatment. In total, 29 mice on the HP diet and 31 mice on the LP diet were infected. Thus, infected and uninfected animals were cohoused in the same enclosures. The remaining mice (15 HP, 13 LP) served as uninfected controls. Nematode infection could not be transmitted between mice assigned to different treatments in the shared enclosures due to the long life cycle of *T. muris* [>28 days to maturity (40)], and the relatively short duration of the experiment [<26 days postinfection (dpi)].

Approximately 3 weeks postinfection (mean ± SE: 21.5 ± 1.7 dpi), range: 19–26 (dpi), all mice were trapped, weighed, and transported back to the laboratory (**Figure 2**). Mice were anesthetized via isoflurane inhalation followed by terminal cardiac puncture. Cardiac blood was drawn from each mouse, spun in a heparinized tube, and plasma was separated and stored at −80°C. During necropsy, mesenteric lymph nodes (MLNs) were excised for cell culture (see below). Finally, to determine carcass weight, all major organs (except the brain) were removed and weighed. The spleen, large intestine (emptied of contents), and cecum were separated and weighed. Ceca were frozen for later dissection and parasite enumeration.

### Quantifying Feeding Behavior

Individual-level feeding behavior was assessed using a custom-built monitoring system (Datasheet S1 in Supplementary Material). To track feeding behavior, Radio Frequency IDentification tags (RFID; 8 mm × 1.4 mm FDX-B "Skinny" PIT Tag, Oregon RFID, Portland, OR, USA) were injected subcutaneously. When in the vicinity of an antenna, the RFID tag emits a unique string of numbers that serves as an individual's identifier (i.e., RFID number). In each enclosure, Feeding Event Tracking Apparatuses (FETA) with RFID antennas were deployed. The FETA transmit a signal to an Event Acquisition and Reporting System (EARS), which automatically compiles the data produced by the FETA (FETA number, RFID number, and timestamp). Two FETA were placed in sequence and connected on one end to the sole entrance to the chow hopper.

The directionality of mouse movement through the sequential FETA boxes was used to determine when mice were in the feeder. Specifically, mice were inferred to be in the feeder for the duration of time between two EARS readings on the FETA nearest the chow hopper. Intervals <5 s in duration (7.5% of visits) were removed because they allowed insufficient time to feed. Intervals greater than 2 h were also removed because mice were likely not eating for such a long duration and this small fraction (<0.1%) of visits skewed feeding time distributions. Intervals when chow was not available in the feeder (i.e., during trapping sessions) were also dropped. Next, total time feeding was summed per mouse, then divided by the number of days each mouse was in the enclosure because some mice were in the enclosures for longer than others (35–40 days depending on when they were trapped). This produced a comparable measurement of time spent feeding (min) per experiment day. Mice that lost their RFID tags (*n* = 6) and two individuals that did not eat at the feeder (working RFIDs but stopped visiting feeder) were excluded from feeding behavior analyses. Chow consumption was monitored by weighing each chow hopper every time it was removed for trapping as well as before and after refilling (i.e., every 1–3 days).

### Immune Response Measurements

Following established protocols (See text footnote 1, 41), MLN tissue was excised during necropsy, prepared into single-cell suspensions at a density of 5 × 106 cells/mL, stimulated with *T. muris* antigen at 5 µg/mL, and cultured for 48 h. Supernatants were collected and analyzed in duplicate for interleukin 13 (IL13), interleukin 10 (IL10), interleukin 17 (IL17), and interferon gamma (IFNg) using half-reactions of Beckton Dickinson Cyometric Bead Array Mouse/Rat Soluble Protein Flex Set system (BD Biosciences, Oxford, UK) and an LSRII flow cytometer (BD Biosciences) (See text footnote 1). Concentrations were analyzed with the FCAP Array software (version 3.0.1, BD Biosciences).

Worm-specific IgG1 titers were measured using sandwich enzyme-linked immunoassays (ELISA). Ninety-six-well plates were coated with *T. muris* excretory-secretory (ES) antigen at a concentration of 5 µg/mL in carbonate buffer (0.06 M, pH = 9.6) with 2% bovine serum albumin (BSA). Plates were then blocked overnight with 2% powdered milk in carbonate buffer. After a 2-h incubation at 37°C, plates were washed three times with *tris*buffered saline-Tween (TBST). Plasma samples were added and serially diluted from 1:50 to 1:6,400 with TBST and incubated for 2 h at 37°C. After washing five times with TBST, horseradish peroxidase (HRP) conjugated IgG1 antibody (1:4,000 in TBST) was added. Following a 1-h incubation at 37°C, plates were washed five times and Substrate ABTS (Sigma-Aldrich) was added. Plates were developed for 20 min at 37°C, then stopped with a 1% sodium dodecyl sulfate solution. Absorbance was measured at 405 nm using a Multiscan™ GO spectrophotometer (Thermo Scientific). Titers were determined by the dilution at which absorbance exceeded 4 SDs above background, defined as the plate-wide average absorbance of the two most dilute concentrations of each sample.

To quantify total IgG, a mouse antibody IgG pair was purchased from StemCell Tech (Catalog no. 01998C). Plates were coated at 1 µg/mL 50 μL per well with carbonate buffer (pH 9.6) overnight at 4°C. Plates were blocked for an hour with TBST 20 with 0.1% BSA at 37°C, then washed. Preliminary assays revealed that IgG levels were quite high, so to fall in the range of the standards, samples were diluted 1:81,920 in TBST containing 0.1% BSA. After incubation for 2 hat 37°C and washing, a secondary antibody was applied at a concentration of 1:1,000 in TBST with 0.1% BSA. Plates were incubated at 37°C for 1 h then washed. Next, we added 100 µL/well of p-nitrophenyl phosphate (pNPP) substrate to all wells. Finally, plates were incubated for 30 min and read at 405 nm. Concentrations were calculated in comparison to a plate-specific standard curve (all *R*<sup>2</sup> > 0.999).

### Parasite Quantification

Caeca were cut open longitudinally and examined under a dissecting microscope by an observer blind to infection status. Worms were isolated by scraping the gut mucosa into a series of clean petri dishes of water. After enumeration, worms were stored in 70% ethanol for subsequent length measurements. Each worm was photographed under a dissecting scope (2–4× power) and ImageJ (version 1.49, NIH, USA) was used to measure total length.

### Condition and Nutritional Measurements

Plasma albumin concentrations were measured colorimetrically. QuantiChrom BCG albumin assay kits (BioAssays) were used following the manufacturer's instructions. Briefly, 5-µL duplicates of each standard and sample (diluted 1:2 with ultrapure water) were mixed with 200 µL of BCG reagent. Plates were incubated at 23°C for 5 min before being read at 620 nm. Concentrations were determined in comparison to a standard curve run in duplicate (*R*<sup>2</sup> = 0.997).

Plasma leptin concentrations were analyzed using a RayBio® Mouse Leptin ELISA kit following manufacturer's instructions (RayBiotech, Norcross, GA, USA). Briefly, 10 µL of plasma samples were diluted 1:10, incubated on a 96-well plate coated with mouse leptin antibody. After washing, a biotinylated anti-mouse leptin antibody was added. Next, the plate was washed, horseradish peroxidase-conjugated streptavidin was added. Following another washing step, a buffered 3,3,5,5′-tetramethylbenzidine solution was used to produce a color change reaction stopped after 30 min with 0.2-M sulfuric acid. Color intensity was read at 450 nm. Concentrations were calculated using a standard curve (standards run in duplicate, *R*<sup>2</sup> = 0.995).

### Characterizing Herbivory on Wild Plants Using DNA Metabarcoding

To compare the wild plant species diversity and composition with which lab mice supplemented their diets, we employed a DNA metabarcoding method that involves sequencing and identifying undigested plant DNA obtained from fecal pellets (42, 43). The plant DNA in fecal pellets is likely to reflect very recent foraging activity because the half-life of ingesta in lab mice is approximately 74 min (44), meaning that <1% of contents are retained after 8 h. Thus, the resulting dietary profiles represent plants eaten over the ~2.5- to 5-h period prior to sample collection (44, 45).

Fecal samples from 26 mice (*n*= 7 HP inf, 9 LP inf, 6 HP uninf, 4 LP uninf) from the end of the study were obtained for dietary analysis. Samples were frozen on dry ice and stored at −80°C to preserve DNA prior to extraction. Total DNA from 1 to 2 pellets (~15–30 mg) per mouse was extracted using a Zymo Xpedition Soil/Fecal DNA mini kit with an extraction blank to monitor for potential cross-contamination. Using PCR, the P6 loop of the chloroplast *trnL*(UAA) marker was amplified with primers *g* and *h* (42). The PCRs were run with unique combinations of the *g* and *h* primers that had been modified with 8-nt multiplex identification (MID) tags in order to enable pooling of amplicons for sequencing using established protocols (46). These PCR products were cleaned and normalized using SequelPrep normalization plates, then pooled for sequencing in a 170-nt single-end run of the Illumina HiSeq2500 at Princeton University's Lewis Sigler Institute following Kartzinel et al. (46).

The resulting DNA sequence data were assigned to samples of origin (i.e., demultiplexed), screened to reduce potential sequencing errors, and identified by comparison to a plant DNA reference library. The sequences were demultiplexed and primers were removed using the *ngsfilter* command in the Obitools software (47). We discarded sequences that contained ambiguous base calls (i.e., non-A, T, C, or G characters) that were <9-nt long or that had mean Illumina quality scores of <32. Unique sequences were merged and tabulated within samples to permit quantification of DNA sequence relative read abundance (RRA), a measure of the proportion of each dietary plant species in each dietary sample. Putative errors were screened by using the *obiclean* command to identify sequences that differed by 1 nt from another sequence in the same sample, but that occurred at <5% of the abundance. These sequences were removed from further analyses. Plant DNA was identified by comparison to a reference library developed using the European Molecular Biology Laboratories (EMBL) (Database release no. 130). From this archive, we extracted 35,776 unique sequences (229,430 entries) with ≤3 mismatches to the *trnL*-P6 primers *g* and *h*. Unique dietary sequences were identified by comparison to this reference database using the *ecoTag* command. Operational taxonomic units (OTUs) were identified by clustering at the >97% level using the uclust algorithm (48). To focus on abundant and well-identified plant taxa, we considered only OTUs with >80% identity to the EMBL database and those representing >5% of reads within each sample. Samples were rarefied to even sequencing depth using phyloseq (49) in R (50).

We inspected taxonomic identifications to identify imprecision that could arise from gaps or misidentified DNA sequences in the reference library. Any OTUs that exactly matched a taxon not known to occur in the region were revised to higher taxonomic levels, and the set of references matching any OTUs that were poorly identified (family or higher taxonomic levels) were scrutinized for taxonomic outliers (Table S2 in Supplementary Material). Only two plant taxa were included in the manufacturing process of the chow provisioned to lab mice in this experiment were corn (*Zea mays*) and beets (*Beta vulgaris*), and DNA from these plant products was expected to be destroyed by irradiation during the chow production process; indeed, no DNA sequences that match either of these plant species were identified in the final set of OTUs from our analysis. We calculated RRA by converting the rarefied number of reads into a proportion of reads per sample (i.e., ranging 0–1). Although RRA is not always a reliable measure of proportional dietary utilization in DNA metabarcoding studies (51), the analysis of RRA based on the *trnL*-P6 protocol employed in this study has been supported at least to the level of plant family and functional group in independent studies from different systems (46, 52, 53).

### NMR of Dietary Metabolites

Nuclear magnetic resonance spectroscopy (NMR) was used to examine dietary metabolites from a random subsample of uninfected mice (12 HP and 8 LP) with sufficient frozen fecal samples from the end of the experiment. Fecal pellets from the mice were weighed, crushed, and diluted 1:3 with phosphate buffered saline (PBS). Samples were vortexed, allowed to dissolve overnight, and re-vortexed. After centrifugation, the supernatant was decanted and brought to ca. 30 µL with PBS in a 1.7-mm OD capillary tube (New Era Enterprises, Vineland, NJ, USA). This capillary was inserted into a 5/2.5-mm OD NMR tube containing small amount of D2O, which served as an external lock material. Samples were analyzed on an 800- MHz Bruker Avance III HD NMR spectrometer equipped with a custom-made 1 H/19F/13C/15N//2 H cryoprobe using Topspin v.3.2. Water suppression was done using the first increment of the NOESY pulse sequence (delay–G1–90o –t1–90o –tm–G2–90o – acquisition) applying weak presaturation during the relaxation delay and the 100-ms mixing time period. Data processing was done offline using MNova v.11.0 (Mestrelab Research, Santiago de Compostela, Spain). After zero filling, apodization with 1-Hz Gaussian broadening, careful phase and baseline correction were applied manually. In order to compensate for variance of concentrations the spectral intensities were normalized to total integral, excluding the small region of the residual water signal. Local peak alignment was applied where necessary and possible using the inherent function in MNova. A few key metabolites were identified based on literature data (54). Out of the 20 samples, two were discarded because of technical problems (poor shimming and water suppression) to maintain statistical integrity of the residual data. The collection of 18 spectra was then exported to a spreadsheet.

### Statistical Analysis

The final sample size of this study was 80 (LP infected = 31, LP uninfected = 10, HP infected = 24, HP uninfected = 15). Several mice eluded capture for over 20 days beyond when the rest were trapped and sampled; these were excluded from statistical analyses. An additional mouse that had a severe congenital uterine defect was also excluded.

We assessed effects of diet and infection on immunity and condition. No uninfected mice had detectable IL13 concentrations, so a Kruskal–Wallis rank sum test was first used to compare infected and uninfected mice. The other cytokines (IFNγ, IL10, IL17) also had highly skewed distributions that could not be normalized. Thus, effects of diet and infection were analyzed using Kruskal–Wallis tests. Next, a general linear model (GLM) was used to test the effect of diet on log-transformed IL13 concentration among infected mice. GLMs were also used to test effects of diet and infection on log-transformed worm-specific IgG1, total IgG, and spleen weight relative to carcass weight. Using GLMs, we next tested for effects of diet and infection on mouse condition, including weight gain per day, plasma albumin concentration, plasma leptin concentration, and carcass weight. Plasma albumin and leptin concentrations were log transformed for the analysis. Diet–infection interaction terms were dropped when non-significant. To account for potential correlation structure among plasma components, we conducted a principal component analysis (PCA) of worm-specific IgG1, total IgG, albumin, and leptin. All components were log transformed and scaled prior to analysis. Next, the relationship between diet and infection status and each principal component (PC) was tested using GLMs. This PCA yielded no new insights beyond the univariate analyses described above, and results can be found in Table S1 in Supplementary Material.

To assess differences in fecal metabolites between the diet treatments, the SIMCA-P v.12.1 software package was used (Umetrics, Umea, Sweden). To determine if there was spectrum-wide variation between the diets, partial least squaresdiscriminant analysis (PLS-DA) was performed and orthogonal atrial least squares-discriminant analysis (OPLS-DA) were performed. Widely used in metabolic phenotyping studies, the PLS-DA analyses are better able to detect clusters than traditional PCA, and OPLS-DA provides even stronger discrimination based on known groupings (e.g., diet) (55). The DA methods were validated by calculating *R*<sup>2</sup> and *Q*<sup>2</sup> . Prior to PLS-DA and OPLS-DA analysis, both UV and Pareto scaling were applied. UV scaling is a better choice for maintaining uniform contribution from all peaks regardless of their absolute intensity, while Pareto scaling provides access to "spectrum-like" loadings plot.

To examine morphological and behavioral changes, we first tested for effects of diet and infection on large intestine size (tissue weight/carcass weight) using a GLM. Similarly, we used a GLM to determine if time spent feeding per day in the enclosures varied with diet or infection status. Third, we examined supplemental foraging on natural plants in the enclosures using fecal DNA metabarcoding data. We tested for significant differences in the composition of plant species eaten according to diet, infection status, and leptin levels (log-transformed) using permutational MANOVA (perMANOVA) in *vegan* (56) in R. For the per-MANOVA, we calculated pairwise Bray–Curtis dissimilarity metrics for each pair of samples, which ranges from 0 to 1 (from completely identical to mutually exclusive dietary compositions, respectively). For the six wild-plant families with the greatest overall RRA, we compared mean RRA between treatment groups (infected vs. uninfected and LP vs. HP diets) and across levels of leptin. Exploratory trend lines were fit to the data using generalized linear modeling. We investigated differences in RRA between genera of the legume family (Fabaceae) in more detail because this family comprised greatest overall RRA.

Lastly, we tested for differences in parasite load and weight gain, corrected for the amount of time each individual spent feeding. Among infected individuals, a non-parametric Kruskal– Wallis test was used to determine if diet affected worm count. In infected individuals, we also tested for an effect of diet on worm length (log transformed for normality) with a GLM. The amount of weight each mouse gained while in the enclosures was divided by the amount of time each spent feeding. Weight gain per minute feeding was compared between diets and infection status using a GLM. All analyses were run in R version 3.1.2.

### RESULTS

### Dietary Quality and Infection Affected Immunity

Measures of immune function were affected by both diet and infection. MLN cells of the uninfected mice produced no detectable IL13 [a cytokine strongly induced by nematodes (57)] in response to stimulation with nematode antigen, making them significantly lower than the infected mice (Kruskal–Wallis χ<sup>2</sup> = 7.52, df = 1, *p* = 0.0061; **Figure 3A**). Among infected mice, those on the high-protein diet had higher IL13 concentrations in culture supernatants than those on the low-protein diet [Est (±SE): 0.56 ± 0.25, *t*= 2.20, *p*= 0.033; **Figure 3A**]. The cytokines IFNγ and IL17 mediate pro-inflammatory responses, primarily to intracellular pathogens and extracellular bacteria and fungi, respectively (58, 59). IL10 mediates anti-inflammatory, regulatory responses (58, 60). Surprisingly, IFNγ, IL17, and IL10 were also higher in infected mice but did not differ between diets (Figure S1 in Supplementary Material). In mice and many other mammals, IgG1 is an important antibody class for fighting GI nematode infection (57). In addition to being higher in infected mice (*p* = 0.00006), *T. muris* specific IgG1 titers were elevated in mice eating the high-protein diet (*p* = 0.032; **Table 2**; **Figure 3B**). Total IgG concentrations did not vary with diet or infection status (all *p* > 0.33; **Table 2**; **Figure 3C**), so we cannot account for the elevated IgG1 titers observed in uninfected mice eating the high-protein diet as a correlate of high total IgG concentrations. Spleen size also did not differ significantly with infection status or diet (all *p* > 0.10; **Table 2**; **Figure 3D**).

G1 (IgG1) were affected by diet and infection, but not (C) total IgG concentration or (D) spleen size (weight/carcass weight). Asterisks denote significant effects of diet or infection (Inf).

*Significant effects are highlighted in bold.*

### Dietary Quality Affected Nutrition and Condition

Dietary protein affected multiple aspects of mouse nutrition and condition, whereas parasite infection did not. Although neither diet (*p* = 0.19) nor infection (*p* = 0.052) significantly affected rates of weight change (**Figure 4A**; **Table 2**), other condition measures were more sensitive. Mice on the LP diet had significantly lower albumin concentrations (*p* = 0.006; **Figure 4B**; **Table 2**) than those on the high-protein chow. The low-protein diet also led to elevated leptin concentrations, a metabolic and immune-regulatory hormone released in proportion to body fat, as well as to higher carcass weights (*p* = 0.044; **Figures 4C,D**; **Table 2**). Infection did not affect any of these other condition measures (all *p* > 0.15; **Table 2**).

### Dietary Quality Affected Fecal Metabolites

Nuclear magnetic resonance spectroscopy provides data and information on metabolites at the molecular level. The average spectrum of the 18 fecal samples, after normalization, peak alignment, and identification of some components (54), is shown in **Figure 5A**. PLS-DA (UV scaled) revealed distinct clustering by diet (**Figure 5B**), with a high proportion of variance explained by PC1 and PC2 (*R*<sup>2</sup> *Y* = 0.993) and decent predictability (*Q*<sup>2</sup> = 0.533). The OPLS-DA (Pareto scaled) coefficient plot (**Figure 5C**) shows many positive and negative intensities, corresponding to metabolites in higher or lower relative concentrations in LP and HP diets, respectively. A great variety of sugars are present in the fecal extract, and aliphatic amino acids are clearly in higher abundance in the LP diet group (**Figure 5C**). The differences in metabolite relative abundances between the diet treatment indicates the LP diet produces clear alterations in the metabolism of the mice and/or their gut microbiota.

### Dietary Quality and Infection Altered Aspects of Morphology and Behavior

The morphology and behavior of the mice was affected by diet and infection. First, infected mice had significantly heavier large intestines (emptied large intestinal tissue, weight relative to carcass weight) than uninfected mice (*p* = 0.022; **Figure 6A**; **Table 2**). Diet and infection status affected cecum size, with larger ceca in infected mice on the HP diet (diet: *p* = 0.017, infection: *p* = 0.050; **Figure 6B**; **Table 2**). Additionally, mice on the LP diet spent more time feeding than those on the HP diet (*p* = 0.0004; **Figure 6C**; **Table 2**).

### Herbivory on Wild Plant Species Varied among Individuals

Our DNA metabarcoding strategy yielded a total of 1,433,421 demultiplexed sequence reads of high quality, 14,521 of which were unique. After removing putative sequencing errors and picking OTUs, our cleaned-up database comprised 3,017 OTUs representing 1,296,521 sequence reads (>90% of the original). Sequence counts per sample were uneven (range = 1,648– 162,180, mean = 49,702), and all were much higher than the extraction blank (*N* = 229). After removing sequences that poorly matched the reference database and that were never >5% of reads in any sample, we were left with 82% of the raw DNA sequence reads (1,167,514; of raw DNA sequence reads). These raw DNA sequence reads represented 26 plant OTUs in subsequent analyses (Table S3 in Supplementary Material). We rarefied samples to an even depth of 1,301 sequences/sample. The six most heavily utilized wild plant families (RRA > 0.05; Table S3 in Supplementary Material) included legumes (Fabaceae, cumulative RRA across all samples = 0.32), grasses (Poaceae = 0.18), wood sorrel (Oxalidaceae = 0.12), roses (Rosaceae = 0.08), violets (Violaceae = 0.07), and asters (Asteraceae = 0.07).

Individuals varied considerably in the composition of wild plants eaten, but the overall composition of wild foods eaten did not differ between treatment groups. The composition of consumed wild plants did not differ by diet quality (perMANOVA: pseudo-*F*1,22 = 0.83, *R*<sup>2</sup> = 0.03, *p* > 0.05) or infection status (pseudo-*F*1,22 = 0.96, *R*<sup>2</sup> = 0.04, *p* > 0.05; **Figure 7A**) of the mice, at least in part reflecting the high inter-individual variation in diet composition (mean pairwise Bray–Curtis dissimilarity = 0.89). Diet compositions were more closely associated with leptin levels (pseudo-*F*1,22 = 1.61, *R*<sup>2</sup> = 0.06, *p* = 0.06; **Figure 7B**), although this trend was marginally non-significant.

Figure 4 | Diet, but not infection status, affected most measures of condition. (A) Weight change over the course of the experiment (corrected for no. of days in the enclosure) was not affected by diet or *T. muris* infection. However, the LP diet led to reduced (B) albumin concentration and increased (C) carcass weight and (D) leptin levels. Asterisks denote significant effects of diet.

Nevertheless, some variation among groups in their proportional utilization of different plant types was apparent. Most strikingly, legumes (family Fabaceae) were virtually unutilized by LP-uninfected lab mice even though the mean RRA of legumes eaten by mice in other treatments ranged from ~0.3 to 0.4 (**Figure 7A**). Within Fabaceae, clover (*Trifolium* sp.) was proportionally more utilized by individuals assigned to all groups except the LP-uninfected group, while beggar's lice (*Desmodium* sp.) was proportionally more utilized by infected individuals (**Figure 7C**). Legumes tended to have higher RRA in the diets of individuals with high leptin levels, and there was a trend of decreasing Oxalidaceae and Violaceae RRA with leptin (**Figure 6B**). Lab mice from the LP-uninfected treatment utilized proportionally more Oxalidaceae (mean RRA ~0.3 vs. <0.1) and Violaceae (mean RRA ~0.2 vs. <0.1; **Figure 7A**). These trends did not, however, reflect a significant relationship (*p* > 0.05).

### Effects of Infection Visible Despite Low-Intensity Infections

Perhaps due to the greater chow consumption of mice on the LP diet and increased large intestine size of infected mice, the net effects of diet and infection reveal the dynamic complexity of scaling-up this host–parasite interaction to the individual level. For example, worm counts did not differ by diet treatment (Kruskal–Wallis χ<sup>2</sup> = 0.27, df = 1, *p* = 0.60; **Figure 8A**). Worm length also did not vary with dietary quality (*t* = 0.57, *p* = 0.57).

Costs of infection were nonetheless detectable if individuallevel feeding behavior was accounted for. Although overall weight change did not vary among treatments (**Figure 4A**), mice on the LP diet spent more time feeding (*p* = 0.0004; **Figure 6C**). Examining weight gain per time feeding revealed that infected mice gained weight at a slower rate than the uninfected mice (*p* = 0.042; **Figure 8B**; **Table 2**), suggesting that they had to invest that food energy into something other than growth (e.g., the immune response) or that parasites usurped it.

### DISCUSSION

Above all, this study shows that a broad perspective on the resource demands of parasites and immunity is needed to understand mammalian defense. In order to gain this understanding as well as insight into the ecology and evolution of host defenses in general, we must incorporate behavioral as well as physiological responses to resource limitation and infection. Using diverse data types, we discovered that within-host dynamics of infection and defense were strongly impacted by the interactions of the host with its wider environment. Most importantly, effects of infection, in terms of reduced weight gain, were only visible after accounting for variation in individual feeding behavior,

with identification of selected metabolites. (B) PLS-DA scores plot (UV scaling); the subgroups of mice on HP and LP diets are clearly separated into distinct clusters. For three components *R*<sup>2</sup> *Y*(cum) = 0.993 and *Q*<sup>2</sup> (cum) = 0.533, showing decent validity of the statistics. The ellipse denotes Hotelling's T2. (C) Loading data along the NMR spectrum (Pareto scaling) reveals that there are a great number of metabolites, which are present in distinct quantity in the separated clusters of samples. All the negative intensities belong to peaks of metabolites, which are present in greater quantity in the cohort on HP diet (green), while the positive intensities depict metabolites in larger concentration in the LP diet group (blue), respectively. Some tentative assignments are shown on the plot. Abbreviations: AAs, amino acids; Ala, alanine; Cho, aldehydes; Gly, glycine; Ile, isoleucine; Lac, lactones; Lys, lysine; nuc, nucleic acids; Phe, phenylalanine; Tyr, tyrosine; succ, succinate; Val, valine.

highlighting the complexity of resource–immune–infection relationships at the individual scale. In the sections below, we discuss each aspect of our results in detail before returning to the broader implications in a concluding section.

### Dietary Quality and Infection Affected Immunity but Not Parasite Loads

The 6% protein (low protein) chow reduced investment in both IgG1 and IL13. The direction of the effect of dietary protein on IgG1 responses to *T. muris* E/S antigen, the predominant antibody response to primary infection (61), was difficult to predict since both higher and lower responses are reasonable given previous studies. In the lab, higher IgG1 antibody concentrations during protein restriction were documented in mice infected with *T. muris* (14). However, mice infected with the nematode *Heligmosomoides polygyrus* experience lower levels of total IgG1 (15) when fed a 3% protein diet. Interestingly, a 7% protein diet permitted total IgG1 levels indistinguishable from those in mice fed a 24% protein diet (15). Similarly, the reduction in IL13 on the LP diet was not a certainty; in a previous laboratory experiment, IL13 concentrations did not vary with dietary protein level in mice infected with *H. polygyrus* (27). Spleen size was not sensitive to dietary protein, a somewhat unexpected result given that low-protein diets reduce spleen size during *H. polygyrus* infection (15, 34). However, differences between these nematodes in their infection sites (small intestine vs. cecum) and the immune responses they typically induce in C57BL/6 mice (Treg vs. Th2) (62), could account for their different effects in protein-limited hosts.

size (emptied of contents, relative to carcass weight), whereas (B) the HP diet was associated with a heavier cecum relative to carcass weight. (C) Mice on the LP diet spent more time feeding than mice on the HP diet. Asterisks denote significant effects of diet or infection (Inf).

The reductions in IgG1 and IL13 could be due to direct effects of protein limitation on T- and B-cell function, which are known to be regulated by host metabolic activity (63, 64). Indeed, trade-offs between markers of protein nutrition and wormspecific antibodies to the nematode *Teladorsagia circumcinta* were visible in a wild Soay sheep population. Moreover, these trade-offs predicted overwinter survival with nutrition increasing survival for older individuals while investment in immunity led to greater survival odds for young sheep (37). Alternatively, the elevated carcass weights and leptin concentrations of the LP mice, due to increased chow consumption, suggest they had a higher body fat content than HP mice. Food intake in uninfected mice may be regulated by dietary protein content (65). Increased chow consumption, weight, and leptin levels have also been detected in mice fed a low-protein diet in previous nematode infection-diet manipulation experiments (27, 66). Carcass weight has been posited as an indicator of tolerance during helminth infection (67), but our data suggest that compensatory feeding may disrupt that association. Obesity is associated with pro-inflammatory responses (68), which may be mediated in part by leptin (6, 69, 70). Thus, the reduced IgG1 and IL13 levels could be a consequence of the fatter LP mice shifting from Th2 responses to pro-inflammatory Th1 responses. Although it is a common finding across multiple nematodes (27, 66), it remains unknown whether protein limitation itself or increased body fat due to limitation-induced overeating drive the pro-inflammatory shift during protein limitation.

Given how integral IgG1 and IL13 are to the development of an effective Th2 response to *T. muris* infection (61, 71), it is surprising—and contrary to our initial hypotheses regarding potential relationships between diet and immunoparasitological outcome—that parasite load did not differ between the dietary treatments. IgG1 levels in uninfected HP mice were much higher than those in uninfected LP mice, which contributed to the overall effect of diet on IgG1 levels.  Yet, among infected mice, IL13 concentrations were over 70 times higher in mice eating the HP diet, so it is unclear why that did not translate to differences in parasite load. Across both diets, infected mice had cleared most of their parasites by the time they were sampled at the end of the study. Due to this clearance rate and the high numbers of individuals with below-detection immune responses, there was insufficient statistical power to examine individuallevel variation in immunity and its relation to parasite load. The treatment-level patterns suggest that perhaps the lower concentrations of IL13 in the mice given the LP diet were also sufficient to reduce *T. muris* survival by that time point, while those on the HP protein treatment had excess expression. Alternatively, the higher IL13 concentration (generated by stimulating MLN cells with *T. muris* antigen), might indicate the strength of the memory response. Thus, mice on the HP diet might be better protected during reinfection, a pattern also seen during protein limitation and infection with *H. polygyrus* (27, 34).

The low worm counts are surprising, given the 5.5-fold higher loads observed at a similar time point in a prior experiment (mean ± SE; Leung et al.: 34.2 ± 7.7 worms, this experiment: 6.1 ± 1.9 worms) in these same enclosures with similar mouse age and weight at infection, enclosure acclimation time, experiment duration, and parasite dosing performed by the same technician (See text footnote 1). The number of dpi was a strong predictor of worm counts, which declined quickly over time. If sampled sooner, differences in burden between diet treatments may have been more detectable. However, this rate of decline was indistinguishable from the prior experiment (See text footnote 1), and thus cannot explain the overall lower worm counts.

Instead, differences in hatching ability between batches of *T. muris* eggs, differential availability or composition of supplementary forage, slightly different sampling time points and/or effects of the different brand and composition of chow could contribute to the disparity in worm counts between experiments. For example, although the chow used by Leung et al. (See text footnote 1) had a similar macronutrient composition to the HP diet in this study, the highly refined ingredients in the specialty diet might have altered the gut microbial flora, reduced *T. muris* hatching,

made the intestines a less hospitable environment for *T. muris*, and/or increased the efficacy of mouse immune defenses. This is an important area of future inquiry. In any case, IgG1 and IL13 data confirm that the mice in the current study were infected for long enough to stimulate an immune defense to *T. muris*.

## Dietary Quality Affected Fecal Metabolites

In the uninfected mice, the dietary quality significantly altered the nutrient environment within the mouse; in infected hosts, *T. muris* would likely experience similar nutrient alterations. Estimating cecum nutrient content from those in feces using NMR revealed that sugar metabolites varied greatly between diets, and amino acids like phenylalanine and alanine were higher in the HP group. This is not surprising given that the difference in protein between the treatments was compensated with a higher percentage of carbohydrates to achieve equal calorie density. Relative abundances of tyrosine, which is involved in the regulation of immune signaling pathways (72), and the aliphatic amino acids tended to be higher among mice on the LP diet. This preliminary exploration demonstrates that NMR metabolite

*Trifolium* (clover) was common in all but the LP-uninfected treatment and an OTU-representing *Desmodium* (beggar's lice) was eaten only by infected mice.

analysis is a powerful, non-species-specific tool for examining host nutrition. Further analysis can include 2D NMR and spiking to identify additional metabolites, quantification of absolute concentrations using reference samples [PULCON method (73)], and metabolic pathway analysis [STOCSY (74)]. Despite differences in the nutrient environment revealed by NMR, *T. muris* counts and length did not differ between the dietary treatments. *T. muris* may feed on intestinal tissues and secretions (75, 76), rather than host blood or ingesta. *T. muris* also strongly affects, and is affected by, the host intestinal microbiome (76–78), so any changes in metabolic profile that alter microbiome composition could potentially have stronger effects on *T. muris* hatching, development, and survival. Future work could explore how host gut microbiome and metabolites effect *T. muris* infection success, and in turn, how they are affected by helminth infection.

### Dietary Quality and Infection Altered Aspects of Morphology and Behavior

Our data support the hypothesis that mice attempted to compensate for differences in the protein composition of the chow by altering their physiology and behavior. Mice spent 30% more time eating the LP diet per day, a significant time investment that could also come with increased predation risk in fully natural settings (79). While this investment could partially close the gap in protein acquisition between treatments, consumption would need to be 400% higher in the LP treatment to achieve similar protein levels to individuals on the HP diet. However, dietary protein does not affect host metabolic rate (16), so maintenance costs (**Figure 1**) are likely similar between treatments. Insects were also present in the enclosures and diet metabarcoding tools could be used in future studies to examine if, and to what degree, lab mice are able to utilize them as a food source. Albumin was reduced on the LP diet, revealing protein limitation within those hosts. Albumin's primary role is maintaining osmotic pressure in the blood (80), and, in wildlife, lower levels can reflect costs of reproduction (81) and indicate reduced survival probability in wildlife (37).

Interestingly, the choice of supplemental wild forage was marginally associated with differences in leptin concentrations. Leptin concentrations tended to be higher in animals that ate proportionally more plants in the legume family, Fabaceae, which includes clover (*Trifolium* sp.) and beggar's lice (*Desmodium* sp.). These plants have a higher protein content and are widely known to be good forage for livestock and wildlife (82), even increasing sheep weight gain by an average of 40% compared with feeding on ryegrass alone (83). Thus, behavioral compensation for immunological or infection costs may explain the trend toward higher consumption of these plants in infected mice. Conversely, lab mice that ate proportionally more plants in the families Oxalidaceae (*Oxalis* sp.; wood sorrel) and Violaceae (Viola sp.) tended to have low leptin levels. Mice consumed similar amounts of chow in the enclosures (g chow/g mouse) as they did in the laboratory setting, so wild plants probably do not represent a replacement food source for most individuals. We cannot quantify the *amount* of wild plant matter eaten by mice in the different treatments using DNA metabarcoding, but these emerging trends are suggestive of compensatory foraging behaviors worth further investigation.

Infected mice found physiological ways to compensate for the costs of *T. muris* infection, rather than following our alternate hypotheses of increased foraging or infection-induced inappetence. The increased cecum size of infected individuals could be a consequence of parasite manipulation to create more habitat space, but the increase in large intestine size with infection is more difficult to explain. Large intestines, emptied of contents and relative to body size, were over 10% heavier in infected mice. This additional weight was not due to the worms themselves, which were located in the cecum. In the average size mouse, this difference translates to a 15-mg increase in colon weight. Hosts could increase colonic tissues to enhance water balance, nutrient reabsorption, or house more commensal microbes to aid in digestion. Similarly, *H. polygyrus*, which resides in the small intestine, is associated with increased investment in small intestine tissues (67). An influx of immune cells or an increase in gut microbe communities could also contribute to these differences; this hypothesis could be examined histologically in future studies. However, given estimates of approximately 1 pg per *E. coli* cell (84) and 2.2 pg per lymphocyte (85), they cannot fully account for the increase in colon size in *T. muris*-infected mice.

To some degree, infected mice may also supplement their food intake with clovers*.* This trend was driven by a lack of beggar's lice in feces of uninfected mice and an absence of clover in LP-uninfected mice. Plants in the Fabaceae family, including clovers and beggar's lice, tend to be high in protein (82), and that could help hosts either resist infection by providing more resources for immune defense or tolerate infection by compensating for costs of infection and immune defense. A larger sample size, particularly of LP-uninfected mice, would help elucidate the efficacy and generality of this potential behavioral compensation mechanism. At this stage, however, we can conclude that infected mice gained less weight per minute of chow feeding than uninfected mice, despite increased large intestine size and use of potentially high-protein wild forage. This in turn suggests that their compensation for the costs of infection (e.g., parasite resource theft, elevation of immune defenses, tissue repair, etc.) was incomplete.

### CONCLUSION

Ultimately, biomedical and evolutionary immunology both aim to explain how resource costs of immunity and infection at the cellular level scale up to the whole organism embedded in its natural environment. Seeking such an explanation—and indeed bridging from biomedical to evolutionary immunology requires altered experimental designs that allow organisms to modify both their behavior and physiology in response to infection. An experimental approach is critical given the plausibility of alternate hypotheses (e.g., increased or reduced investment in immunity with a high-protein diet; increased or reduced foraging in response to parasite infection). Laboratory mice in semi-natural enclosures such as those studied here provide such an opportunity. The enclosures provided an environment with natural thermal regimes and space for activity (e.g., foraging, digging burrows) that could generate stronger energetic demands, and therefore stronger trade-offs than under typical laboratory conditions. Given the importance of gut microbiomes in *T. muris* infection (76–78) and in nutrition–immune interactions in general (19, 86, 87), the more diverse gut microbiomes generated by the enclosures (See text footnote 1) likely provide more realistic results than laboratory settings would. Additionally, this system shows potential for future studies of microbiome-helminth-diet interactions pertinent to the increasing rates of diseases linked to nutrition and immune dysregulation in developed nations (e.g., diabetes, obesity, autoimmune disorders). Our custom-built feeding monitoring system and dietary DNA metabarcoding allowed detection of compensatory feeding on chow and wild plants, respectively, that otherwise masked effects of infection.

By investigating lab mice, we were able to utilize a large suite of physiological and immunological tools, which proved useful given their variable responses to protein and infection. Much work remains to deploy such tools and new experimental designs to definitively dissect the mechanisms of the host–parasite interaction. In semi-natural enclosures, we observed high interindividual variation that reduces statistical power and, therefore, requires much larger sample sizes than traditional laboratory studies. For example, the relationship between infection and weight gain was just above the significance threshold (*p* = 0.052), but our power to detect an effect with total sample size of 80 individuals was only 0.45, far below the ideal power of 0.8. Thus, although we failed to detect an effect of infection on weight gain, we cannot conclude that *T. muris* does not affect mouse weight gain. Indeed, once we corrected for variation in weight gain due to time spent feeding, an effect of infection on weight gain was detectable (*p* = 0.042). Additionally, to prevent uncontrolled disease transmission and "contamination" of the enclosures with parasite eggs, we ended the experiment before *T. muris* could develop into adults, which may be a more energetically costly parasite life-stage. Longer term studies and those with trickle infections (i.e., small doses over time) will provide additional insight into this host–parasite interaction. Finally, variation in individual movement and thermoregulatory behavior is difficult to monitor in the enclosures, but may contribute to overall energy budgets and weight gain. With some weatherproofing alterations, remote activity monitoring systems such as those developed to study the activity of barn mice (88), plus temperature-sensing chips, could be a useful addition to such field enclosure studies.

Nonetheless, examining the interactions among diet, nutrients, immunity, and parasites in a realistic context revealed the central role of feeding behavior in infection outcomes and the complexity of interactions among environmental resources and within-host dynamics. Future experiments must therefore account for behavioral heterogeneity among individuals if we are to elucidate costs of parasitism and defense. Moreover, in the wild, altering feeding behavior in response to infection is a strategy available to individuals, but it may come with costs in terms of energy spent foraging, predation risk, and less time available for other behaviors (e.g., reproduction). Housing laboratory mice in outdoor enclosures provided new insights into the resource costs of immune defense to helminth infection and how hosts modify their feeding behavior to compensate for those costs.

### ETHICS STATEMENT

This research was conducted in accordance with animal care protocols approved by the Princeton University Animal Care and Use Committee (Protocol no. 1982-14).

### AUTHOR CONTRIBUTIONS

SB and AG designed the study, with intellectual contributions from CC and AL. SB, CH, and AG carried out the experiment and performed immune and condition measurements. QC and RG designed and built the feeding monitoring system. Diet metabarcoding was performed by CH and TK, and analyzed by SB and TK. NMR was performed and analyzed by IP. SB analyzed the data (unless otherwise noted) and wrote the manuscript. All coauthors provided editorial feedback.

### ACKNOWLEDGMENTS

We thank Sriveena Chittamuri, William Craigens, Justine Hamilton, Jackie Leung, Daniel Navarrete, Rebecca Neill, Hannah Priddy, Rob Pringle, Ross Pringle, and Alison Salamandran for assistance in the lab and field. P'ng Loke and Rowann Bowcutt kindly provided *T. muris* parasites. We also thank J. Leung for comments on the draft manuscript. Funding was provided by the Princeton University's Department of

### REFERENCES

Ecology and Evolutionary and Center for Health and Wellbeing, and National Science Foundation Research Experience for Undergraduates program (DBI—1358737). License for the SIMCA software was generously provided for student use by Umetrics (Umea, Sweden).

### SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at http://www.frontiersin.org/articles/10.3389/fimmu.2017.01914/ full#supplementary-material.

Figure S1 | Interferon gamma (IFNg), interleukin 10 (IL10), and interleukin 17 (IL17) concentrations were higher in infected mice than uninfected mice (Wilcoxon tests; IFNg: *W* = 506, *p* = 0.046, IL10: *W* = 483, *p* = 0.010, IL17: *W* = 441, *p* = 0.0036), but did not vary with diet (IFNg: *W* = 843, *p* = 0.67, IL10: *W* = 929, *p* = 0.13, IL17: *W* = 899, *p* = 0.28).

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

*Copyright © 2018 Budischak, Hansen, Caudron, Garnier, Kartzinel, Pelczer, Cressler, van Leeuwen and Graham. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*

# The immune and non-immune Pathways That Drive chronic gastrointestinal helminth Burdens in the Wild

*Simon A. Babayan1 \*, Wei Liu1 , Graham Hamilton2 , Elizabeth Kilbride1 , Evelyn C. Rynkiewicz3 , Melanie Clerc3 and Amy B. Pedersen3*

*<sup>1</sup> Institute of Biodiversity, Animal Health & Comparative Medicine, University of Glasgow, Glasgow, United Kingdom, 2Glasgow Polyomics, Glasgow, United Kingdom, 3 Institute of Evolutionary Biology and Centre for Immunity, Infection and Evolution, School of Biological Sciences, University of Edinburgh, Edinburgh, United Kingdom*

Parasitic helminths are extremely resilient in their ability to maintain chronic infection burdens despite (or maybe because of) their hosts' immune response. Explaining how parasites maintain these lifelong infections, identifying the protective immune mechanisms that regulate helminth infection burdens, and designing prophylactics and therapeutics that combat helminth infection, while preserving host health requires a far better understanding of how the immune system functions in natural habitats than we have at present. It is, therefore, necessary to complement mechanistic laboratory-based studies with studies on wild populations and their natural parasite communities. Unfortunately, the relative paucity of immunological tools for non-model species has held these types of studies back. Thankfully, recent progress in high-throughput 'omics platforms provide powerful and increasingly practical means for immunologists to move beyond traditional lab-based model organisms. Yet, assigning both metabolic and immune function to genes, transcripts, and proteins in novel species and assessing how they interact with other physiological and environmental factors requires identifying quantitative relationships between their expression and infection. Here, we used supervised machine learning to identify gene networks robustly associated with burdens of the gastrointestinal nematode *Heligmosomoides polygyrus* in its natural host, the wild wood mice *Apodemus sylvaticus*. Across 34 mice spanning two wild populations and across two different seasons, we found 17,639 transcripts that clustered in 131 weighted gene networks. These clusters robustly predicted *H. polygyrus* burden and included wellknown effector and regulatory immune genes, but also revealed a number of genes associated with the maintenance of tissue homeostasis and hematopoiesis that have so far received little attention. We then tested the effect of experimentally reducing helminth burdens through drug treatment on those putatively protective immune factors. Despite the near elimination of *H. polygyrus* worms, the treatment had surprisingly little effect on gene expression. Taken together, these results suggest that hosts balance tissue homeostasis and protective immunity, resulting in relatively stable immune and, consequently, parasitological profiles. In the future, applying our approach to larger numbers of samples from additional populations will help further increase our ability to detect the immune pathways that determine chronic gastrointestinal helminth burdens in the wild.

Keywords: wild immunology, *Apodemus sylvaticus*, transcriptome, machine learning applied to immunology, *Heligmosomoides polygyrus*, anthelminthics

#### *Edited by:*

*Lluis Tort, Universitat Autònoma de Barcelona, Spain*

#### *Reviewed by:*

*Katherine Buckley, George Washington University, United States Jorge Galindo-Villegas, Universidad de Murcia, Spain*

> *\*Correspondence: Simon A. Babayan*

<EMAIL_ADDRESS>
#### *Specialty section:*

*This article was submitted to Comparative Immunology, a section of the journal Frontiers in Immunology*

*Received: 07 October 2017 Accepted: 09 January 2018 Published: 05 February 2018*

#### *Citation:*

*Babayan SA, Liu W, Hamilton G, Kilbride E, Rynkiewicz EC, Clerc M and Pedersen AB (2018) The Immune and Non-Immune Pathways That Drive Chronic Gastrointestinal Helminth Burdens in the Wild. Front. Immunol. 9:56. doi: 10.3389/fimmu.2018.00056*

### INTRODUCTION

Chronic helminth infections challenge our understanding of how the immune system functions. In natural populations, while individuals are continually exposed to helminth parasites, there is substantial variability in infection burdens, with some consistently showing no active infection (1, 2). This suggests that individuals can potentially control infection. In addition, while anthelminthics are commonly employed to reduce helminth burdens, after drug clearance, infections typically return to their initial burdens (3–5), suggesting that the immune system does not or cannot easily acquire complete protection against helminths. Such observations have lead to the hypothesis that hosts must balance the costs of helminth infection with the immunopathology and/or protein-energy-related costs associated with eliminating these parasites (6–8). If this trade-off exists, it suggests that the host can regulate the intensity of its immune attack on the parasites, and, alternatively, that these parasites can avoid and suppress host immune responses (9, 10). While these ideas have received substantial theoretical and experimental support, the mechanisms that underlie chronic helminth infection dynamics "in the real world," and whether or how the balance of resistance and susceptibility might change over the lifetime of an individual, either naturally or in response to vaccination, remain very poorly understood. We suggest that addressing these questions is necessary for the development of more effective and sustainable treatment and immunization against parasitic helminths, which are the leading cause of productivity loss in livestock (11–13) and remain a substantial agent of poverty in the developing world (14, 15). It is, therefore, paramount to study hosts in their natural environment, outside of the controlled laboratory, where the data may be messy, but where relationships between hosts and their natural parasite communities are the result of millions of years of coevolution and more likely to be biomedically relevant.

Identifying protective immune mechanisms in inherently variable natural populations warrants relatively large sample sizes, the possibility to manipulate parasitological, physiological, and environmental parameters, and the ability to monitor individuals repeatedly over time. The wood mouse *Apodemus sylvaticus* has been extensively studied by disease biologists since the 1960s (16, 17) and is now increasingly used for immunological studies (18, 19). It offers many of the features that originally made small rodents attractive experimental models (i.e., affordable, easy to maintain and handle, prolific breeding, and short time to maturity) and adds features that make it an ecologically and epidemiologically sound model for mammalian, including human, immunology, and parasitology (i.e., high diversity and prevalence of parasites, good trapability, large population sizes, genetic relatedness to *Mus musculus*). Indeed, *A. sylvaticus* harbors a diverse and prevalent community of parasite and pathogens that closely resembles those found in larger mammals, including humans (20–24), and domestic animals (25). Importantly, unlike *Mus musculus* which is naturally infected by very few gastrointestinal helminth parasites (26, 27), *A. sylvaticus* is the natural host of the nematode *Heligmosomoides polygyrus*, which is routinely found in >50% of wild wood mice (28) and is closely related to *H. bakeri*

(29), a species extensively studied as a model of human gastrointestinal helminths that is known for its ability to suppress the immune system of its host (30). In this study, we therefore used *H. polygyrus* as a model for chronic endemic helminthiasis, and aimed to identify immune networks that regulate infection burden in a natural host–helminth system. Due to the preponderance of the laboratory mouse model in immunology, there are few, if any, reagents developed and optimized for non-model organisms (31). We, therefore, utilized a transcriptomics approach. While the genome itself provides invaluable information about immune resistance to infection (32–34), messenger RNA expression profiles provide a time- and context-sensitive picture of the immune system as it responds to antigenic stimuli. Moreover, unlike PCR and other candidate gene-driven approaches, transcriptomics allow discovery of unexpected correlates of immune protection and inclusion of physiological processes in the discovery of determinants of resistance to disease that would otherwise be overlooked in a purely immunology-focused approach.

In this study, we aimed (i) to identify gene networks that were either positively or negatively associated with *H. polygyrus* infection burden, (ii) to use those genes to identify immune pathways that promoted immunity to *H. polygyrus*, and (iii) to test whether those protective pathways were affected by sex and anthelminthic drug treatment. We used a newly published *A. sylvaticus* genome for the assembly and analysis of RNAseq samples generated from the spleens of 34 wild-caught wood mice from two distinct populations showing natural variation in helminth burdens, and then randomly treated a subset with anthelminthics to experimentally reduce their nematode infection burdens. Identifying transcripts predictive of infection burdens in a natural system, and maximizing the chances of these associations to generalize across populations requires detecting potentially weak signals among many variables, while simultaneously avoiding spurious variation (false positives). We addressed this challenge by reducing transcripts into co-expression gene networks, and identifying therein robust predictors of immunity to *H. polygyrus* using supervised machine learning—a class of statistical models that can integrate multiple variables that each carry weak signal into a stronger predictor by mapping them to a response variable, here *H. polygyrus* burdens. These types of supervised machine learning approaches hold promise for ecological and immunological studies (35, 36). We then applied a similar supervised approach to a narrower set of transcripts explicitly associated with immunity to assess how the immune system regulates its response to chronic helminth burdens in both male and female mice. To ensure that our models were not specific to a single population, time point, and sequencing platform, and thereby to increase the generalizability of our conclusions, we split samples spanning two wood mouse populations repeatedly at random into training sets on which predictive models were generated, and corresponding test sets on which the predictions of trained models were validated against observed values. Finally, given the widely reported propensity of anthelminthic-treated individuals (humans included) to return to their initial helminth burdens within weeks, we assessed whether the correlates of protection identified by our supervised learning approach were impacted by anthelminthic treatment.

### MATERIALS AND METHODS

### Ethics Statement

All procedures on animals were approved by the University of Glasgow ethics committee and the UK Home Office (PPL60/4572) and conducted in accordance with the Animals (Scientific Procedures) Act 1986.

### Wood Mouse Field Treatment and Sampling

To minimize false discovery and ensure that our findings could generalize to other wood mouse populations, we included samples from two geographically and temporally distinct *A. sylvaticus* populations at two woodland sites: Haddon Wood (N 53.16°, W 3.1°, hereafter referred to as HW) in north-west England in 2011, and Callendar Wood (N 55.99°, W 3.78° hereafter referred to as CW) in the Central Lowlands of Scotland in 2015.

Mice were trapped using Sherman live traps (H. B. Sherman 2 × 2.5 × 6.5 inch folding trap, Tallahassee, FL, USA) baited with grain, carrot, and bedding material. Two traps were placed every 10 m in a 70 m × 70 m (total 128 traps per grid). At first capture, each individual was microchipped (AVID microchips, Lewes, UK). In CW only, mice were allocated to one of two sexbalanced groups: anthelminthic drug-treated or control. Age and reproductive status, assessed by the body mass, color of the coat, position of testes, occlusion of the vagina, and visible signs of lactation and pregnancy were also recorded. The drug treatment removed gastrointestinal nematodes, which are typically the most abundant parasites in wood mice (4). Drugs were given as a single oral dose of 2 µl g<sup>−</sup><sup>1</sup> of body weight of a mixture of 9.4 mg kg<sup>−</sup><sup>1</sup> Ivermectin and 100 mg kg<sup>−</sup><sup>1</sup> Pyrantel (IVM + PYR), and controls consisted of oral delivery of 2 µl g<sup>−</sup><sup>1</sup> of body weight of water (H2O). Mice that were recaptured 14 ± 3 days post first capture were sacrificed on site. After cervical dislocation and exsanguination by cardiac puncture, the spleen of each individual was extracted and immediately transferred to a tube containing 4 ml RNA later, followed by whole removal of the intestine and storage in phosphate-buffered saline for subsequent dissection and parasite identification.

### Parasitology

*Heligmosomoides polygyrus* are ingested at their third larval (L3) stage, penetrate the submucosa of the small intestine within 24 h, and migrate to the muscularis externa, where they develop into L4-stage larvae. 8–10 days after infection, adult worms begin to emerge into the lumen of the intestine and attach to the intestinal villi within ~14 days, where they mate and release eggs. We measured adult *H. polygyrus* intestinal burdens of all sacrificed mice to assess immune resistance to infection. Worm burdens are preferred over fecal egg counts because they have much lower variability and likely reflect interactions with the host's immune system during larval development. While it was not possible to measure individual level differences in exposure to parasite infective L3s in the wild, we selected mice of similar ages within both sexes, and assumed lifetime exposure levels were similar across all selected mice. Worm burdens were compared between sites (CW or HW), sexes and treatment groups using negative binomial generalized linear models because of the significant leftskew coupled with high worm prevalence, over a zero-inflated distribution.

### RNA Sequencing, Assembly, Annotation, and Quantification

A transversal segment of each spleen was cut under sterile conditions, weighed, and processed following the RNeasy kit (Qiagen). RNA quantity and quality were assessed using a Tapestation (Agilent Technologies), and stored at −80 until RNA sequencing. All RNA samples were evaluated on a Bioanlayzer (Agilent Technologies) immediately prior to Poly-A selection of messenger RNA. Ten wood mice from Haddon Wood were sequenced on an Illumina Solexa 454 with 100 bp paired end reads by Edinburgh Genomics (2011). Twenty-four wood mice from Callendar wood were sequenced on an Illumina NextSeq 500 with 75 pb paired end reads by Glasgow Polyomics (2015). Raw reads generated by the sequencers were quality-checked using FastQC v. 0.11.5 (37) and sequences of low-quality and sequencer adapters trimmed using cutadapt v. 1.14 (38). The resulting trimmed reads were then quality-checked using FastQC and MultiQC (39) as above (see MultiQC report in Supplementary Data Sheet 1), and aligned to a reference transcriptome and quantified using Kallisto v. 0.43.1 (40). An *A. sylvaticus* transcriptome generated from a recently published genome (assembly ASM130590v1, https://www.ncbi. nlm.nih.gov/nuccore/LIPJ00000000.1) was used as the reference for assembly and read counts.

### Transcriptome Analysis Dimensionality Reduction

To select within the full transcriptome only the genes that had a potential role in resistance or tolerance to *H. polygyrus* across the full wood mouse transcriptome, we first grouped highly interconnected transcripts that may form distinct biological pathways by applying a weighted correlation network analysis (41) with the R package WGCNA (42). This entailed constructing gene networks (or modules) using a co-expression similarity measure defined as: