Buffalo Gut Microbes May Affect the Host Th2 Responses During Fasciola Gigantica Infection via SCFAs Metabolism and TLR Signaling Pathway

Background: Helminth-induced Th2 responses are essential to modify the structure and diversity of gut microbes. However, observations have come mainly from studies of helminth-infected humans or rodent models. Very little research has been conducted in veterinary animals. Methods: In this study, we searched for links between microbiota and Th2-biased responses during the time course of Fasciala gigantica infection in buffaloes.16S rRNA gene amplicon and metagenome sequencing were applied to analyze the structure and function of the gut microbiota. Results: Both alpha and beta diversities decreased during infection, and gut microbes differed considerably across different sections of the gut at different stages. Immune responses changed when the microbiota traverses the gut wall into the peritoneal cavity, in line with the changes in Th2 response induced by F. gigantica infection. We found that the order Coriobacteriales was greatly decreased at the early stages in which the Peptostreptococcaceae and Family_XIII families are closely linked to the upregulation of IgG1 and IL4, respectively. The F. gigantica infection signicantly reduced short-chain fatty acid (SCFAs)-producing microbes, reduced the concentrations of gut SCFAs and downregulated the SCFAs-producing metabolic pathways. In addition, The microbes associated with TLR2 increased and showed similar trend to the TLR2 and Th2 cytokine production during infection, suggesting that bacteria ligands might recognize TLR2 and subsequently induce a Th2-biased response. Conclusions: Our data show that buffalo gut microbes may affect the host Th2 response during F. gigantica infection via the SCFAs metabolism and TLR signaling pathway. These ndings provide new insights into the relationship between F. gigantica–microbiota-host, which may provide new potential therapeutic targets for prevention and control Fasciolosis.


Background
Ruminants comprise a signi cant proportion of the livestock species in China. Ruminants and gut microbiota, counting at least 10 14 species of microbes, have coevolved over millions of years, and the microbes play important roles in supporting the function of the immune system [1,2]. It is reasonable that both microbes and helminths interact with each other, and such interaction may affect host homeostasis.
For example, reduced Enterobacteriaceae has been attributed to the upregulation of Th2 cytokines [3]. Also, attenuation of the gut microbiota composition and IL17 expression were observed when STAT6 and IL13 de cient mice were infected with Nippostrongylus brasiliensis [4].
As the most common de nitive host of Fasciola gigantica, buffaloes are also the most economically important ruminant species [5]. Infection of buffaloes with F. gigantica is common in southern China [6]. The outbreaks of F. gigantica-caused human disease in Yunnan, China, suggested that F. gigantica infection is a serious public health problem [7]. However, F. gigantica ukes mainly parasitize the liver of buffaloes, and the ukes are e cient immune-modulators in that they release immune effectors that can exploit the host immune responses. This makes it highly di cult to control such helminth diseases. The most recent data from buffaloes rather than rodent models showed that during early F. gigantica colonization in buffaloes, a modest increase of Th2 cytokines appeared, but the immune system was suppressed during chronic F. gigantica infection [8]. The results showed a mixed Th1/Th2 cytokine pro le in the middle infection stage and a biased Th1/Treg response in the late infection stage, suggesting a balance of in ammatory and immune-regulatory mechanisms in the liver of buffaloes during F. gigantica infection [9].
In the animal gut, TLRs recognize microbe-associated molecular patterns (MAMPs) that are shared by many types of microbes [10]. Previous studies demonstrated that activation of TLR2 induces a Th2biased response via production of IL10 and IL13 [11][12][13][14][15][16][17][18]. Previous study also showed that the expression of Toll-like receptor (TLR) and NOD-like receptor (NLR) genes in F. gigantica-infected buffaloes altered to evade host immune defenses [8]. These ndings suggest that TLRs function as multiplayers in the regulation of adaptive immune responses. However, few studies on the relationship between microbes and Fasciola-induced immunity has been documented for now. Thus, we asked that whether microbes in different gut may affect F. gigantica-associated immunity.
In this study, we aimed to obtain a much clearer picture of the functional links between host and microbiota during the time course of F. gigantica infection in buffaloes. Both 16S rRNA gene amplicon sequencing and metagenome sequencing were applied to investigate changes in abundance and functional associations of microbes. Correlation networks were generated applying amplicon sequence and immune-related gene expression data in order to identify links between F. gigantica-associated host immunity and gut microbiota.

