The right response at the right time: Exploring helminth immune modulation in sticklebacks by experimental coinfection

Abstract Parasites are one of the strongest selective agents in nature. They select for hosts that evolve counter‐adaptive strategies to cope with infection. Helminth parasites are special because they can modulate their hosts’ immune responses. This phenomenon is important in epidemiological contexts, where coinfections may be affected. How different types of hosts and helminths interact with each other is insufficiently investigated. We used the three‐spined stickleback (Gasterosteus aculeatus) – Schistocephalus solidus model to study mechanisms and temporal components of helminth immune modulation. Sticklebacks from two contrasting populations with either high resistance (HR) or low resistance (LR) against S. solidus, were individually exposed to S. solidus strains with characteristically high growth (HG) or low growth (LG) in G. aculeatus. We determined the susceptibility to another parasite, the eye fluke Diplostomum pseudospathaceum, and the expression of 23 key immune genes at three time points after S. solidus infection. D. pseudospathaceum infection rates and the gene expression responses depended on host and S. solidus type and changed over time. Whereas the effect of S. solidus type was not significant after three weeks, T regulatory responses and complement components were upregulated at later time points if hosts were infected with HG S. solidus. HR hosts showed a well orchestrated immune response, which was absent in LR hosts. Our results emphasize the role of regulatory T cells and the timing of specific immune responses during helminth infections. This study elucidates the importance to consider different coevolutionary trajectories and ecologies when studying host‐parasite interactions.


SI.1 Information on Diplostomum pseudospathaceum sampling and use
Limnea stagnalis were sampled by hand or using a small dip net. In the laboratory, each snail was subsequently rinsed with filtered lake water and individually placed in a plastic cup (Bioware 200ml,Huhtamaki) with filtered lake water. After two hours of direct light exposure, snails were screened for trematode infections by inspecting the shed cercariae in the water of the cup. The snails were kept in 16L aquaria at 18°C, with 16 hours of light per day. We used clone mixes from a pool of snail hosts in every infection round in order to overcome strong influences of D. pseudospathaceum genotype specificities: following a recovery period of at least two weeks post sampling, D. pseudospathaceum positive snails were individually placed in plastic cups with filtered lake water and exposed to direct light for 60 minutes. After verification of infection status and snail viability, 10 snails shedding the largest number of cercariae of the day were transferred to new plastic cups with fresh lake water and exposed to direct light for another 60 minutes. Cercariae from this supernatant were used to create a pool of D. pseudospathaceum cercariae of similar age.

SI.2 Further information on reverse transcription
The reverse transcription protocol was modified by using 0.2 µL Qiagen RNAse inhibitor (instead of 1 µL). The manufacturer ensured that 0.2 µL is sufficient due to differences in effective inhibitor concentrations.

SI.3 Further information on direct sequencing
PCR conditions were the same in all sequencing attempts. All PCR products were checked on a gel for the right size and amplification specificity. 5 µL aliquot of the completed PCR reaction were mixed thoroughly. 2 µL of Illustra ExoStar 1 Step were added to the reaction mix and incubated at 37°C for 15 minutes. Incubation at 80°C for 15 minutes inactivated the enzymes. Afterwards, the cycle sequencing was prepared as follows:  Toll-like receptor 2; Germline-encoded pattern-recognition receptor (Zhu et al., 2013;Brunner et al., 2017) Amplification efficiency of primer product was not within acceptable range p40phox Component of NADPH oxidases (Stutz et al., 2015) Product sequencing revealed amplification of unspecific target vegfa1

SI.4 Further information on gene expression analyses
Stimulates macrophage and monocyte migration (Brunner et al., 2017) Unspecific primer products, ambiguous PCR products ly75 Reduces B-lymphocyte proliferation (Brunner et al., 2017) Unspecific primer products, ambiguous PCR products cmip Signaling protein in Th2 pathway (Robertson et al., 2015) Unspecific primer products, ambiguous PCR products SI.5 Gene expression targets, gene references and primer sequences Promotor of granulocyte and neutrophil migration, required for activation of the innate immune response (Leemans et al., 2004;Rhodes et al., 2009)

