Systems biology analyses of the dynamic host response to Toxoplasma gondii infection in a murine model

Summary

Toxoplasmosis affects a third of the global population and is of particular concern for immunologically compromised individuals. Toxoplasmosis induces host physiological events ranging from immunological to metabolic responses across multiple biological compartments. To understand the sequence of host responses during acute and chronic Toxoplasma gondii infection, eight male BALB/c mice were infected with 2000 T. gondii ME49 tachyzoites with a further eight uninfected mice used as controls. Plasma cytokines status, urinary metabolic profiling and fecal microbial profiles were characterized to monitor temporal variation related to T. gondii infection. The results showed elevated serum interferon-γ (IFN-γ), interleukin-12p40 and  necrosis factor-α during acute phase of infection with concomitant perturbation in host energy metabolism and host-gut microbiome co-metabolism of phenolics and a shift in microbial composition. However, the differences were less pronounced during the putative chronic phase of infection with elevated IFN-γ, differences in urinary N-acetyls and O-acetyls of glycoproteins with no shift in gut microbial composition. Structural equation modelling on the current data showed host immune responses as the main driver for changes observed in urinary metabolites and gut microbial composition. Such an approach can be applied to other models of infectious diseases to aid understanding of host–pathogen interactions and potential biomarker discovery.

INTRODUCTION

Toxoplasmosis, one of the most common parasitic diseases worldwide, is caused by the protozoan Toxoplasma gondii. It is a parasite with no specific host and has more than one obligatory host in its life cycle (Tenter et al. 2000). The only known definitive host of T. gondii is the Felidae family mainly wild and domestic cats. T. gondii has several potential routes of transmission within or between host species that can be either vertical (transplacentally via tachyzoites) or horizontal (orally via ingestion of T. gondii at various stages of its life cycle) transmission (Tenter et al. 2000). Humans are usually infected through ingestion of oocysts from contaminated food or water or by contact with cats. Toxoplasma gondii is prevalent in most parts of the world with an estimated one third of the global human population infected (Montoya and Liesenfeld, 2004). Toxoplasma gondii infection may cause miscarriage, stillbirth or congenital disease in its hosts, thus it is of importance to the medical and veterinary world.

Metabolic profiling approaches were first introduced into parasitology in an initial study on Schistosoma mansoni-infected mice (Wang et al. 2004). Since then, metabolic profiling strategy has been used to investigate other rodent host–parasite models including Schistosoma japonicum-hamster (Wang et al. 2006), and Trypanosoma brucei brucei-mouse (Wang et al. 2008) with the aim of characterizing the metabolic responses of hosts to these parasites at the molecular level. Although such single matrix/platform approach allowed for comprehensive characterization of host-parasite interaction from the metabolic perspective, a combination of different sample matrices, analysed using different platforms, will enable a more holistic understanding of such complex biological interactions involving multiple levels of biological complexities. Correlating immune response to an invading organism as well as to changes in the host metabolic profile can give new insights into immunological processes (Lamour et al. 2015) and may contribute to new developments in early diagnosis, therapy and immune-prevention. Monitoring of changed biochemical pathways will help to discover the factors responsible for either disease susceptibility (Lamour et al. 2015) or resistance (Canuto et al. 2014).

Cell-mediated immune (CMI) response is of particular importance in infection with intra-cellular parasites and is the major mechanism of resistance against T. gondii (Wong and Remington, 1993). In addition, resistance in mice after T. gondii infection was dependent on both CD4⁺ and CD8⁺ T lymphocytes (Hakim et al. 1991). In vivo treatment with either recombinant interleukin (IL)-12, interferon-γ (IFN-γ), necrosis factor-α (TNF-α), IL-1α and IL-1β, all cytokines of the Th1 subset, has had a protective effect against T. gondii infection (Gazzinelli et al. 1993, 1994). Synthesis of macrophage IL-10 has been shown by several studies to enable the establishment of chronic infections by different pathogens, including T. gondii, by avoiding the induction of an excessively strong CMI and thus allowing the parasite to persist in the vertebrate host (Fiorentino et al. 1991). Furthermore, T. gondii infection results in significant mucosal response with IgM detected 1–2 months after infection followed by IgA and a vigorous IgG antibody response was observed in patients with acute toxoplasmosis (Turunen et al. 1983). Despite enormous progress in toxoplasmosis research and understanding the genome (Boyle et al. 2006) and proteome (Zhou et al. 2004) of T. gondii, more work is needed to ascertain the complex developmental changes within different hosts once T. gondii is established. The host metabolic response to T. gondii infection and the ensuing systematic and metagenomic responses to parasitic nutrient demand and survival will require systems biology approaches to delineate the effect and ascertain the role of each biological compartment to parasitic insult. Correlation between the host metabolic, metagenomic and immune responses to parasitic infection may lead to new insights into host–parasite interactions (Lamour et al. 2015). A plethora of cofactors involved in signal transduction and epigenetic changes of T. gondii are essentially by-products of metabolic pathways such as amino acids and lipid pathways (Kim and Weiss, 2008).