Metacercariae, animal infection and sample collection
Thirty-ve swamp buffaloes between the ages of 8 to 10 months were purchased from a local buffalo farm in the Guangxi Zhuang Autonomous Region, PR China. Infective metacercariae were prepared as in our previous protocol, and animal infection experiments using liver uke-free buffaloes for F. gigantica infection were conducted as described in our previous published paper [9]. Brie y, buffaloes were kept in separate protective fence and provided with commercial feed and clean drinkable water during the research period. All buffaloes were con rmed as liver uke negative by fecal examination and F. gigantica-speci c IgG-antibody-based ELISA. To eliminate any potential uke infection, all buffaloes were treated with a single dose of triclabendazole (5% w/v) 30 days before the experimental F. gigantica infection. All buffaloes were classi ed into 7 groups (5 buffaloes/group). Buffaloes in Group I was mocked with PBS. Each buffalo in Groups II-VII was incoulated with 500 metacercariae via oral gavage.
Buffaloes were sacri ced and all faeces in rumen (Ru), duodenum (Du), jejunum (Je), colon (Co) and rectum (Re) and tissue samples were collected at 3, 10, 28, 42, 70, and 98 days post infection (dpi). The liver in each animal was used to examine pathological lesions and the uke burdens. Fluke eggs were examined by ltering bile uid via a 0.15 mm size mesh. Only the samples with pathological lesions on liver were collected. Liver tissues were kept in 10% PBS-buffered formalin, and para ns were mounted onto glass slides, and stained with hematoxylin and eosin (H&E). Stained sections were examined under microscope at 400 × magni cation and imaged (Zeiss Axio Imager). The liver samples for RNA extraction were kept in RNA store buffer (Tiangen Biotech, Beijing, China), stored at -80 °C until experimental use. For the faeces samples, A total of 120 samples from three stages (stage PCL refers to 3 dpi, when worms migrate from the gut and peritoneal cavity to the liver surface; stage LP refers to 10−28 dpi, when worms live in liver parenchyma; stage BD refers to 42−98 dpi, when worms live and lay eggs in the bile duct) and ve gut sections (Ru, Du, Je, Co and Re) were prepared. In each stage and each gut section, three samples were from control buffaloes, and ve samples were from F. gigantica-infected animals (details see S1 Table and Fig. 1). We divided the entire infection into three stages so we could attempt to understand whether the rst contact when newly excysted juveniles cross the gut wall into the peritoneal cavity (stage PCL) would affect the gut microbiota and immune responses, and whether this would differ during parasite migration into the liver capsule and parenchyma (stage LP) or in the bile ducts where they develop into sexually mature adults (stage BD).

Liver histopathological evaluation and immunohistochemical study
Samples of liver tissue showing pathological lesions were collected from each animal, and the hematoxylin and eosin (H&E) staining was carried out according to the previous method [9]. DNA damage has been reported to be associated with parasite infection, including O. viverrini [19], trypanosome [20] and Leishmania [21]. To investigate the DNA damage to liver tissues post F. gigantica infection, 8nitroguanine (2 µg/ml) in the liver tissues was assessed by an immuno uorescence labeling method as described previously and was done by Professor Ma personally [22]. To detect the co-expression of NF-B and 8-nitroguanine, immunohistochemical staining was performed using anti-NF-B mouse monoclonal antibody (1:300, DAKO, Glostrup, Denmark) and anti-8-nitroguanine rabbit polyclonal antibody as described previously [22]. The scores were evaluated for each specimen based on the degree of staining: -, negative; +, less than 25%; ++, 25-50%; +++, 50-75% and ++++, more than 75% of the cells in liver tissues. Sections were blindly analyzed by two anatomists who were unaware of the status of the liver tissues.  expression of TLRs  (TLR1, TLR2, TLR3, TLR4, TLR5, TLR6, TLR7, TLR8 and TLR9) by qPCR, which was performed with the SYBR® Green Master Mix (Life Technologies, Carlsbad, CA, USA) in a 96-well optical plate on an ABI 7500 PCR instrument following the manufacturer's protocol. All primers are listed in S13 Table. DNA extraction and 16S rRNA sequencing A 200 mg faeces sample was taken from each buffalo's gut section (Ru, Du, Je, Co and Re) contents. Bacterial genomic DNA was extracted at Novogene Bioinformatics Technology Co., Ltd by using a Tiangen kit in accordance with the manufacturer's protocol. The quality and integrity of each DNA sample was determined by electrophoresis in 1% agarose gels with Tris-acetate-EDTA buffer, and the DNA concentration was quanti ed with a NanoDrop ND-2000 spectrophotometer. Sequencing libraries were generated by using a TruSeq® DNA PCR-Free Sample Preparation Kit (Illumina, USA) according to the manufacturer's protocol. The library quality was evaluated on the Qubit@ 2.0 Fluorometer (Thermo Scienti c) and Agilent Bioanalyzer 2100 system. The library was sequenced on an Illumina HiSeq2500 platform. The 16S rRNA community pro les were characterized by primer-speci c PCR ampli cations that were performed by targeting the V3-V4 region of the 16S rDNA gene (341F: CCTAYGGGRBGCASCAG, 806R: GGACTACNNGGGTATCTAAT) [23]. All PCR reactions were conducted with Phusion® High-Fidelity PCR Master Mix (New England Biolabs).

Taxonomic analyses of sequenced reads
Paired-end reads were merged using FLASH software. Raw data were ltered using QIIME quality lters [24,25]. The UPARSE pipeline was used to pick operational taxonomic units (OTUs) at an identity threshold of 97%, and chimeras were removed using the UCHIME algorithm [26,27]. The RDP classi er tool was used to assign taxonomic data to each representative sequence of OTUs. Alpha diversity (community richness) was evaluated based on the rarefaction curves generated on metrics, and the number of OTUs in each sample was determined (Chao1 metric and Shannon index). The beta diversity and taxon composition were analyzed by QIIME (version 1.7.0) for calculating both weighted and unweighted UniFrac. The unweighted UniFrac phylogenetic distance metric was analyzed by using a principal coordinate analysis (PCoA) and unweighted pair group method with arithmetic mean (UPGMA) clustering. Cluster analysis was also preceded by principal component analysis (PCA) using the FactoMineR package and ggplot2 package in R software (Version 2.15.3). The PICRUSt (Phylogenetic Investigation of Communities by Reconstruction of Unobserved States) [28] was used to generate the KEGG (Kyoto encyclopedia of genes and genomes) pathway, and all functions were categorized at levels 2 and 3.