tnfr1
Tumor necrosis factor receptor 1; functions in regulation of inflammation, mediates cellular apoptosis and differentiation (Zhu et al., 2013;Brunner et al., 2017) ENSGACT00000013502 AACTACTACAGAGCCAAGGGCAAG ACGGCACTCAGCGGTACAATTC ENSGACT00000011232 CTCAGCCACAGTTCCAACCGTTC GTCGGATGTTCTGGACCTCGAGT tcr-β T-cell receptor β-chain; function in binding of MHC-peptide ligands to initiate adaptive immune response (Yanagi et al., 1984;Smith-Garvin et al., 2009;Stutz et al., 2015) ENSGACT00000016457 GAGGGCAAAAACTTCACCTG TAGGAGAATCTGGCCGTTTG tgf-β Transforming growth factor β; cytokine with functions in cell growth, migration, differentiation and proliferation of T and B-cells (Zhu et al., 2013;Robertson et al., 2015) ENSGACT00000016968 TCCCGCTTCGTCACCAACCA ACGTCTGTCTGGCCACATTCAC The following was pipetted into each well of a 96 well plate --3.9 µL Sample Pre Mix --3.1 µL sample (preamplified) --vortexed 20s, spun down 30s Assay Pre-mix -prepared in a 1.5 ml tube: total 665.3 µL (for 96 reactions, includes overage) -369.6 µL 2X Assay loading Reagent (Fluidigm, PN 85000736) -295.7 µL low TE buffer The following was pipetted into each well of a 96 well plate (7 µL per well) -6.3 µL Assay Pre-mix -0.35 µL from each of the 100µM primers (fwd and rev ) or 0.7 µL from the mix After priming of the chip, Sample Pre Mix and Assay pre Mix were loaded according to the manufacturer's instructions and the chip was run under cycler protocol: "GE Fast 96x96 PCR+Melt v2". We tested effects on S. solidus infection rates by using a GLMM with the origin of the host and the origin of the parasite as well as their interaction as fixed effects and fish family as random term. The model did not differ significantly from the Nullmodel (likelihood ratio test: Χ 2 3 = 4.2365, p = 0.237).

SI.9 Further information on S. solidus growth and parasite index
S. solidus growth differed significantly between the two parasite populations. HG parasites grew faster and larger than LG parasites (Type III Wald chisquare tests: H:P:T three-way interaction: Χ 2 4 = 24.8413, p < 0.0001).  week Weight of the worm HG worm (LR fish) HG worm (HR fish)

SI.10 Further information on D. pseudospathaceum infection rates
LG worm (LR fish) LG worm (HR fish)

SI.11 The effect of S. solidus weight on D. pseudospathaceum infection rates
Using S. solidus weight as a covariate in the statistical model did not improve the model fit in week 3, but did so at later time points, namely for data from LR fish in week 6 (likelihood ratio test: Χ 2 2 = 10.01, p = 0.0067) and data from both fish origins in week 9 (likelihood ratio test: Χ 2 2 = 13.37, p = 0.0013). The model fit for HR data of week 6 was not improved (likelihood ratio test: Χ 2 2 = 4.82, p = 0.0897). Due to very large eigenvalues, we z-transformed the weight of the worm in week 6 and week 9.

