Toxoplasma infection in intermediate host species closely associates with inflammation. This association has led to suggestions that the behavioural changes associated with infection may be indirectly driven by the resulting sustained inflammation rather than a direct behavioural manipulation by the parasite. If this is correct, sustained inflammation in chronically infected rodents should present as widespread changes in the gastrointestinal microbiota due to the dependency between the composition of these microbiota and sustained inflammation. We conducted a randomized controlled experiment in rats that were assigned to a Toxoplasma-treatment, placebo-treatment or negative control group. We sacrificed rats during the chronic phase of infection, collected their cecal stool samples and sequenced the V3-V4 region of the 16S rRNA gene to characterise the bacterial community in these samples. Toxoplasma infection did not induce widespread changes in the bacterial community composition of the gastrointestinal tract of rats. Rather, we found sex differences in the bacterial community composition of rats. We conclude that it is unlikely that sustained inflammation is the mechanism driving the highly specific behavioural changes observed in Toxoplasma-positive rats.

Infection of laboratory rats

We obtained 50 eight-week old male and female Wistar Hannover rats (HanTac:WH; 24 male rats, 24 female rats and 2 control rats) from InVivos Private Limited, Singapore. All rats were housed in pairs with ad libitum supply of food and water (12:12 h light-dark cycle, lights on at 0700 h). Rats were allowed to acclimatize to the vivarium for one week prior to the start of experiments. We randomly assigned rodent cages to two experimental groups based on randomly generated numbers[PT1] . Rats assigned to the treatment group received five million Toxoplasma gondii tachyzoites in 500 µl buffered saline via the intraperitoneal route[PT2] . Those assigned to the control group received only saline via the same route. Toxoplasma parasites were maintained in monolayers of immortalized human foreskin fibroblasts (ATCC # CRL-4001) and harvested through syringe lysis with a 25-gauge needle, as previously described (Tan, Tong, & Vyas, 2020; Vyas, Kim, Giacomini, Boothroyd, & Sapolsky, 2007). A type II Prugniaud (Pru) Toxoplasma gondii strain was used in the experiments. We used two of the 50 rodents as negative controls (i.e. rats were not inoculated with anything but instead used to generate sampling blank controls). Both negative control rats were subject to equivalent experimental conditions for the duration of our laboratory experiment and sample processing. All animal experimental procedures were conducted at the Nanyang Technological University and reviewed and approved by the institutional Animal Care and Use Committee (# A18074).

DNA extraction, PCR and sequencing

All rodents were euthanised at eight weeks post-infection to ensure they had reached the chronic phase of the Toxoplasma gondii infection (J. P. Dubey, Speer, Shen, Kwok, & Blixt, 1997). Rodents were then presented on sterile boards and their fur sprayed with ethanol to minimise contamination of internal tissues. Rodent ceca were then harvested under a flow hood, snap-frozen in liquid nitrogen, and stored at -80^oC. Total DNA was extracted from 250 mg of stool collected from within each cecum sample using the QIAamp PowerFecal Pro DNA Kit (Qiagen #51804), following instructions from the manufacture. DNA samples were then stored at -80^oC prior to shipping to Adelaide, Australia, on dry ice. The experimenter was blind to the infection status of all rodents during the isolation of DNA, and the samples were processed in a random order.

DNA amplification was performed at the Australian Genome Research facility (AGRF). AGRF amplified the bacterial 16S rRNA V3-V4 gene region using forward primer 341F (5’- CCTAYGGGRBGCASCAG- 3’) and reverse primer 806R (GGACTACNNGGGTATCTAAT). A modified Illumina two-stage PCR library preparation technique was then used, adapted from (Illumina, 2014). Primary PCR amplicons were generated using region-specific primers coupled with nextera stubby forward (5’-TCGTCGGCAGCGTCAGATGTGTATAAGAGACAG- 3’) and reverse (5’-GTCTCGTGGGCTCGGAGATGTGTATAAGAGACAG- 3’) adapter sequences and Platinum SuperFi (ThermoFisher, Australia). PCR reactions were cycled at 95 °C for 7 min; 29 cycles of 94 °C for 30 s, 50 °C for 60 s and 72 °C for 60 s; and final extension at 72 °C for 7 min. The first stage PCR was cleaned using a 1.8x ratio of Sera-Mag Magnetic Speedbeads (Merck, Australia), and samples were visualised on 2% SYBR E-Gel (Thermo-Fisher). A secondary PCR was performed with PrimeSTAR Max DNA Polymerase (Takara, Japan) to add sample indices and Illumina P5/P7 sequencing adaptors to the first stage amplicons using unique-dual indexing. Second stage PCR reactions were cycled at 98 °C for 10 s; 8 cycles of 98 °C for 10 s, 50 °C for 5 s and 72 °C for 10 s; and holding at 10 °C. The resulting amplicons were cleaned again using a 1.8x ratio of magnetic beads, quantified by fluorometry (Promega QuantiFluor) and then normalised. The eqimolar pool was cleaned a final time using a 1.5x ratio of magnetic beads to concentrate the pool and then measured using a High-Sensitivity D1000 TapeStation assay (Agilent, USA). The pool was diluted to 5nM and molarity was confirmed again using a HighSensitivity D1000 TapeStation assay (Agilent, USA).

Sequencing was performed at the AGRF using a MiSeq with V3 600 cycle chemistry (300 base-pair paired-end), 25% PhiX and clustered at 7 pM.

Sequence analysis

Amplicon sequence variants (ASVs) were generated using the QIIME2 (version r2021.11) analysis pipeline (https://doi.org/10.1038/s41587-019-0209-9). Paired end 300bp Illumina reads were first trimmed on primers using the cutadapt trim-paired QIIME2 plugin with the default parameters. Sequences were denoised, overlapped, dereplicated, and filtered for chimeras using the DADA2 QIIME2 plugin, with the following parameters: --p-trunc-len-f 240, --p-trunc-len-r 220, --p-max-ee-f 2, --p-max-ee-r 4, --p-n-reads-learn 1000000, and --p-chimera-method consensus (Martin, 2011). For classifier training the 341F and 806R amplicon region was extracted from the QIIME compatible SILVA SSU database release 138  using the qiime feature-classifier extract-reads command (Quast et al., 2012). The classifier model was then fit and trained using the commands qiime feature-classifier fit-classifier-native-bayes and qiime feature-classifier classify-sklearn with the default parameters. Run scripts and conda env files used in this analysis can be found at the projects git repository https://bitbucket.org/dpi_data_analytics/202202taggart/src/master/.   

Data cleaning

We discarded ASVs that were not identified as bacteria, were unidentified at the phylum level, assigned as chloroplast or mitochondria, or did not occur in at least two samples. We identified and removed contaminating sequences using the decontam package in R via the isContaminant() function as recommended in the decontam user guide for high biomass sample types (Davis, Proctor, Holmes, Relman, & Callahan, 2018; R-Core-Team, 2020). We used the prevalence method which identifies contaminants by increased presence in negative controls. Otherwise, default settings were used throughout. Subsequently, all taxa identified as contaminants and the negative control samples were removed from further microbiota data analysis. After these data cleaning steps, a total of 829 ASVs from 48 samples (2,474,934 total sequences) remained.

Statistical analysis

We used R software and the microbiome data analysis framework of the phyloseq package (McMurdie & Holmes, 2013; R-Core-Team, 2020). Based on rarefaction curve plots two samples with very low sequence counts (5AB/11BF) were discarded and ASV abundance data were rarefied to the minimum sample sequence read depth of the remaining samples (17,577 sequences) for alpha and beta diversity analyses described below (Figure S1). We estimated alpha diversity in each sample via the exponential transform of Shannon Index values to derive the effective number of ASVs (Jost, 2006). These data were visualised using boxplots, a heat map and barplot grouped by treatment and gender. We tested for differences in alpha diversity between the groups using the Kruskal-Wallis rank sum test in base R (version 4.1.2), followed by the Dunn test for multiple comparisons to determine which groups differed from other groups, with Bonferroni-adjusted threshold P values. The multiple comparison testing used default two-sided P values and an alpha = 0.05 nominal level of significance.

 We visualised differences between sample microbiota (beta diversity) using non-metric multidimensional scaling (NMDS) ordination of Bray-Curtis distances. Two outliers samples, 2BF (female Placebo) and 7BF (female Toxoplasma), were excluded from the beta diversity analysis due to outlying high NMDS1 and NMDS2 scores - i.e. microbiota data with the outlier removed (total 829 ASVs; 44 samples; 2,377,578  sequences) was then rarefied again to the lowest sample depth (18165 sequences) for this beta diversity analysis. We used permutational multivariate analysis of variance (PERMANOVA), as implemented in the adonis() function within the vegan package, to test for compositional differences by treatment and sex between microbiota sample groups (Oksanen et al., 2020). This was followed by testing for homogeneity of group dispersions via the betadisper() function (also from vegan).

 We also assessed differentially abundant ASVs (again using microbiota data with the two outlier samples 2BF and 7BF removed) between the Toxoplasma vs placebo treatments using the ANCOM-BC package (version 1.4; 10.1038/s41467-020-17041-7). Differential ASV abundance was calculated for both treatment and sex using the ancombc function with the default options, except formula=Treatment+Sex, group=Treatment, lib_cut=1000, neg_lb=TRUE, and conserve=TRUE. P-values were adjusted for multiple comparisons using the default Holm-Bonferroni method and <=0.05 was used as a significance cut-off (Holm, 1979).