Metagenomic sequencing, gene catalog construction and annotation
Twelve samples from different stages and lumen sections were chosen based on the 16S rRNA clustering results. All samples were paired-end sequenced on the Illumina platform at the Novogene Bioinformatics Technology Co., Ltd. The reads aligned to the buffalo genome (alignment with SOAP Aligner, parameters [29]: identity ≥ 90%, -l 30, -v 7, -M 4, -m 200,-x 400) were removed. The assembly of reads was executed using SOAP de novo (Version 2.04, parameters: -d 1 -M 3, -R, -u, -F, -K 55) [30,31]. We mapped the clean data against scaffolds using SOAPAligner (Version 2.21, parameters: identity ≥ 90%, -m 200 -x 400). Genes (longer than 100 bp) were predicted on scaftigs (longer than 500 bp) using MetaGeneMark (version 2.10). A non-redundant gene catalog was constructed with CD-HIT (version 4.5.8, parameters: -c 0.95, -G 0, -aS 0.9, -g 1, -d 0) [32,33] using a sequence identity cut-off of 0.95. Reads were realigned to the gene catalog with SOAPAligner (version 2.21) using the following parameters: identity ≥ 95%, -m 200 -x 400. Only genes with ≥ 2 mapped reads were present in a sample [34]. The abundance of genes was calculated by counting the number of reads and normalizing by gene length. Genes were annotated using DIAMOND (V0.7.9, parameters: blastp, -e 1e-5) [35]. The Unigenes and all bacteria, fungi, archaea and viruses in the NR database were aligned. LCA was used to determine the species annotations [36]. Based on the LCA results, gene abundance in each classi cation hierarchy (kingdom, phylum, class, order, family, genus, species) was determined, and relative abundance was obtained by Krnoa analysis [37]. All genes in our catalog were aligned to the KEGG database using DIAMOND (V0.7.98, parameters: blastp, -e 1e-5). Each protein was assigned to the three databases by the highest scoring annotated hits (at least over 60 bits) [38]. The abundance of KEGG orthology was determined by summing the abundance of genes annotated to the same feature.

Immunological parameters
Immunological data of serum IgG and spleen cytokines were used to perform the statistical analysis of linear models. Redundancy analysis (RDA) and canonical correspondence analysis (CCA) were taken from our original experiments. Most of these data were included in a paper titled "Th2-related cytokines responsible for Fasciola gigantica infection and evasion in the natural host swamp buffalo" by our research team [39] that had been accepted by Veterinary Parasitology journal. Brie y, qPCR was used for detection of gene expression.
qPCR program was as follows: 95 °C for 5 min, then 40 cycles at 95 °C for 5 s, and 60 °C for 30 s. The qPCR data was exported using Bio-Rad CFX Manager. The 2-ΔΔCq method was used to calculate the gene expression the level, and the bubalus glyceraldehyde-3-phosphate dehydrogenase (GAPDH) gene was used to normalize gene expression data. For antibody detection, the anti-ESP isotype antibodies IgG, IgG1, and IgG2 were detected using a 96-well ELISA plate (JET BIOFIL, Guangzhou, China) coated with 100 µL F. gigantica ESP (2.5 µg/mL). Two hundred microliter per well of 1% gelatin in PBS was used to block plates and 100 µL calf serum diluted (1:200) in PBST was added. Then, 100 µL horseradish peroxidase (HRP)-conjugated sheep anti-bovine IgG, IgG1, and IgG2 (AbD Serotec, UK) polyclonal antibodies (1:40,000 in PBST) were used to detect the antibody levels. Etramethylbenzidine TMB (100 µL/well; Solarbio, Beijing, China) was used to induce development, and 2M H2SO4 (50 µL/well) was added to stop the reaction. The absorbance at 450 nm was recorded by iMARK Microplate Absorbance Reader (Bio-Rad Laboratories, Hercules, CA, USA). All of the qRT-PCR and ELISA data were analyzed by GraphPad Prism 5 (La Jolla, CA, USA).
We re-analyzed these data with our 16S rRNA sequencing data to determine the relationships between immunological outcomes and diversity indices.

Determination of SCFAs in feces
Feces were collected as stated above from individual buffaloes. Fecal samples (50 mg) were added to 2 ml water (1:3 phosphate water solution, pH 2-3), shocked and resuspended for 2 min. One milliliter of diethyl ether was added and incubated for 10 min before the solution was centrifuged at 2,000 g for 10 min. The supernatants were mixed with internal standard cyclohexanone solution in ether (0.1 ml of 1000 mg/L) and ltered by a 0.45 µm microporous membrane. Samples were analyzed by GC-MS (Thermo Scienti c ISQTM LT) within 12 h. One microliter of sample was injected into GC-MS that was equipped with a TG-Wax column (30 m × 0.25 mm × 0.25 µm). Helium was used as the carrier gas at a ow rate of 1.0 mL/min. The injection temperature was 240 °C, and the GC temperature program was as follows: begin at 100 °C, increase to 150 °C at 5 °C/min, then increase to 240 °C at 30 °C/min and hold for 30 min. The split ratio was 75:1, and the ion source temperature was 200 °C. Concentrations of SCFAs were analyzed by using single ion monitor (SIM) scan mode (quantitative ion 60, 73) and calculated by using the internal standard method.