SI.12 Further information on host condition and immunological parameters
The condition factor (CF) differed significantly between host populations (F 1,4 = 25.027, p = 0.0075) and according to an interaction between time point and treatment (F 10,170 = 2.543, p = 0.007). FDR-corrected post hoc comparisons confirmed significantly higher condition of HR than LR hosts at all time points (LMMs; each p < 0.0001). Treatment had no significant influence; the CF increased between week 3 and week 9 in D. pseudospathaceum infected HR fish. The hepatosomatic index (HSI) was significantly affected by an interaction between treatment and time point (F 10,170 = 4.102, p < 0.0001). LR controls had higher HSIs than HR controls in week 9; the HSI increased between week 3 and 9 in LR controls (LMMs; each p < 0.001). In week 9, LG infection was associated with a smaller HSI than in controls in LR fish; infection with D. pseudospathaceum correlated with significantly higher HSI in comparison to co-infection with HG S. solidus in LR sticklebacks (LMMs; each p < 0.001). Splenosomatic indices (SSI) and head kidney indices (HKI) were not affected by experimental factors. The statistical models were based on log10-transformed calibrated normalized relative quantities (CNRQ values). The weight of the fish was included as covariate to account for size related effects. Non-parametric permutational multivariate analyses of variance (PERMANOVA) were calculated on Euclidean distances and 10,000 permutations that were constrained within fish family. PERMANOVA results were FDR corrected. If significant (marked in bold letters), single genes were analyzed with linear mixed models (LMMs). Statistics for differences between host types or interactions are mentioned whenever significant. T: time point (week 3, 6, or 9); P: parasite type (low growth, LG; high growth, HG); H: host type (low resistance, LR; high resistance, HR); all genes: data from all 23 genes; innate: 11 genes (cd97, csf3r, il-1β, marco, mif1, mst1ra, nkef-β, p22 phox , saal1, sla1, tnfr1); adaptive: nine genes (stat4,cd83,igm,stat6,foxp3,mhcII,; complement: three genes (cfb, c7, c9); Th1: two genes (stat4, tnfr1), Th2 covers three genes (stat6, cd83, igm); Treg covers three genes (il-16, foxp3, tgf-β). The statistical models were based on log10-transformed calibrated normalized relative quantities (CNRQ values). The weight of the fish was included as covariate to account for size related effects. Data from genes from significantly differentially expressed functional gene groups was analyzed with linear mixed models (LMMs; function lme() from nlme) and analyses of variance (ANOVAs). Conditional pseudo R 2 values (Nakagawa & Schielzeth, 2013;Johnson, 2014) were calculated with sem.model.fits() from piecewiseSEM (Lefcheck, 2015). No gene was significantly differentially expressed after FDR correction.. The statistical models were based on log10-transformed calibrated normalized relative quantities (CNRQ values). The weight of the fish was included as covariate to account for size related effects. Data from genes from significantly differentially expressed functional gene groups was analyzed with linear mixed models (LMMs; function lme() from nlme) and analyses of variance (ANOVAs). Conditional pseudo R 2 values (Nakagawa & Schielzeth, 2013;Johnson, 2014) were calculated with sem.model.fits() from piecewiseSEM (Lefcheck, 2015). Differentially expressed genes are marked in bold letters if significant after FDR correction. The statistical models were based on log10-transformed calibrated normalized relative quantities (CNRQ values). The weight of the fish was included as covariate to account for size related effects. Data from genes from significantly differentially expressed functional gene groups was analyzed with linear mixed models (LMMs; function lme() from nlme) and analyses of variance (ANOVAs). Conditional pseudo R 2 values (Nakagawa & Schielzeth, 2013;Johnson, 2014) were calculated with sem.model.fits() from piecewiseSEM (Lefcheck, 2015). Differentially expressed genes are marked in bold letters if significant after FDR correction. The statistical models were based on log10-transformed calibrated normalized relative quantities (CNRQ values). Nonparametric permutational multivariate analyses of variance (PERMANOVA) were calculated on Euclidean distances and 10,000 permutations that were constrained within fish family. The weight of the fish was included as covariate to account for size related effects. D effect: effect of D. pseudospathaceum infection. In this case, gene expression was only affected by an interaction between host type and time. The statistical models were based on log10-transformed calibrated normalized relative quantities (CNRQ values). The weight of the fish was included as covariate to account for size related effects. Non-parametric permutational multivariate analyses of variance (PERMANOVA) were calculated on Euclidean distances and 10,000 permutations that were constrained within fish family. PERMANOVA results were FDR corrected. If significant (marked in bold letters), single genes were analyzed with linear mixed models (LMMs). Statistics for differences between host types or interactions are mentioned whenever significant. The statistical models were based on log10-transformed calibrated normalized relative quantities (CNRQ values). The weight of the fish was included as covariate to account for size related effects. Data from genes from significantly differentially expressed functional gene groups was analyzed with linear mixed models (LMMs; function lme() from nlme) and analyses of variance (ANOVAs). Conditional pseudo R 2 values (Nakagawa & Schielzeth, 2013;Johnson, 2014) were calculated with sem.model.fits() from piecewiseSEM (Lefcheck, 2015). Differentially expressed genes are marked in bold letters if significant after FDR correction. The statistical models were based on log10-transformed calibrated normalized relative quantities (CNRQ values). The weight of the fish was included as covariate to account for size related effects. Non-parametric permutational multivariate analyses of variance (PERMANOVA) were calculated on Euclidean distances and 10,000 permutations that were constrained within fish family. Statistics for differences between host types or interactions are mentioned whenever significant. PERMANOVA results were FDR corrected. In this case, no result remained significant after FDR correction. The statistical models were based on log10-transformed calibrated normalized relative quantities (CNRQ values). The weight of the fish was included as covariate to account for size related effects. Non-parametric permutational multivariate analyses of variance (PERMANOVA) were calculated on Euclidean distances and 10,000 permutations that were constrained within fish family. Statistics for differences between host types are mentioned when significant and infection effects where then tested for HR and LR host types separately. PERMANOVA results were FDR corrected. In this case, no results remained significant after FDR correction. D: Diplostomum pseudospathaceum infection.