The current study was conducted to characterize the complex host–microbiome–parasite interactions of toxoplasmosis in an experimental murine model. Here a systems biology approach was used to investigate the acute and chronic phases of T. gondii infection by integrating immune, metabolic and gut microbial matrices of the host particularly selected peripheral cytokines, metabolic profiles from urine, and fecal bacterial composition.

MATERIALS AND METHODS

Ethics

Animal care and handling procedures used in this study followed the guidelines and requirements of local and national authorities for the care and use of animals for scientific purposes in Malaysia. The study was approved by the Joint Committee on Research and Ethics, International Medical University, Malaysia (IMU-R-093/2012).

Animal handling and sampling

The toxoplasmosis murine model used in this study was established according to the protocol used by (Chew et al. 2012a , b ). All animals used in this study were male BALB/c mice (4–5 weeks of age) purchased from a university-approved experimental animal supplier (Selangor, Malaysia). A total of 16 animals (eight controls and eight infected) were housed in groups of 4 in macrolon cages and maintained at controlled environmental conditions (temperature: 26 °C, relative humidity: 60–70%, light/dark cycle: 12/12 h) with free access to water and standard rodent chow. Mice were acclimatized for 14 days prior to the start of experiment, in order to minimize potential stress-related variation in biological profiles. For the experimental group, eight mice were infected intraperitoneally with 1 mL of saline solution containing 2000 T. gondii tachyzoites (ME49 strain) each whilst the remaining eight control mice were each given 1 mL of saline solution intraperitoneally.

Urine, blood and fecal samples were collected one day before infection (pre-infection) and on days 7, 14, 21, 28, 35 and 42 post-infection. Sample collection took place between 0800 and 1100 h to minimize diurnal variation in concentrations of biofluids. At least 50 µL of urine and eight fecal pellets were collected from each mouse into Petri dishes by gently rubbing the abdomen. Sample which do not achieved targeted volume/count was excluded from the subsequent analyses. To avoid cross-contamination, new gloves and Petri dish was used for each separate animal. Urine and fecal samples were collected for one-dimensional (1D) proton nuclear magnetic resonance spectroscopy (¹H NMR)-based metabolic profiling and terminal restriction fragment length polymorphism (TRFLP) analyses, respectively. Small amounts of tail blood samples (50–80 µL) were collected for cytokine profiling. All samples were collected into 1·5 mL microcentrifuge tubes, and immediately stored at −80 °C. Animals were weighed at each time point, using a weighing balance (Precisa, Switzerland). Cardiac puncture was performed pre-terminally on all animals 42 days post-infection (dpi) under deep anaesthesia for collection of blood. The collected blood samples were left at room temperature for clotting for at least 30 min and subjected to centrifugation at 1000  g for 10 min. The sera was carefully separated from the red blood cells using a micropipette and stored at −80 °C until further use. Brain and other organs were harvested subsequently at autopsy.

Histopathology

Mouse brain tissues harvested on day 42 were fixed in 10% neutral buffered-formalin and embedded in paraffin wax. The 4 µm thick sections were stained with haematoxylin and eosin stain (Sigma Aldrich, USA) for subsequent histopathological assessment and toxoplasma cyst identification.

Cytokine multiplex quantification

A ten-plex cytokine immunoassay customized for IFN-γ, tumour necrosis factor-α (TNF-α) and a panel of IL; IL-1β, IL-2, IL-4, IL-5, IL-6, IL-10, IL-12p40, IL-17 was carried out using the Mouse Cytokine/Chemokine Magnetic Bead Panel Kit in a 96-well plate # MCYTOMAG-70 K (Millipore, USA) according to manufacturer's instructions, using a hand-held magnet for the washing steps. The plate was read using the Luminex 200TM Milliplex Analyser (Luminex, USA) with the Luminex Software xponentv3·1 (Luminex, USA) for data acquisition. The Median Fluorescent Intensity (MFI) data were analysed using a five-parameter logistic or spline curve-fitting method for calculating cytokine concentrations in the samples. The values were multiplied by two which was the dilution factor used for sample dilution for calculation of concentrations. The experiment was carried out in duplicates for each sample.

¹H NMR spectroscopy of urine

Urine samples were prepared by mixing 20 µL of urine with 30 µL of a phosphate buffer containing 90% deuterated water (D₂O) and 0·25 mm 3-trimethylsilyl-1-[2,2,3,3-²H₄] propionate (TSP) and left to stand for 10 min. The resulting mixtures were then transferred into 5 mm NMR microbore tubes with 1·7 mm stem (Norell, USA).