Statistical analysis
Statistical analysis of qPCR measurements was performed in GraphPad Prism 5 (Graph Pad Sofware, La Jolla, CA, USA). One-way analysis of variance (ANOVA) and Tukey's multiple-comparison test were conducted to assess differences between groups. Two-way ANOVA and Bonferroni post-tests were performed to analyze differences between the groups during the time course of infection. Adonis, Anosim, MRPP and Amova were used to test for signi cant shifts in microbial community composition. T-test, MetaStat and LEfSe were used to test for signi cant differences among groups. Pearson rank correlations were calculated to determine the relationships between immunological outcomes and diversity indices. The species data matrix was rst analyzed with the detrended correspondence analysis (DCA) to determine which model would be most suitable for correlation analysis [40]. Based on the results of the DCA, RDA or CCA was used for bacterial community analysis [41] to analyze the general relationships between each microbial community and immunological outcomes. Statistical analyses were performed using GraphPad Prism 5. Results were considered to be statistically signi cant with p < 0.05 (*), < 0.01 (**) and < 0.001 (***).

Competing interests Ethics approval and consent to participate
All experiments were evaluated and approved by the Animal Ethics Committee of Guangxi University, China. Animals used in this study were treated in line with good animal practices as demanded by the Animal Ethics Procedures and Guidelines of the People's Republic of China. The protocol assigned by the IACUC/ethics committee approved our animal experiments.

Con rmation of F. gigantica infection in buffaloes
In this study, 35 buffaloes were classi ed into 7 groups (5 buffaloes/group). Buffaloes in Group I was mocked with PBS. Control buffaloes were euthanized at the start of the experiment to obtain baseline values for subsequent investigations. Buffaloed in Groups II-VII were infected with 500 metacercariae via oral gavage. F. gigantica infection was processed as diagrammed in Fig. 1A and was con rmed in all challenged buffaloes by observing pathological lesions in livers and F. gigantica eggs in feces (Fig. 1B).
In ltration of in ammatory cells, red blood cells and fatty degeneration was observed during infections, and no such pathological changes were found in controls. Fluke eggs were detected in the feces of infected animals at 98 dpi, and no F. gigantica eggs were observed in controls. Since DNA damage has been reported to be associated with uke infection, such as Opisthorchis viverrini [19], we wandered whether this would happen during F. gigantica infection. Our results showed that DNA damage occurred and progressed during the F. gigantica infection (Fig. 1C). Speci cally, the 8-nitroguanine staining on 3 dpi, 10 dpi and 42 dpi was mainly observed in the nucleus of the same in ammatory cells and signals were stronger than in controls, and the co-expression of NF-B and 8-nitroguanine also increased along with the infection time course. Livers of the uninfected buffaloes appeared negative for 8-nitroguanine immunoreactivity. These results con rmed the successful infection of F. gigantica during our experimental time course.
TLRs expression during F. gigantica infection In this study, we divided the entire infection process into three stages, which provided us chances to address whether the rst contact when newly excysted juveniles cross the gut wall into the peritoneal cavity (stage PCL) would affect the gut microbiota and immune responses, and whether this would differ during parasite migration into the liver capsule and parenchyma (stage LP) or in the bile ducts where they develop into sexually mature adults (stage BD). We rstly investigated the TLRs expression in mesenteric lymph nodes (MN), hepatic lymph nodes (HLN), liver, and spleen samples during infection. The results revealed that TLR1, 2, 3, 4, 6, 7 and 9 in MN, spleen and HLN increased signi cantly only at stage PCL compared with the control group. TLR1, 3, 4, 6, 7 and 9 in these tissues were sharply decreased at stages LP and BD compared with stage PCL ( Fig. 2A These results suggest that stronger TLRs responses were induced by F. gigantica at stage PCL rather than at stages LP and BD in MN, spleen, HLN and liver.

Characteristics of gut microbiota post F. gigantica infection
In this study, a total of 120 samples from three stages (stage PCL, LP, and BD) and ve gut sections (Ru, Du, Je, Co and Re) were collected. In each stage and gut section, 3 samples were from control buffaloes, and 5 samples were from F. gigantica-infected buffaloes (S1 Table and Fig. 1). 16S rRNA sequencing was applied to reveal the microbiome differences in the various gut sections at different infection stages after F. gigantica infection. The 16S rRNA sequencing reads revealed for each microbiota on average 2,751 operational taxonomic units (OTUs) from an average of 10,063 reads (S2A Fig; S2 and S3 Tables). The OTUs numbers revealed that microbiota in F. gigantica-infected buffaloes was signi cantly decreased in comparison with the non-infection group. Speci cally, Clostridiales was the most abundant microbe in all groups, and this taxon was slightly increased in infection groups compared with controls. The second abundant microbe was Coriobacteriales, which decreased at stages PCL and BD but increased at stage LP in comparison to controls (S2A, B Fig).
Alignments and phylogenetic assignments of 16S tags were performed at 97% similarity, and the analysis identi ed 45 phyla, 337 families and 716 genera in the bacterial domain (S4-6Tables). Among the phyla, the Firmicutes, Actinobacteria and Bacteroidetes were the major bacterial phyla, followed by the Proteobacteria, Verrucomicrobia, Euryarchaeota, Cyanobacteria and others (Fig. 3A). The microbial distribution pattern at the phylum level was not signi cantly different in the same gut section at different infection stages, and similar trends were observed in different gut sections at the same infection stage. However, a signi cant difference was observed at the genus level (Fig. 3B, S7 Table). For example, Ruminococcaceae_UCG-005 in the Co section at PCL and LP stages was signi cantly higher postinfection compared with the non-infection group (S7 Table).
Taxa clustering demonstrated that 37 families from 10 phyla were changed between control and infected groups, between different gut sections, and between different infection stages (Fig. 3C, S8 Table). The alpha diversity (richness and evenness) was determined, and similar trends were observed as assessed  Table). Microbiota changes in different gut sections at different infection stages were also observed. For example, microbiota in PCL. Ru. Inf was signi cantly decreased compared with the BD. Je. Inf group (S9 Table). The ower gure also displays the changes listed above (S2B Fig). PCoA based on the weighted UniFrac metric revealed a good separation of trends among all groups without any speci c classi cation (Fig. 3D).

