The effect of helminth infection on the microbial composition and structure of the caprine abomasal microbiome

Haemonchus contortus is arguably the most injurious helminth parasite for small ruminants. We characterized the impact of H. contortus infection on the caprine abomasal microbiome. Fourteen parasite naive goats were inoculated with 5,000 H. contortus infective larvae and followed for 50 days. Six age-matched naïve goats served as uninfected controls. Reduced bodyweight gain and a significant increase in the abosamal pH was observed in infected goats compared to uninfected controls. Infection also increased the bacterial load while reducing the abundance of the Archaea in the abomasum but did not appear to affect microbial diversity. Nevertheless, the infection altered the abundance of approximately 19% of the 432 species-level operational taxonomic units (OTU) detected per sample. A total of 30 taxa displayed a significantly different abundance between control and infected goats. Furthermore, the infection resulted in a distinct difference in the microbiome structure. As many as 8 KEGG pathways were predicted to be significantly affected by infection. In addition, H. contortus-induced changes in butyrate producing bacteria could regulate mucosal inflammation and tissue repair. Our results provided insight into physiological consequences of helminth infection in small ruminants and could facilitate the development of novel control strategies to improve animal and human health.


Results
Haemonchus contortus infection changed the microbial habitat of the caprine abomasum. The primary infection by H. contortus in young goats appeared to negatively impact appetite and resulted in a decreased feed intake (Supplementary data). Consequently, infection led to a significantly reduced body weight gain (from 14.7 kg in uninfected control goats to 12.6 kg in infected goats) during the 50-day experiment period (Fig. 1). The infection significantly increased abomasal luminal pH to approximately 4.50 from a normal level of pH 2.93 in control uninfected goats (Fig. 1). Abomasal pH is known to affect parasite egg production when the abomasal pH reaches pH 4.0 to 4.5 9 . Thus, it is conceivable that parasite-induced changes in the microenvironment affect the survival and reproductive capacity of invading parasites.
Rarefaction curve analysis suggested that the sequencing depth in this study was adequate (Fig. 3), Common microbial diversity indices, such as Chao1 (Fig. 3), richness, Pielou's evenness, and Shannon and Simpson indices (Table 1), were evaluated. Our results suggested that H. contortus infection in goats did not appear to significantly affect species-level microbial diversity in the abomasal microbiome. For example, the Shannon index (mean ± sd = 4.10 ± 0.28 for control and 4.50 ± 0.88 for infected goats) did not differ significantly (P > 0.05 based on a modified t-test). In addition, species richness (91.96 ± 9.24 for control and 94.31 ± 18.21 for infected goats) and Phylogenetic Diversity (PD whole tree) (0.36.50 ± 3.45 for control and 36.41 ± 4.91 for infected) were indistinguishable between uninfected control and infected goats (P > 0.05).
The principal component analysis (PCA) results, however, showed a distinct difference in the microbial composition of the abomasal microbiome between uninfected control and infected goats (Fig. 4).
Haemonchus infection altered the microbial composition of the abomasal microbiome. The difference in the microbial composition (relative abundance) of the abomasal microbiome between goats in the uninfected control and infected groups was compared using the Linear Discriminant Analysis (LDA) Effect Size (LEfSe) algorithm 11 . 81 species-level OTU had a significant difference in relative abundance based on the Wilcoxon non-parametric t-test corrected for multiple hypothesis testing (P < 0.05). Of them, a total of 30 taxa displayed a significant difference in their abundance between goats in the uninfected control and infected groups at a stringent cutoff value (absolute LDA score log 10 ≥ 2.0). These taxa were depicted in Fig. 5. The relative abundance of bacteria was significantly higher in goats from the infected group, indicating that the total bacterial load in the abomasum was likely increased by infection. On the other hand, the abundance of the domain Archaea, exclusively from the phylum Euryarchaeota, was significantly higher in goats from the uninfected control group (Fig. 6). The relatively high abundance of Archaea in uninfected control goats was solely due to significantly higher abundance of uncultured species in a candidate family Methanomassiliicoccaceae. At the species level,  a total of 44 OTU had a significant difference in their abundance between the two groups (absolute LDA score log 10 ≥ 2.0). Of them, the 20 most abundant OTU were shown in Table 2. Eight of the 20 OTU contributed to the core abomasal microbiome. For example, Selenomonas ruminantium ( Fig. 6C) was approximately 37 fold more abundant in the abomasum of infected goats than uninfected control goats (LDA score 4.18).
Biological pathways and functional categories inferred from the 16S data. Microbial phylogeny and biological function are strongly linked. As a result, gene families and biological pathways can be predicted using 16S rRNA gene sequences. In this study, we analyzed the 16S rRNA gene sequence data to estimate the functional potential of H. contortus infection using the Phylogenetic Investigation of Communities by Reconstruction of Unobserved States (PICRUSt) algorithm 12 . After adjusting for copy number variations of 16S rRNA genes, a total of 5,404 KEGG gene families (mean ± sd = 4,283.1 ± 334.7 per sample) were predicted from the OTU table generated using the QIIME closed reference protocol. Of them, six gene families had significantly different abundance between the goats in the control and infected groups by LEfSE (absolute LDA log 10 scores > 2.0). For example, both α (KEGG#: K00174) and β (K00175) subunits of 2-oxoglutarate ferredoxin oxidoreductase had significantly higher abundance in the infected goats whereas asparagine synthase (K01953), peptide/nickel transport system ATP-binding protein (K02031), simple sugar transport system permease protein (K02057), and  methyl-galactoside transport system substrate-binding protein (K10540) were more abundant in the uninfected control group. The OTU contribution analysis suggests that the higher abundance of 2-oxoglutarate ferredoxin oxidoreductase subunits were likely due to the higher incidence of several OTU assigned to the genus Prevotella in the infected goats. Likewise, the predominance of simple sugar transport system permease protein (K02057) in goats from the control group may be attributed to the OTU belonging to Clostridia, especially those from the family Lachnospiraceae and the genus Butyrivibrio. Among the 30 major OTU contributors to the gene family K02057, 21 belonged to Clostridia, including at least six OTU assigned to Butyrivibrio. After classifying >5,400 predicted gene families into higher KEGG functional categories (pathways), LEfSE analysis identified 8 pathways displaying significantly different abundance between goats in the control and infected groups. As shown in Fig. 7, infection may have an effect on a broad range of biological functions, especially on ABC transporters, carbohydrate metabolism (TCA cycle), and amino acid metabolism (such as tyrosine metabolism). Intriguingly, hits assigned to biodegradation of naphthalene, one of the polycyclic aromatic hydrocarbons, were significantly higher in the infected goats, suggesting that H. contortus infection may have the potential to modify xenobiotics metabolism in the caprine gut.

Discussion
The barber's pole worm Haemonchus contortus is the most economically important parasitic nematode of small ruminants and remains the major determinant of profitability and sustainability in sheep and goat production for small farmers worldwide due to increased treatment cost and reduced productivity. Furthermore, the rapid emergence of anthelmintic drug resistant strains in this parasitic species is becoming a serious concern. The parasite elicits a strong immune response; and the first commercial vaccine is available to control the infection 13 . Our recent unpublished results suggested that infection of goats with H. contortus activated numerous biological pathways, including cytokine-mediated signaling and type I interferon signaling pathways (P < 10 −16 ) in the mucosa of the pyloric abomasum. Indeed, immune-mediated pathology rather than direct effects of the parasite itself may be directly responsible for clinical manifestations , such as reduced appetite, weight loss, and diarrhoea 3 .
It has been known that infection with helminth parasites induced a significant change in the structure and function of the host gut microbiome 7,14,15 . A 21-day infection of pigs with the whipworm Trichuris suis significantly altered approximately 13% of genera in the proximal colon microbiome and ~26% of metabolic pathways 14 . Furthermore, the effect of a primary infection with T. suis on the porcine proximal colon microbiome appeared to be long-lasting and independent of worm burden 15 . Similar effects on the gastrointestinal microbiome have been observed in hamsters infected with the liver fluke Opisthorchis viverrini 16 . A recent study showed that helminth-colonized humans had higher microbial species richness 17 . We performed LEfSe analysis to identify taxa displaying significant differences in abundance between goats infected with H. contortus and uninfected control goats. The mean number of species-level OTU identified in the caprine abomasum in this study was approximately 432 while 43 OTU were shared by the abomasal microbiome of all goats tested regardless of infection status. 81 OTU (i.e., ~19% of all OTU identified in individual abomasal microbiomes) were classified as significantly discriminative features before internal Wilcoxon test. Of them, 44 discriminative OTU were detected with an absolute LDA score (log 10 ) > 2. However, our results did not show any significant difference in various microbial diversity indices (P > 0.05; Table 1), including richness, Pielou's evenness, and Shannon, Simpson, and Brillouin indices at both the species and genus levels. A recent report also did not detect any effect of infection with the hookworm Necator americanus on community structure and microbial diversity in the human fecal microbiome 18 . Diet is known to be the most dominant determinant of the gut microbial composition. In our study, both uninfected control and H. contortus-infected goats were fed the same diet and were raised under the same environmental conditions. In contrast, a study of indigenous humans from Malaysia that were infected with  helminth parasites showed increased microbial species richness 17 but the diets of the two human cohorts studied were uncontrolled and could have confounded indicators of microbial diversity. While H. contortus infection in goats did not appear to affect microbial diversity at the site of infection in the abomasum, it did result in a significant alteration to the microbial composition of the abomasal microbiome, including altering the abundance of approximately 19% of all species-level OTU. These changes may be attributed to elevated pH in the microbial habitat because total viable bacteria have been shown to increase 2 to 3 orders of magnitudes to 10 8 -10 9 /ml in the abomasum of sheep as luminal pH rises 19 . Concurrent changes in abomasal anaerobic bacterial densities and pH values due to parasitic infection were also observed in sheep infected with T. circumcincta 8 . Our findings showed that H. contortus infection in goats tended to increase bacterial abundance with a concomitant reduction in Archaea abundance, which could have an important pathophysiological consequence, especially in host protein metabolism. While depressed appetite resulting from the infection likely reduced nutrient availability to the host and affected weigh gain, increased protein/nitrogen of microbial origin could compensate for this effect. More detailed nutritional determinations are needed to resolve this issue.
Haemonchus contortus infection increased the abundance of the genus Prevotella as well as the family Prevotellaceae. The mean abundance of Prevotella in the abomasal microbiome of the uninfected goats was 16.65% (± 5.47; sd) but increased to 25.35% (± 11.29; P < 0.05) in the infected goats. A total of 148 of the 1,663 OTU detected in this study were positively assigned to Prevotella, including P. ruminicola. While Prevotella abundance varies widely among cows, it appeared to be readily inducible and was the most abundant genus in the rumen 20,21 , especially in its liquid fraction 22,23 . In humans, the abundance of Prevotella was associated with higher carbohydrate or fiber-rich diets 24 . Reduced Prevotella incidence was observed in autistic children 25 as well as in obese Ossabaws pigs 26 . The species consisting of the genus Prevotella have versatile biological function, which is generally considered non-cellulolytic but saccharolytic 27 . In ruminants, Prevotella species play a crucial role in ruminal protein degradation, especially during oligopeptide breakdown, because they possess a rate-limiting dipeptidyl peptidase type IV (DPP-IV) activity, which is responsible for oligopeptides cleavage 28 . Furthermore, 14 of the 25 major OTU contributors to the gene family 2-oxoglutarate ferredoxin oxidoreductase β (K00175) were assigned to the genus Prevotella. 2-oxoglutarate ferredoxin oxidoreductase [EC:1.2.7.3] belongs to a family of oxidoreductases that participate in the citric acid (TCA) cycle, suggesting the H. contortus-altered abundance of Prevotella affects energy metabolism in the caprine gut microbial ecosystem. During helminth infection in ruminants, feed intake is significantly reduced while endogenous protein loss accelerates due, in part, to the processes that enhance mucin excretion and the synthesis of effector molecules for host immune responses 29,30 . As a result, protein availability and quality in the intestine play a larger role as the infection progresses and serves as a main determinant of host nitrogen utilization efficiency 31 . The increased abundance of Prevotella in H. contortus-infected goats would have functional importance in host protein metabolism. Indeed, Prevotella is more abundant in nutrient-inefficient bulls and associated with positive residual feed intake 21 .
The pathophysiological relevance of parasitic infection-induced changes in microbial composition of the gut microbiome is still unclear. It is conceivable that both passive and active mechanisms may be involved. Parasitic infections disrupt microbial habitats. In the caprine abomasum, luminal pH values are significantly higher due to the infection. As a result, survival and proliferation of certain microbial species becomes favored 32 . For example, a substantially elevated level of Mucispirillum was observed in the luminal contents of the proximal colon of pigs infected by T. suis (whipworms) 14 . Mucispirillum species colonize the intestinal mucus and have a predilection for mucosal surfaces. Mucosal disruption induced by invading T. suis was likely responsible for dislodging mucus native microbial species to the gut lumen. More importantly, the alteration to the gut microbiome induced by parasitic infection is part of host immune responses 7 and may represent an adaptive mechanism. In this study, sequences assigned to Pasteurellales (and the family Pasteurellaceae), represented by a group of pathogenic Gram-negative bacteria, were absent in the uninfected control goats but significantly increased in H. contortus-infected goats (Fig. 6B). A well studied species from this group is Pasteurella multocida, which is associated with the clinical syndromes defining bovine respiratory diseases (BRD) complex, including pneumonia. P. multocida-induced BRD symptoms are related to environmental stresses as well as concurrent viral or bacterial infections 33 . Our data provided further evidence that parasitic infections increased the risk of subsequent bacterial infections. Furthermore, as mentioned above, altered Prevotella abundance may be involved in the compensation by the host for its impaired protein metabolism. Parasite infection-induced changes in the composition of the gut microbiome could explain special nutritional needs during the infection that include requirements for the sulfur-containing amino acid methionine and other essential amino acids. Furthermore, it has long been known that helminth parasites have evolved abilities to dampen inflammation for their survival in the host 34 . However, the mechanism by which parasites evade host immune surveillance remains elusive. We hypothesize that helminth infection modulates the gut butyrate biosynthesis by altering the abundance of butyrate-producing bacteria. Butyrate is a potent inhibitor of inflammation 35 . In this study, we observed a significant difference between uninfected control and H. contortus-infected goats in the relative abundance of the genus Butyrivibrio, known to harbor butyrate-producing bacteria. At the species level, at least two OTU positively assigned to the genus Butyrivibrio, GreenGene ID# 292105 and #847896, had 84.5 and 4.1 fold higher abundance in goats from the uninfected control than H. contortus-infected groups, respectively (LDA scores > 2.0). Motivated by this observation, we plan to examine the effect of parasitic infection on the faction of butyrate-producing bacteria of the gut microbiome in host-parasite models and their potential role in regulating mucosal inflammation.
The abomasum of ruminants shares similar structure and function with the human stomach. Gastric pH is essential for protein digestion and absorption of nutrients like calcium, iron, and vitamin B12, and gastric acidity also plays a critical role in the pathogenesis of diseases, such as abomasal ulceration, abomasitis, abomasal bloating, and gastric tumors. Low acidic environment in the abomasum functions as a potent barrier against bacterial infection and growth. However, uncontrolled proliferation of bacteria also results in several abomasal diseases in ruminants 36,37 . Previous studies suggested that the number of microbes in the human stomach was related Scientific RepoRts | 6:20606 | DOI: 10.1038/srep20606 to luminal pH changes induced by diseases and therapy 38 . An in vitro model has been established to simulate the effect of pH on gastric bacteria 32 . Nevertheless, the goat-H. contortus infection model could offer a unique opportunity to understand the pathogenesis of human gastric disorders that are affected by changes in pH and the related effects on microbial communities and function.

Methods
Animals and parasitology. Twenty male Alpine dairy kid goats were obtained locally within 72 h after birth and were subsequently raised indoors on a concrete floor throughout the experiment. All goats were housed together in the same building separately by pens. The goat kids were fed milk replacer with 25% protein and 28% fat until weaning at approximately two months of age. After weaning, the young goats were fed brome grass hay ad libitum and had free access to water. When the goats reached approximately three months of age, 14 randomly selected goats were orally infected with a single dose of 5,000 H. contortus infective third-stage larvae (L3) and maintained for 50 days (primary infection). Six goats were not infected and remained parasite naive throughout the experiment and served as uninfected controls. At 50 days post infection (dpi), both infected and uninfected control goats were sacrificed. The live bodyweight of each animal was recorded immediately prior to the primary infection and prior to necropsy at 50 dpi. Fecal egg counts were monitored four times during the experiment, immediately prior to the experimental inoculation, and at 4 and 6 weeks of the primary infection, and immediately prior to necropsy at 50 dpi, as previously described 39 . Adult worms as well as immature larvae from both contents and the tissue of the abomasum were isolated and counted 39 . Abomasal pH values were measured using a hand-held pH meter. The luminal contents of the abomasum were then sampled and snap frozen in liquid nitrogen prior to storage at − 80 °C until total DNA was extracted. During the experiment, animals were handled according to a protocol (Protocol #12-025) approved by the BARC Animal Care and Use Committee; and Institutional Animal Care and Use Committee (IACUC) guidelines were strictly followed. The experimental procedures were carried out in accordance with the approved guidelines.
DNA extraction, 16S rRNA gene amplicon preparation, and sequencing. Total DNA was extracted from the abomasal luminal contents 14 . Briefly, DNA was extracted using a QIAamp DNA stool kit (Qiagen, Valenica, CA) with some modifications. An eight-minute incubation at 95 °C was used to replace the 70 °C lysis recommended in the standard protocol. DNA integrity was verified using a BioAnalyzer 2000 (Agilent, Palo Alto, CA). DNA concentration was then quantified using a QuantiFluor fluorometer (Promega, Madison, WI). The hypervariable V3 -V4 regions (E. coli position 341 to 805) of the 16S rRNA gene, which has been frequently targeted for interrogating bacterial communities 40 , were directly amplified from 10 ng of total DNA with PAGE-purified Illumina platform-compatible adaptor oligos that contain features such as sequencing primers, sample-specific barcodes, and 16S PCR primers (forward primer, 341/357F: NNNNCCTACGGGNGGCWGCAG; reverse primer, 805/785R: GACTACHVGGGTATCTAATCC). The PCR reaction included 2.5 units of AccuPrime Taq DNA Polymerase High Fidelity (Invitrogen, Carlsbad, CA) in a 50-μ l reaction buffer containing 200 nM primers, 200 nM dNTP, 60 mM Tris-SO 4 , 18 mM (NH4) 2 SO 4 , 2.0 mM MgSO4, 1% glycerol, and 100 ng/ul bovine serum albumin (New England BioLabs, Ipswich, MA). PCR was performed using the following cycling profile: initial denaturing at 95 °C for two min followed by 20 cycles of 95 °C 30 sec, 60 °C 30 s, and 72 °C 6 s. Amplicons were purified using a Agencourt AMPure XP bead kit (Beckman Coulter Genomics, Danvers, MA) and quantified using a BioAnalyzer high-sensitivity DNA chip and a QuantiFluor fluorometer. The purified amplicons from individual samples were pooled in equal mass (molar) ratios. The purified amplicon pool was further spiked with approximately 25% of whole-genome shotgun libraries prepared using an Illumina TruSeq DNA sample prep kit with a compatible adaptor barcode to enhance sequence diversity during the first few cycles of sequencing for better cluster differentiation. The concentration of the final library pool was quantified using a BioAnalyzer High-sensitivity DNA chip kit (Agilent). The library pool was sequenced using an Illumina MiSeq Reagent Kit on an Illumina MiSeq sequencer. Raw sequence data are deposited to the MG-RAST server (http://metagenomics. anl.gov) and available to the public (ProjectID# 13390; MG-RAST ID 4629311.3 to 4629350.3).
Bioinformatics and data analysis. The data were preprocessed using MiSeq Control Software (MCS) v2.4.1. Raw sequences were first analyzed using FastQC version 0.11.2 to check basic statistics, such as GC%, per base quality score distribution, and sequences flagged as poor quality. The four maximally degenerate bases ("NNNN") at the most 5′ end of Read#1 of the read pair, which were designed to maximize the diversity during the first four bases of the sequencing run for better identification of unique clusters and improve base-calling accuracy, were then removed. The presence of forward and reverse PCR primers at the 5′ and 3′ ends of each sequence read was scanned; the reads without primers were discarded. Chimeric reads were also removed. Approximately 96.9% of raw reads were retained after these initial quality control steps. The processed pair-end reads were then merged using PandaSeq v2.8 to generate representative complete nucleotide sequences (contigs) using default parameters 41 . The overlapping regions of the pair-end read were first aligned and scored; and reads with low score alignments and high rate of mismatches were discarded. As a result, high-quality contigs can be formed from approximately 92% of raw read pairs (the mean number of contigs per sample ± sd = 208,241.00 ± 51,142.79; N = 20). Contig sequences were initially analyzed using BLAT against the SILVA database 42 . Principal component analysis (PCA) was performed using the ade4 package in R as previously described 15 . QIIME pipeline (v1.9.1) was used to analyze the 16S rRNA gene sequences 10 . A "closed-reference" protocol in the pipeline was used for OTU picking. The default QIIME parameters were used, except that the quality-filtering based on OTU abundance threshold was lowered to 0.0001%. The latest version of GreenGene database (v13.8) was used for taxonomy assignment (greengenes.lbl.gov). PyNAST (v1.2.2) was used for sequence alignment. In addition, the greedy heuristic clustering algorithm CD-HIT-OTU 43 (v4.5.5) was used to pick OTU for the comparison purpose (Supplementary data). Alpha (α )-diversity at the species-level was estimated using QIIME and Primer programs. Beta diversity was also estimated using PCA in R. The OTU relative abundance values were analyzed using the LEfSe algorithm 11 to identify taxa and KEGG gene families and/or pathways that display significant differences between two biological conditions. The algorithm first used the non-parametric factorial Kruskal-Wallis sum-rank test to detect taxa with significantly different abundance with respect to the treatment class, followed by unpaired Wilcoxon rank-sum test to detect biological consistency between uninfected and H. contortus-infected goats, and then used LDA to estimate the effect size of each differentially abundant feature. Furthermore, PICRUSt (v1.0.0) 12 , a software package designed to predict metagenome functional contents from marker gene (e.g., the 16S rRNA gene) surveys, was used with default parameters to predict gene contents and metagenomic functional information based on the OTU table generated using the closed-reference protocol in QIIME. Briefly, the OTU table was first normalized by dividing each OTU by the known/predicted 16S copy number by using the PICRUSt workflow normalize_by_copy_number.py. The gene contents or the abundance of KEGG Orthologs were predicted from the normalized OTU table using the workflow predict_metagenomes. py. The predicted metagenome function was further analyzed by collapsing thousands of KEGG Orthologs into higher functional categories (pathways) (categorize_by_function.py). In addition, specific OTU contributing to a given function or pathway was identified by using the workflow metagenome_contributions.py.