A standard 1D ¹H NMR spectrum was acquired for each sample with a pulse [recycle delay (RD)-90°-t ₁-90°-t _m-90°-acquire free induction decay (FID)] on a Bruker AVIII 600 MHz spectrometer (Bruker Biospin, Fallenden, Switzerland) with a 5 mm BB(F)O broadband probe operating at 600·13 MHz (ambient probe temperature 27 °C). Samples were automatically delivered to the spectrometer by a Bruker SampleCase^™. The field frequency was locked on D₂O solvent. The water peak was suppressed by gradient water pre-saturation during the relaxation delay (RD) of 4 s and mixing time (t _m), of 0·01 s. The 90° pulse length was adjusted to ~10 µs and an acquisition time of 2·65 s was used. In total, 32 transients were collected into 64 K data points using a spectral width of 20 ppm. An exponential line broadening function of 0·3 Hz was applied to the free induction decay prior to Fourier transformation.

Fecal DNA extraction and TRFLP analyses

DNA was extracted using QiAmp DNA stool mini kit (Qiagen, USA) according to the manufacturer's protocol. The TRFLP was carried out as described in (Chong et al. 2012) with minor modifications. Briefly, PCR was conducted using universal primers (27F and 1492R) targeting bacterial 16S rRNA gene region. Both primers were tagged with the fluorescence dye phosphoramidite fluorochrome 5-carboxyfluorescein (FAM) and 6-carboxy-hexachlorofluorescein (HEX) at 5′, respectively. The amplicons were digested using a four-base cutter MSP-1 before subjected to electrophoretic separation using ABI 3100 and ABI 3730XL genetic analysers (Applied Biosystems, USA). ROX labeled GeneScan 500 control was used as size standard. The resulting fragment profiles were scored, aligned and noise filtered using the web-based programme T-REX (Culman et al. 2009). Noise filtered procedure was carried out by removing peaks, which are lower than the overall standard deviation. In addition, alignment was carried out by binning the peaks into TRFs with the clustering threshold of 0·5 bp, starting from the smallest peak among all the profiles (Smith et al. 2005).

Data processing and multivariate statistical analysis

For metabonomics analysis, ¹H NMR spectra of urine samples were manually phased and baseline-corrected using Bruker TopSpin (Version 3.1, Bruker Biospin, Fallenden, Switzerland). All spectra were referenced to the TSP resonance at δ 0·00. The complete spectra (δ 0·0–10·0) were digitized into 7 K data points using an in-house developed MATLAB script (O. Cloarec, Imperial College London). The region containing the water resonance (δ 4·5–6·5) was removed from each spectrum to eliminate baseline effects of imperfect water suppression. In addition, the regions δ 0·0–0·5 and δ 9·2–10·0 containing only noise were removed. For each spectrum, normalization to the total sum of the residual spectrum was carried out prior to pattern recognition analyses followed by scaling of the data to unit variance. Both unsupervised (principal component analysis, PCA) and supervised multivariate data analysis methods (projection to latent structures-discriminant analysis, PLS-DA and orthogonal partial least squares discriminant analysis, O-PLS-DA (Trygg and Wold, 2002; Cloarec et al. 2005)) were employed to visualize and interpret experimental differences. The statistical significance and validity of subsequent results were calculated using a permutation test (number of permutations = 1 000). The covariance plot was used to aid interpretation of the significance of each metabolite from the permutation tests. The colours projected onto the spectrum indicate the significance of the metabolites with blue indicating no significant difference at P > 0·05 confidence levels and red indicating high significant difference at P < 0·05. In this study, only metabolites with P < 0·01 will be considered significant.

The TRFLP data obtained from T-REX was exported into PERMANOVA + add-on of the PRIMER6 multivariate data analysis package (Plymouth Marine Laboratory, UK) for analyses. Alpha-diversity indices including Shannon diversity index (H′) and Pielou's measure of species evenness (J′) were calculated using the DIVERSE option. For multivariate-analyses, the data was pre-treated with unit variance normalization procedure before generating a Euclidean distance matrix, which was used to perform permutational multivariate analysis of variance (PERMANOVA) and PCA. The PERMANOVA was carried out to test for gut microbial compositional difference between the groupings, while the PCA was used to visualize the distribution of the sample across the sampling points. The significance of the PERMANOVA was calculated based on 999 permutations and corrected using Monte Carlo correction.

Canonical analysis of principal coordinates (CAP) is a ‘constrained’ ordination method which uses a priori hypotheses for dimension selection (Anderson and Willis, 2003). In this study, CAP was use to illustrate dimensions with strongest separation in relation to infection and time points.

Structural equation modelling (S.E.M.)