Impact of F. gigantica infection on gut microbial composition at different infection stages
In this study, buffaloes were infected with F. gigantica metacercariae, and samples were taken from PCL, LP and BD stages and ve gut sections at each stage. PCoA was applied to analyze the relative abundance of microbiota in groups at different infection stages from 16S rRNA transcripts. The results revealed that microbiota was highly shifted at stage BD (Fig. 4A). Generally, at the family level, the greatest changes occurred in Lachnospiraceae, Methanosaetaceae, Pseudomonadaceae, and Vibrionaceae families during infection. As for Lachnospiraceae, no such bacteria were observed in controls or at PCL and BD infection stages. Methanosaetaceae was only observed in the BD infection stage. Pseudomonadaceae was decreased dramatically at all infection stages, while Vibrionaceae was increased signi cantly at stage PCL and then decreased sharply at stages LP and BD (Fig. 4B).
We then investigated whether the F. gigantica infection affected microbiota composition between infection stages. Beta diversity index (S4A Fig), PCoA (S4B Fig) and NMDS (Fig. 4B) showed a signi cant difference in the microbiota composition among all groups. Furthermore, the beta diversity box demonstrated signi cant differences between control and stage PCL, control and stage LP, as well as between PCL and LP, LP and BD, and PCL and BD infection stages (Fig. 4C). We then analyzed the community structure differences between groups during infection using the Anosim method, most changes appeared to happen at stages PCL and BD. At stage PCL, signi cant differences were found between Ru and Co, Du, Je and Re. At stage BD, Je intestine showed microbiota changes compared with Co, Re and Ru sections. Signi cant differences were also observed between Du and Re and between Du and Ru at stage LP. This suggested that microbiota in the same intestine section at different infection stages exhibited different pro les when F. gigantica infection occurred (Fig. 4D).
To investigate compositional dissimilarity among different stages, T tests were used based on the fact that the relative abundances of the all analyzed samples demonstrated signi cant differences between the control and the corresponding infection stage. By comparing with the control group, eight genera at stage PCL, six genera at stage LP and one genus at stage BD were found to be signi cantly changed (S5 Fig). Among these, Acetitomaculum and Akkermansia were signi cantly changed from stage PCL to LP, and this trend continued from stage LP to BD, while Ruminiclostridlum_6 and Desulfolibrio were signi cantly altered between stages PCL and BD, and the trend continued from stage LP to BD (S10 Table). Furthermore, at both family and genus levels, abundance shu ing and replacement of community members caused the community shifts. At the genus level, microbiota in Co, Je, Re, and Ru sections at stage PCL was strongly changed in comparison to the control group. A similar phenomenon was observed in Du and Ru at stage LP. For one speci c intestine section, when comparing different infection stages, we found that the relative abundances and community structure were signi cantly altered. Regarding the Co, Je, Re and Ru intestine sections, the microbiota changed between PCL and LP, and PCL and BD stages, while in the Du and Ru intestine sections, the microbiota was only altered at the LP and BD stages.
We investigated whether the bacteria changed along with infection stages in the above intestine sections. For example, in the Co intestine section, the genera Alistipes and Anaerosporobacter were signi cantly increased at the PCL and LP infection stages compared with the non-infection group. Compared with the BD and PCL stages, we found that Ruminococcaceae_UCG-005, Bacteroides, and Tyzzerlla_4 were increased, while the Ruminococcaceae_NK4A214_ group, Famuly_XIII_AD3011_ group, Pseudobutyrivibrio, Saccharofermentans, and Lachnospiraceae_UCG-010 genera were sharply decreased. There were 31 genera signi cantly altered between PCL and BD stages in Co, Je, Re or Ru intestine. All information concerning changes of microbiota at different infection stages is listed in S11 Table. Associations of F. gigantica-induced immune responses with gut microbiota Based on the above statement, F. gigantica infection has a signi cant impact on the microbial community. However, the F. gigantica infection also induces Th2-biased and TLRs-related immune responses during infection as previously reported [8,9]. Therefore, we assumed that the microbiota community might be associated with Th2 and TLRs responses. Data for serum IgG and spleen cytokines (our previous reported data) were collected and re-examined together with the TLRs expression data obtained in this study to perform the CCA/RDA analysis. The results revealed that IgG1 appeared to be closely related to Oceanospirillales, Verrucomicrobiales, Burkholderiales, and Clostridiales orders ( Fig. 5A), among which Oceanospirillales was only changed between stage PCL and control. Verrucomicrobiales was altered between stages LP and BD, and between stages PCL and BD. Burkholderiales was found to be changed between stage PCL and control, stages LP and BD, and stage LP and control. Importantly, order Clostridiales was the most common and abundant order that underwent large changes during stages PCL, LP, and BD (S10 Table). Although Lachnospiraceae, Coriobacteriaceae, Family_XIII, Peptostreptococcaceae and Ruminococcaceae families belonging to order Clostridiales were signi cantly changed between infection stages (S6 Fig, S10 Table), the RDA analysis revealed that only the family Peptostreptococcaceae was closely associated with IgG1 (Fig. 5B). The Th2 cytokines IL4 and IL5 from the spleen were not closely linked to any bacterial orders, but at the family level, the spleen IL5 was closely linked to the Ruminococcaceae, while IL4 was closely related to Family_XIII ( Fig. 5B; S6 Fig).
We further investigated the associations between microbiota and TLRs innate responses. The CCA/RDA analysis revealed that at the order level, TLR1, 2, 3, 4, 6, 7 and 9, expressed in both HLN and spleen, appeared to be closely linked to the orders Bacteroidales, Burkholderiales, and Pseudomonadales, while TLR5 expressed in spleen was likely related to Clostridiales (Fig. 5C). At the family level, the TLRs 1, 2, 3, 4, 6, 7 and 9 in both HLN and spleen appeared to be closely associated with Lachnospiraceae, and Ruminococcaceae, while TLR5 in the spleen was likely related to Coriobacteriaceae (S6 Fig). At both order and family levels, the spleen TLR5 expression was opposite to that of TLR1, 2, 3, 4, 6, 7 and 9. In MN tissue, TLR2 and TLR4 appeared to be closely linked to order Methanobacteriales and family Methanobacteriaceae, but neither was signi cantly changed. TLR5 was tightly linked to order Oceanospirillales and family Prevotellaceae. TLR1 and 9 were likely related to Bacteroidales and Burkholderiales, respectively (Fig. 5D, S6 Fig). We also investigated the innate responses in liver (Fig. 5E)  To further examine the different functions of the microbes during the F. gigantica infection stages based on the changes in OTUs among all the samples, 12 samples from stages PCL (5 samples from Ru and Co), LP (4 samples from Ru, Re and Co) and BD (3 samples from Ru and Re) were selected to perform the metatranscriptome sequencing (Fig. 3, S2 Fig). NMDS analysis showed a good degree of difference between samples at both family and genus levels (stress < 0.2, S8 Fig). Gene prediction revealed that 61,125 to 231,996 genes in the selected samples were annotated, and the unigene CDS lengths were distributed from 200 bp to 1600 bp (Fig S8). All the unigene annotations are listed in S12 Table. Based on these annotations, KEGG analysis was applied to identify possible pathways involved in host immune responses. As shown in Fig. 6A, in line with the PICRUSt prediction, metatranscriptome sequencing showed that KEGG pathways involved in metabolism were the most abundant. These included carbohydrate metabolism, amino acid metabolism, lipid metabolism, and nucleotide metabolism. Among these metabolic pathways, energy metabolism, lipid metabolism, and carbohydrate metabolism have been reported to be associated with SCFAs production [42,43,44]. Moreover, SCFAs have been documented to be associated with host immune responses [45,46,47].
We then analyzed the KEGG pathway difference between infection stages. The results showed that energy metabolism, lipid metabolism, and carbohydrate metabolism pathways involving SCFAs production signi cantly differed between PCL and BD, and LP and BD. Although slightly decreased between stage PCL and control, and stage LP and control, these pathways were mainly decreased at stages BD and LP compared with stage PCL (S12 Table). Moreover, we found that several SCFAsproducing bacteria including Coriobacteriales, Lachnospiraceae, Prevotellaceae, Coriobacteriaceae, and Ruminococcaceae in the Re intestine were also altered between stages PCL and LP, and LP and BD (S10 and 11Tables). We investigated the production of SCFAs among infection stages by GC-MS methods, and the results showed that the main SCFAs including acetate, propionate, butyrate, and valerate were signi cantly decreased at the PCL, LP and BD infection stages in comparison to the control (Fig. 6B).
However, these SCFAs showed different pro les between different stages. Acetate was sharply increased at stage BD compared with stage PCL. Propionate was increased at stage LP compared with stage PCL but was decreased at stage BD. Butyrate decreased at stage LP compared with stage PCL and then increased at stage BD. The production of valerate showed a similar trend to that of butyrate (Fig. 6B). In line with the SCFAs pro ling, we found that the SCFAs-producing microbes Coriobacteriales, Coribacteriaceae, Pseudomonadales, Pseudomonadaceae, Enterobacteriales and Family_XIII decreased at stage PCL compared with control (S2B Fig). In addition, previous research revealed that SCFAs could regulate type 2 innate immune responses via histone deacetylase (HDAC) inhibition [48] and could induce pro-in ammatory cytokine production alone or in combination with toll-like receptor ligands [46,49]. We, therefore, supposed that the SCFAsproducing bacterial changes between stages PCL and LP, and LP and BD might affect the host Th2 responses via inhibition of HDAC. Our qPCR data showed that HDAC2 in spleen was signi cantly increased at stage PCL compared with controls and the other infection stages, while HDAC3 was increased at stage PCL and LP but was decreased at stage BD (Fig. 6C). These results imply that the SCFAs-producing bacteria may have released the inhibition of HDACs to mediate host Th2 responses. The Th2 response induced by F. gigantica has been supported by previous studies of Th2 cytokine expression in serum and PBMC [8]. These data suggest that SCFAs metabolism may regulate host Th2 responses by inhibiting HDACs.
Microbe-related TLR ligands may mediate the host Th2 response via TLR signaling As stated above, we found that TLRs in HLN, spleen, MN, and liver were associated with the alteration of gut microbes during the F. gigantica infection. Previous reviews delineated the roles of TLRs in immune activation and crosstalk between gut microbes and TLRs [50,51,52]. Generally, in the gut, various microbes carry types of TLR ligands such as lipoproteins, LPS, agellin and single-and double-stranded RNAs. These ligands are secreted into the gut and subsequently interact with MN or liver through the portal vein. These ligands can interact those TLRs [50] expressed on the surface of immune cells such as stromal cells, macrophages and hepatocytes to produce Th1-, Th2-or Treg-related cytokines and induce in ammatory or anti-in ammatory responses [51]. In this study, not only TLRs in the liver but also those TLRs in HLN, MN and spleen were related to speci c types of bacteria. These bacteria mainly belong to Bacteroidales, Burkholderiales and Pseudomonadales, Clostridiales, Oceanospirillales and Verrucomicrobiales, all of which carry many kinds of TLR ligands. We thus assumed that the microbes that were associated with the TLRs might affect host Th2 responses via secretion of TLR ligands by TLR signaling. To test this hypothesis, we analyzed the Th2-related cytokines IL4, IL13 and IL10 in MN tissues, where the TLR ligands from microbes in the gut may interact with TLRs. We found that IL4 and IL10 in MN were signi cantly increased at stage PCL compared with control, while IL13 was sharply increased at stage BD compared with control (data not shown). Simultaneously, the expression of TLRs in MN also increased dramatically at stage PCL (Fig. 2, S1 Fig). Similar results were observed previously in the liver.
For example, Zhang et al. [8]found that at 3 dpi (stage PCL) and 10 dpi (stage LP), the expression levels of TLR1, TLR2, TLR4, TLR6, TLR8, and TLR10 in the liver were increased, and Shi et al. [9] reported that IL4 expression in liver was increased at 28 dpi (stage LP), whereas IL13 was signi cantly increased at 28 dpi (stage LP) and 70 dpi (stage BD). More importantly, Shi et al. found that macrophages from buffalo liver secreted Th2-biased cytokines (cytokine production rather than mRNA expression) during F. gigantica infection (unpublished data). The ndings above together with the model from Miura K et al. [51] suggest that Th2 cytokines produced during F. gigantica infection may be induced by microbe-related TLR-ligands mediated TLR signaling.

Discussion
To our best knowledge, this is the rst report of the relationships between gut microbiota alteration, F. gigantica infection and host immune responses in buffaloes. Unlike gastrointestinal parasites, ukes do not reside in the gut tract for very long. Whether ukes would affect the host gut microbiota as in gut helminths was unknown. Previous reports concerning the relationships among gastrointestinal parasite infection, gut microbiota changes and immune responses (reviewed in ref. [53,54] may not be applicable to uke infection because of different parasite species and different host immunological environment [4,[55][56][57][58][59][60][61]. Moreover, previous reports mostly investigated such relationships using rodents that differ from the domestic animals that actually display immunological responses against uke infections. Previous studies indicated that helminth infection can change both the species abundance and composition of the gut microbiome [55,62]. Our results are consistent with these reports with signi cantly changed alpha and beta diversity indexes in infected buffaloes compared with controls. However, our ndings are opposite to the ndings in humans and primates that were naturally or experimentally infected with gut helminths, such as Necator americanus, Clonorchis sinensis, and Trichuris trichiura [63][64][65][66][67]. This might due to the differences in animal hosts, parasite species, infection stages, time of sampling or experimental setups (reviewed in ref [68]).
As F. gigantica parasites traverse the gut wall into the peritoneal cavity (stage PCL) and migrate into the liver capsule and parenchyma (stage LP) before they inhabit the bile ducts (stage BD), we supposed that at each stage these parasites would alter the host gut microbiome because of their close contact with host tissues and consequent immunological changes. This is evidenced by our results demonstrating that the gut microbes in F. gigantica-infected buffaloes were not only greatly changed at different stages but also in different gut sections in comparison to the control and between infection stages. For instance, in our study, most microbe changes happened at stages PCL and BD, implying that F. gigantica parasites affect buffalo gut microbes most when they traverse the gut wall into the peritoneal cavity and when they inhabit the bile ducts. This might be associated with some immunological or environmental factors such as Th1, Th2/Tregs balance, the gut pH value or nutrition [4,69]. In addition, our ndings demonstrated that bacterial abundances changed in different gut sections. At both family and genus levels, microbe abundance was altered in different gut sections. Till now, very few studies have investigated the changes in gut microbiota in different gut sections during the entire infection process. Because of the different experimental setups and parasite-animal models, we cannot yet reach a consensus conclusion. However, it is certain that ukes, despite not living for a long time in the host gut tract, do impact the gut microbiota during infection.
Recent studies have provided evidence that parasite-induced Th2 responses are essential to modi cations of the host gut microbiota, since such modi cations were lacking in Th2-de cient mice [4,[70][71][72]. Our previous studies demonstrated that a Th2-biased response was induced at the early stage of F. gigantica infection [8,9]. However, the relationships between the changes in gut microbiota at different stages and F. gigantica-induced Th2 responses in buffaloes were unclear. In this study, we revealed that gut microbiota changes are closely associated with Th2 responses during F. gigantica infection. CCA/RDA analysis revealed that the family Peptostreptococcaceae was strongly related to IgG1, the spleen IL5 was closely linked to the Ruminococcaceae, while IL4 was closely related to Family_XIII. This suggested that these microbes in the F. gigantica-infected buffalo gut are closely associated with Th2 immune responses. How microbes affect the F. gigantica-induced Th2 responses? Inspired by the previous study on the impact of SCFAs on the innate Th2 immune responses [48], we supposed the SCFAs might be the clue. In fact, we found signi cant decreases of butyrate, propionate, valerate, and acetate occur at PCL, LP and BD infection stages compared with control. In line with this pro ling, SCFAsproducing microbes were decreased at stage PCL in comparison to control. We also found that the decrease of such SCFAs-producing bacteria was related to the increase of Th2 responses. Moreover, our metagenome data also revealed that the KEGG pathways involved in SCFAs-production down-regulated during F. gigantica infection. These results suggest that F. gigantica infection may cause a decrease of SCFAs-producing microbes, which may subsequently induce a Th2-biased response via a certain pathway.
Previous research revealed that SCFAs could regulate type 2 innate immune responses via histone deacetylase (HDAC) [48]. Our results suggest that SCFAs might mediate host Th2 response via release of the inhibition of HDACs (mainly the spleen HDAC2 and HDAC3). This was consistent with Thio et al [48] who found that butyrate inhibited HDACs, which in turn decreased the production of Th2 cytokines and decreased GATA3 expression. However, we did not observe a signi cant increase of HDACs in HLN, liver or other tissues. More research concerning this topic is needed in the future. Previous studies have shown that immune responses might be modulated by SCFAs via activation of free fatty acid (FFA) receptors types 2 and 3 (FFAR2 and FFAR3) and G protein-coupled receptor 109A (GPR109A) (review in Ref. [73]).
We did not nd such activation of FFARs in the current study, implying that the Th2 responses might be associated with HADCs instead of FFARs, or perhaps the concentration of SCFAs produced by the SCFAsproducing microbes was not high enough to activate FFARs [74] in our study.
Besides Th2 responses, TLRs may also link with F. gigantica-induced microbiota changes. In our study, we found that F. gigantica-induced microbiota changes are closely related to TLR1, 2, 3, 4, 6, 7 and 9 at both order and family levels (Fig. 5C, D, E; S6 Fig). These related microbes carry many kinds of TLR ligands that recognize speci c TLRs and subsequently induce adaptive immune responses. Regarding the Th2 response, the most probable associated TLR is TLR2 as mentioned above. In our study, TLR2 mainly recognized the ligands from the orders Bacteroidales, Burkholderiales and Pseudomonadales, Verrucomicrobiales and Oceanospirillales. These bacteria were increased during infection stages compared with control, in line with the increase of TLR2 expression during infection. More importantly, the Th2-related cytokines IL4, IL13 and IL10 were signi cantly increased at stage PCL in MN and liver.
TLRs function as sensors for innate signals that can direct the adaptive immune response [75]. In the animal gut, TLRs recognize MAMPs that are shared by many types of microbes [10]. Previous studies demonstrated that activation of TLR3, 4, 5 and 9 skews T cell differentiation to the Th1 phenotype through IL12 produced by DCs, while activation of TLR2 induces a Th2-biased response via production of IL10 and IL13 [11,12,13]. TLR2 de ciency might result in an enhanced immunological phenotype bene cial for the induction of resistance during Schistosoma japonicum and Schistosoma mansoni infection ,showed decreased egg burdens, along with extended Th1 responses and signi cantly suppressed activation of Tregs cells [76,77]. Meanwhile, TLR2 and TLR2/4 KO mouse showed signi cantly in uence on oral microbial compositions, TLR2 and TLR2/4 KO mouse with increased Gramnegative species, while TLR4 KO mice up-regulated Gram-positive species [78]. These ndings suggest that TLRs function as multiplayers in the regulation of adaptive immune responses, and their corresponding ligands in certain conditions may react with certain TLRs and lead to the release of different cytokines that participate in different immunological processes. These results suggested that bacteria ligands were recognized by TLR2 and subsequently induced the production of Th2 cytokines, which in turn induced a Th2-biased immune response. This is consistent with previous reports concerning TLRs and the Th2 response [13,14].

Conclusions
In summary, our data revealed that buffalo gut microbes may affect the host Th2 responses during F. gigantica infection via via the SCFAs metabolism and TLR signaling pathway. Although such phenomena require further research to elucidate the detailed mechanisms involved, we suggest the mechanistic studies may better be carried out using F. gigantica-mouse models focusing on the early infection stage in which the Th2-biased immune response dominates considering the high cost of the buffalo. and Guangxi Science and Technology project (AA17204057, AB16380106, AB18221074, nycytxgxcxtd-09-05, . The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Availability of data and material
The metagenomic sequence data generated from this study can be found under European Nucleotide Archive (ENA) under the bio-project number PRJEB30859.

Ethics approval
The study design was reviewed and approved by the Animal Ethics Committee of Guangxi University, China.