s.e.m. was carried out using ‘Lavaan’ and visualized with ‘semPlot’ package under R programme.(Epskamp, 2015)^, (Rosseel, 2012) A total of three models were built to represent the hypothesized relationships among the measured variables under each category which make up of cytokines (immunology), NMR profiles (metabolic composition) and TRFLP profiles (gut microbial composition). Difference in the model fit was compared using chi-square different test. Model showing lowest reading for Chi-square value, Akaike information criterion (AIC) and Bayesian information criterion (BIC) scores were selected as the best-fitted model.

Distance-based linear modelling (DISTLM)

DISTLM were carried out as described in (Yap et al. 2015) Briefly, the parameters showing the most consistent responses to the corresponding overall profiles (NMR, TRFLP or cytokines) was selected using stepwise selection under the second-order bias-corrected AIC (Yanagihara et al. 2011).

RESULTS

Histology and physiological monitoring

The mean bodyweight between control and infected mice was significant across time point (F = 11·69, P < 0·01) but insignificant across treatments (i.e. control vs infection, F = 1·46, P = 0·20) (see Supplementary Materials Fig. S1, available from http://journals.cambridge.org/PAO). Putative chronic infection was established by day 42 as observed from the haematoxylin and eosin stained paraffinized brain tissues (Fig. 1). No brain cyst was observed from brain tissue collected from day 14.

Fig. 1. Brain histopathology of control (A) vs infected (B) BALB/c mice 42 dpi. Toxoplasma gondii cysts (arrows indicate bradyzoite tissue cysts) were observed in the brains of BALB/c mice 42 dpi (B). Footnote: The data are representative of mice in the whole experiments [n = 16 (eight controls; eight infected)]. All tissues were H and E-stained and viewed at ×400 magnification.

Characterization of host immune responses to parasite infection

Pattern in overall host immune responses derived from cytokine multiplex (see Supplementary Materials Table S1, available from http://journals.cambridge.org/PAO) was analysed using CAP plot (constrained analysis). In brief, the procedure selects axes showing greatest discrimination against treatments (i.e. control vs infected) and time-points (days 0, 1, 7, 14, 21, 28, 36 and 42). A strong separation between control and infected animals was apparent in the first component (CAP1 – black – control vs red – infected) (P = 0·001, Q ² = 46·3%) (Fig. 2A). The demarcation between the control and infected animals was mainly contributed by the upregulation of IL-1β in the former while upregulation of IL-12p40 was observed during the acute infection phase and overall upregulation of TNF-α and IFN-γ throughout the infection (vector overlay – Fig. 2A). Further, among the infected group, samples from days 7 and 14 formed a distinct cluster at the bottom of the second component (CAP2). A follow-up PERMANOVA analysis exhibited a significant infection effect (pseudo-F = 9·43, P = 0·001) and significant interaction between infection and time-points (pseudo-F = 2·115, P = 0·022) (Table 1). This suggested that the time-course immunological response for both groups were significantly distinct to each other. However, no significant difference was detected when comparing the pooled cytokines profiles (both control and infected groups) between different days (pseudo-F = 1·624, P = 0·087).

Fig. 2. (A) Canonical analysis of principal coordinates (predictive power, Q ²: 46·3%) scores and corresponding loadings (blue subplot) plot of cytokine measurements from infected (red) and control animals (black); (B) PCA time trajectory plots derived from serum cytokine data of infected animals showing differentiation of pre-infection (PI) from post infection at each time point (days 7, 14, 21, 28, 36 and 42). (C) PCA scores plot derived from urinary ¹H NMR spectra of male BALB/c mice from infected and control animals. (D) Subsequent PCA mean trajectory plot derived from the urinary ¹H NMR spectra of infected animals only. The error bar indicates the standard error. Key: , PI; , Day 7; , Day 14; , Day 21; , Day 28; , Day 36; , Day 42.

Table 1. Permutational multivariate analysis of variance (PERMANOVA) results based on 1000 times permutation for immunological, metabolic and microbial profiles

df, degree of freedom; SS, sum of square; MS, mean square.

To understand the dynamics of host cytokine responses to T. gondii infection, a PCA time trajectory plot was constructed for the infected mice (Fig. 2B). The principal component time trajectory plot showed clear partitioning between samples from days 7 and 14 and that of pre-infection. However, samples from day 21 onwards converged with the cytokine levels detected during pre-infection.

Characterization of host metabolic responses to parasite infection

Initial analysis of the urinary NMR data using CAP showed clear separation between the control and infected animals (Fig. 2C) with the pre-infected animals overlapping with the controls along the first principal component. Subsequently, principal component mean trajectory analysis was performed on the urinary data derived from the infected groups to identify key time points where infection trajectory changes over time. The resulting mean trajectory scores plot (Fig. 2D – infected group) showed that, upon infection, the metabolic trajectory moved along the first principal component from pre-infection to day 7. The metabolic trajectory moved further along the first principal component on day 14. However, by day 21, the metabolic trajectory changed direction along the first principal component, whereby the metabolic composition of samples collected from days 28, 36 and 42 clustered in the same region. The metabolic composition for day 42 was found to ordinate close to the pre-infection (Fig. 2D). Nevertheless, the corresponding PLS-DA scores plot derived from the urine NMR data of days 14 and 42 samples (Fig. 3A and B, respectively) showed clear separation between the infected and control animals indicating that the infected animals never recovered 42 dpi. The model predictive power, Q ², was 95 and 77·3% respectively, indicating the statistical validity of the model.

Fig. 3. PLS-DA scores plot derived from urinary ¹H NMR spectra of infected and control animals on (A) day 14, and B) day 42. The corresponding covariance plots showing the colour-coded significance of urinary metabolite profiles calculated using permutation test between infected and control animals on (C) day 14, and (D) day 42. Key: 2OG, 2-oxoglutarate; DMG, dimethylglycine; GAA, guanodinoacetate, NAG, N-acetyls of glycoproteins; OAG, O-acetyls of glycoproteins; PAG, phenylacetylglycine; TMAO, trimethylamine N-oxide

Significant metabolites partitioning according to both treatments and time-points were established based on PERMANOVA (pseudo-F = 1·66 and 3·44 respectively, P = 0·001). In addition, contrasting temporal response pertinent to the status of infection (control vs infected) was also detected (interaction – pseudo-F = 1·53, P = 0·005) (Table 1). Statistical significance of the urinary metabolites was calculated using permutation test (number of permutations = 1000). Only metabolites with a P-value of 0·01 and below were selected and summarized in Table 2. Figure 3C and D show the corresponding covariance plots of control vs infected animals derived from urinary samples obtained from days 14 and 42. Urine collected from mice on day 14 showed elevated α-hydroxy N-valerate, 2-oxoglutarate, β-aminoisobutyrate, citrate, acetylcholine, fumarate, guanidinoacetate, hippurate, methionine and N-methyl nicotinate in the control group, whereas the infected group showed elevated urinary creatine, dimethylglycine, N-acetyl glutamate, phenyalacetylglycine, taurine and trimethylamine N-oxide. In contrast, changes in urinary metabolites were less pronounced on day 42 with the control group showing higher levels of urinary α-hydroxy N-valerate, acetate, guanidinoacetate, N-acetyls- and O-acetyls of glycoproteins, whereas the infected group showed higher levels of urinary fumarate, glycine, hippurate, phenyalacetylglycine and trimethylamine N-oxide.

Table 2. ¹H NMR-derived metabolites that differ significantly between control and infected animals on days 14 and 42

All metabolites were significant at the level of P < 0·01 calculated using permutation test.

s, singlet; sb, singlet broad; d, doublet; t, triplet; m, multiplets.

*, tentatively assigned as.

+, relatively higher in the control group.

–, relatively lower in the control group.

Characterization of host fecal microbiota composition in response to parasite infection

Fecal bacterial richness (Shannon diversity index) and evenness (Pielou's evenness) of the experimental animals were inferred using TRFLP (see Supplementary Materials Table S2, available from http://journals.cambridge.org/PAO). Sample with more TRFs is generally considered ‘richer’ than sample with low TRFs count. On the other hand, community that harbour TRFs with similar intensity is considered ‘more even’ than community containing large number of low-intensity peaks but dominated by only a few high dominant peaks. Based on Shannon Diversity index (H′), no significant difference in bacterial richness was detected (T = 1·88, P = 0·064) between the control and infected animals (H′_control = 3·30 ± 0·92, H’_infected = 3·06 ± 0·92). However, significant lower Pielou's evenness (J′) was detected (T = 2·67, P = 0·009) in the infected arms in comparison to the controls (J′_infected 0·77 ± 0·01, J′_control = 0·80 ± 0·01). When the gut microbial distribution pattern was ordinated on the CAP plot (Figure 4A), partitioning in accordance to infection and control groups was observed on the vertical axis (CAP2 – P = 0·01, Q ² = 30·4%). Further, among the infection arm, separation of samples collected on day 14 was apparent. However, the same pattern was not detectable under the PCA trajectory plot for the infection group (Figure 4B). The significance of the sample partitioning was validated using PERMANOVA (Table 1). Based on the analysis, significant effects of infection, time-points and interactions were detected (pseudo-F = 1·66–4·12, P = 0·001).

Fig. 4. (A) Canonical analysis of principal coordinates (predictive power, Q ²: 30·4%) scores plot derived from fecal TRFLP profiles of infected and control animals. (B) PCA time trajectory plots derived from TRFP data of infected animals showing differentiation of pre-infection (PI) from post infection at each time point (days 7, 14, 21, 28, 36 and 42). (C) Causal Model derived using structural equation modelling. A total of three models were built based on the sequelae of responses from cytokines, bacterial and metabolic profiles. Among the three models, model 3 was selected as the best fit model. Key: , PI; , Day 7; , Day 14; , Day 21; , Day 28; , Day 36; , Day 42.

Hypothesized biological causal model for Toxoplasma infection

s.e.m. was used to evaluate and compare alternative models of causal influences between immunological responses, metabolites variation and gut microbial compositional changes in relation to Toxoplasma infection. A total of three models were hypothesized: Model 1 assumed that the responses to infection followed a fixed chronological sequence starting with immune response, follow by metabolic response and lastly microbial response. Model 2 assumed that the chain of event starts with immune response; follow by microbial response and lastly metabolic response. Model 3 assumed that immunological responses modulate both microbial and metabolic responses while the latter two cross influencing each other.

From the three hypothesized models (Figure 4C), model 3 exhibited significantly (P < 0·05) lower Chi-square value (0) in comparison with models 1 (17·627) and 2 (12·409). The AIC and BIC scores of model 3 were recorded at 4225·80 and 4246·80, respectively. The hypothesized model 3 suggested a causative relationship between immunological parameters vs gut microbial and urine metabolic compositional changes, while a reciprocal relationship occurred between the latter two biological systems.

Correlation between host immune, gut microbial and metabolic responses to T. gondii infection

Distance based linear modelling was carried out to identify the parameters (within each systems) that exhibited strongest correlation to the corresponding responses. For instance, when we modelled the metabolic pattern using immunological variables, IFN-γ was selected as the best explanatory variable which accounted for approximately 5% of the total variation in the urine metabolic profiles of the experimental animals. However, no immunological parameter was selected when we modelled using gut microbial compositional pattern. On the other hand, changes in formate, glycine and lactate were consistent with the responses observed in gut microbial profiles and explained ~45% of the variance. Three TRFLP peaks were selected as the parsimonious model for the urine metabolic profiles. However, the taxonomic identity of the three selected TRFs was unknown as no match was obtained when comparing the size of the TRFs with the reference generated through in silico digest of RDP 16S data (Table 3).

Table 3. Distance-based linear modelling of predictors to the multivariate data cloud derived from gut microbial composition and metabolic profiles

^a Corresponding TRF showing highest correlation to the selected peak.

^b Cross-matching of virtual fragments derived from the digestion of RDP-10 16S rRNA gene database.

DISCUSSION

The present study characterized cytokine expression levels, metabolic profiles and fecal microbiota profiles of the mammalian host in response to T. gondii infection over a period of 42 days. To our best knowledge, this work presented one of the longest longitudinal cross-biological systems monitoring for T. gondii infection in a mouse model. Similar previous works were mainly targeting the responses during acute phase (cf. 7–9 days) (Heimesaat et al. 2006; Benson et al. 2009; Molloy et al. 2013) and only a few have extended their studies to understand the host systemic and microbial dynamics in the putative chronic phase (cf. > 20 days) (Arsenijevic et al. 1997; Fenoy et al. 2015). Even when available, the studies were usually focussed on host immunological responses with reduced emphasis given to host–microbiome systemic integration. Using a systems biology approach, we are able to identify the links between the immune–metabolic–microbial patterns. Additionally, we also proposed a causal model for T. gondii infection using the structural equation modelling technique.

Histopathology of haematoxylin and eosin stained paraffinized brain tissues from the current study showed that putative chronic infection was established by day 42 represented by brain cysts formation. This is in agreement with previous published works which showed that chronic infection can be established in a mouse model 2–4 weeks post-infection with T. gondii ME49 strain (Ferguson et al. 1994a , b ; Chew et al. 2012b ). However, as bradyzoite-specific genes such as BrayBAG1 and LDH2 were not detected here, the precise development stage is difficult to be ascertained. Nevertheless, the size of the brain cyst (ca. 15 µm) detected in this study is consistent with the chronic phase of other BALB/c mice Toxoplasmosis model (e.g. Chew et al. 2012a , b , 15–20 µm at 2–15 weeks post-infection; Melzer et al. 2010. 7–20 µm at day 30 post-infection) (Chew et al. 2012b ; Melzer et al. 2010).

A myriad of physiological responses had been found to be associated with different phases of T. gondii infection. For instance, comparing the energy balance and cytokine gene expression of murine model between acute and chronic T. gondii infection, Arsenijevic et al. (1997) showed that elevated TNF-α, IL-1β, IL-5 and IFN-γ coupled with increased energy expenditure were detected as early as day 1 post-infection during acute phase. In contrast, recovery/reduction of energy metabolism and cytokines expression were detected except TNF-α and IL-10. Separately, miR-132, which altered dopamine pathway of the infected mice was also found to be elevated during acute phase (Xiao et al. 2014), while a reduction was detected in the later chronic state (Li et al. 2015).

Based on our study, apparent overall immunological modulation was detected on days 7 and 14 post-infection, while recovery to the pre-infection stage was exhibited between days 21 and 42. As observed in the trajectory plot (Figure 2B), the immune response trajectory changes after day 14 indicated a change in host response, indicating the start of chronic infection. This observation is in agreement with previous work conducted in our laboratory which showed that acute infection of T. gondii in BALB/c mouse-model lasted for 2 weeks prior to the onset of chronic infection from 2 weeks onwards post infection (Chew et al. 2012b ). It also showed the reproducibility of the T. gondii ME49-infected BALB/c mouse model. Specifically, elevated serum TNF-α, IL-12p40 and IFN-γ was prominent during acute infection phase. Among the three cytokines, continued overexpression was only detected for IFN-γ in the putative chronic stage. Acute infection of T. gondii has been shown to invoke strong host Th1 immune responses (Olguin et al. 2015). In addition, consistent report by Gazzinelli et al. (1993, 1994) showed that IL-12p40 was detected as early as 2 days post T. gondii infection followed by a concomitant increase in IFN-γ . It has been suggested that TNF-α induces natural killer cells activity and IFN-γ synthesis during T. gondii infection (Hunter et al. 1994). On the other hand, Benson et al. (2009) suggested that the production of IL-12 and IFN-γ in the dendritic cells of infected animals is stimulated by the gut microflora (Benson et al. 2009). Activation of IL-12 dependent IFN-γ response during acute T. gondii infection was reported to induce the loss of paneth cells and triggers dysbiosis (Raetz et al. 2013). Further, IFN-γ was also found to facilitate the delay in T. gondii cells growth by inducing the degradation of tryptophan. Further, the cytokine also facilitates protective innate and T cell immune responses to control cerebral T. gondii infection (Sa et al. 2015).

The host immune response to infection can be observed in the urinary metabolic profile represented by changes in urinary glycoproteins specifically N-acetyl glutamate. Elevated urinary glycoproteins had been implicated in various diseases (Vivekanandan-Giri et al. 2011) and are currently used as markers of autoimmunity and malignancy. This is not surprising given that glycoproteins are involved in important cell signalling pathways and cell-to-cell interactions (Lowe, 2001). Thus, a change in host physiology due to a disease state will invariably affect these pathways, which can be observed in the urinary metabolic composition. Although only urinary N-acetyl glutamate was observed to change in the early phase infection, we postulated the possibility of a temporal delay between the expressions of host immune response translated into host urinary metabolic composition as it involved crossing several biological compartments (Yap et al. 2015). This is further supported by the results obtained on termination time point where urinary N-acetyls and O-acetyls of glycoproteins were lower in the infected group.

Reduced urinary excretion of 2-oxoglutarate, citrate and fumarate indicated perturbed host energy metabolism in response to infection. The tachyzoites of T. gondii can infiltrate and proliferate in any mammalian hosts nucleated cells. Like all other intracellular pathogens, T. gondii tachyzoites rely on host carbon sources and essential nutrients to survive (Polonais and Soldati-Favre, 2010). MacRae et al. (2012) showed that T. gondii tachyzoites utilizes glucose as carbon source to generate energy via oxidative tricarboxylic acid cycle, particularly during acute phase infection. In addition, elevated levels of urinary α-hydroxy N-valerate and β-aminoisobutyrate are indicative of perturbation in host amino acids catabolism in response to infection. Urinary excretion of these short-chain fatty acids suggested a transient shift in energy metabolism towards oxidation of both gluconeogenic and ketogenic amino acids such as isoleucine, leucine and valine (Harper et al. 1984). Changes in urinary short-chain fatty acids have been reported in various disease states, including cytomegalovirus infection (Fanos et al. 2013) and Leishmania infection (Lamour et al. 2015). Further, the increased urinary excretion of creatine with concomitant increase in urinary taurine in the infected animals also indicated a shift in host muscle energy homeostasis (Wang et al. 2008; Fanos et al. 2013).

In the current study, the evenness of fecal microbial composition of the infected animals was lower as compared to the controls during acute infection. Indeed, several authors had reported dysbiosis in gut during acute phase of oral T. gondii infection. Possible mechanisms that underlie the changes in gut bacterial community include intestinal pathology (i.e. formation of intestinal casts and loss of paneth cells) and adaptive immunity actions against commensals (Hand et al. 2012; Molloy et al. 2013; Raetz et al. 2013). However, such a difference was not observed in the samples collected on the terminal time point. This could be due to the reduced interactions between T. gondii and host gut microbiome as the protozoans begin to encyst in host tissues. Parasitic infections have been shown to affect host gut microbiome diversity (Kay et al. 2015; Lamour et al. 2015). Holzscheiter et al. (2014) showed, in a schistoshome-infected murine model, that excretion of schistosome egg was reduced when the host gut bacteria was depleted with concomitant alteration in immune and inflammatory responses as well as gut pathology. Alteration of host gut microbiome composition after T. gondii infection in this study can also be indirectly represented through perturbation in the metabolism of mammalian–microbial co-metabolites namely urinary hippurate, phenyacetylglycine, trimethylamine N-oxide, dimethylglycine and acetate. Toxoplasma gondii-induced changes in mouse urine included elevated dimethylglycine and trimethylamine N-oxide, which suggested alteration in host gut microbial activity since dimethylglycine and trimethylamine N-oxide are intermediates of dietary carnitine and choline, which are degraded to trimethylamine by gut microbe before being transported to the liver and oxidized to trimethylamine N-oxide (Smith et al. 1994). In addition, T. gondii infection also induced elevated urinary phenylacetylglycine and decreased urinary hippurate. Both phenylacetylglycine and hippurate are known mammalian–microbial co-metabolites (Yap et al. 2015). Nicholls et al. (2001) showed that both hippurate and phenylacetylglycine were absent in the urine of germ-free rats and were only observed several days later after exposure to laboratory environment implicating their link to microbial metabolism. Hippurate is formed in the liver through glycine conjugation of benzoate produced from bacterial metabolism of plant phenols (Schwab et al. 2001). Further, the gut microbiota is involved in breakdown of aromatic amino acids such as tryptophan and tyrosine (Diaz et al. 2001). The significant changes in urinary mammalian–microbial co-metabolites observed in the current study indicated that T. gondii infection can induce change in bacterial metabolism of phenolics in mammals.

It is noteworthy that several metabolites had shown a consistent trend in both acute and the putative chronic phase. For instance, upregulation of α-hydroxy N-valerate and guanidinoacetate were detected in the urinary metabolite profiles on days 14 and 42. Guanidinoacetate has been found to increased acethylcholinesterase activity and brain dysfunction (Zugno et al. 2008). This may indicate the involvement of the metabolite in symptomatic behavioural changes of the disease. Conversely, a consistent reduction in gut microbial–mammalian co-metabolites including phenylacetylglcine and trimethylamine N-oxide was detected (Zheng et al. 2011). The strong host metabolic–microbial interaction corroborated with the CAP ordination of TRFLP and NMR profiling where consistent microbial and metabolic clustering at day 14 were recorded (Figures 2C and 4A).

When compared across Figs. 2, 3 and 4, it is notable that some data points were inconsistent with each other (e.g. cytokines profiles vs TRFLP profiles). The discrepancy was due to overlapping of some of the data points displayed in the respective scores plots with each other under the dimension projected. In addition, during multivariate statistical analyses, some of the outlying samples were removed due to spurious instrumental reading in order to avoid compromising statistical model predictivity and robustness. The actual sample sizes used for different matrices (i.e. metabolic profiles, cytokines profile and TRFLP profiles) are given in Table 2, and Tables S1 and S2 in the Supplementary Materials (available from http://journals.cambridge.org/PAO).

One major objective in the study of disease progression is to establish sequelae of infection from acute to the putative chronic stages. By comparing the immune-metabolic-microbial patterns, we were able to show that conspicuous immunological alteration occurred earlier (day 7) than host metabolic and microbial modulation (day 14) during acute infection phase. We further established the directionality of the putative causal relationship using structural modelling equation. Based on our model, host metabolism and gut microbial compositional changes were likely to be a function of immunological status (Model 3, Fig. 4C). It is interesting to note that the role of gut microflora in immune regulation towards infection has been well studied (Benson et al. 2009; Caballero and Pamer, 2015). However, the reverse had received little attention. We acknowledged that the model will be more complete if parasites related factor such as cyst size, gene expression pattern is included in the model. However, since that the parasitic factors are equal across the three proposed models (Fig. 4C), Model 3 is likely to remain as the best model even if additional factors are included. Using multivariate regression, we also exhibited close interaction between gut microbial and metabolic dynamics. For instance, the 45% of the variance in TRFLP profiles can be explained by NMR profiles, while 12·2% was explained vice versa (Table 3).

Concluding remarks

The systems biology approach used here to understand the effect of T. gondii infection in mammalian host sheds new light into host metabolic as well as host gut microbe responses to infection. Using structural equation modelling, we were able to show statistically the host immune system as the driver for changes observed in host metabolism and gut microbe composition. Clear differences can be observed between acute and the putative chronic phases of T. gondii infection represented by host immune responses and host energy metabolism. Furthermore, gut microbe compositional change was observed only during acute phase of infection but not during the chronic phase of infection although significant changes in host metabolism can still be observed during the chronic phase of infection. The metabolites identified during the chronic phase infection warrant further investigations on their potential as diagnostic markers for chronic T. gondii infection.