Strongyle Infection and Gut Microbiota: Profiling of Resistant and Susceptible Horses Over a Grazing Season

Gastrointestinal strongyles are a major threat to horses' health and welfare. Given that strongyles inhabit the same niche as the gut microbiota, they may interact with each other. These beneficial or detrimental interactions are unknown in horses and could partly explain contrasted susceptibility to infection between individuals. To address these questions, an experimental pasture trial with 20 worm-free female Welsh ponies (10 susceptible (S) and 10 resistant (R) to parasite infection) was implemented for 5 months. Fecal egg counts (FEC), hematological and biochemical data, body weight and gut microbiological composition were studied in each individual after 0, 24, 43, 92 and 132 grazing days. R and S ponies displayed divergent immunological profiles and slight differences in microbiological composition under worm-free conditions. After exposure to natural infection, the predicted R ponies exhibited lower FEC after 92 and 132 grazing days, and maintained higher levels of circulating monocytes and eosinophils, while lymphocytosis persisted in S ponies. Although the overall gut microbiota diversity and structure remained similar during the parasite infection between the two groups, S ponies exhibited a reduction of bacteria such as Ruminococcus, Clostridium XIVa and members of the Lachnospiraceae family, which may have promoted a disruption of mucosal homeostasis at day 92. In line with this hypothesis, an increase in pathobionts such as Pseudomonas and Campylobacter together with changes in several predicted immunological pathways, including pathogen sensing, lipid metabolism, and activation of signal transduction that are critical for the regulation of immune system and energy homeostasis were observed in S relative to R ponies. Moreover, S ponies displayed an increase in protozoan concentrations at day 92, suggesting that strongyles and protozoa may contribute to each other's success in the equine intestines. It could also be that S individuals favor the increase of these carbohydrate-degrading microorganisms to enhance the supply of nutrients needed to fight strongyle infection. Overall, this study provides a foundation to better understand the mechanisms that underpin the relationship between equines and natural strongyle infection. The profiling of horse immune response and gut microbiota should contribute to the development of novel biomarkers for strongyle infection.


INTRODUCTION
Grazing horses are infected by a complex community of parasitic helminths, mainly strongyles (Ogbourne, 1976;Bucknell et al., 1995;Corning, 2009;Kuzmina et al., 2016). Like other worms, strongyles have a direct life cycle in which they can survive outside of its host in pastures as well as in the horse's intestines (Taylor et al., 2007). Infective larvae (L3 stage) are usually ingested and migrate to their preferred niche in the small or large intestine. After two molts, they will eventually become sexually mature adults and will lay eggs that are passed onto the pasture via feces (Taylor et al., 2007).
Equine strongyle species are classified into Strongylinae and Cyathostominae, which differ, among other criteria, by their respective size (Lichtenfels et al., 2008). Strongylus vulgaris is the most pathogenic of the large strongyles as a result of the intestinal infarction that larval stages can cause during their migration. S. vulgaris prevalence has been drastically reduced since the release of modern anthelmintics (Nielsen et al., 2012). In contrast, nearly all horses are infected by small strongyles or cyathostomins throughout the world (Ogbourne, 1976;Bucknell et al., 1995;Lyons et al., 1999). Small strongyles infections are responsible for milder symptoms, including weight loss or poor hair condition (Love et al., 1999). Infections are more common in young animals (Love et al., 1999) largely due to immunological hyporesponsiveness, meaning there is likely an immune component to infection susceptibility (Lyons et al., 1999). In fact, the prevalence of infections has been found to be lower in horses older than 5 years (Relf et al., 2013;Kornaś et al., 2015). However, larval stages encyst into the colonic mucosa as part of their life cycle and their massive emergence (larval cyathostominosis) can result in abdominal pain, diarrhea (Love et al., 1999;Matthews, 2014) and eventually death (Giles et al., 1985).
Drug inefficacy reports have accumulated worldwide over the recent years, and resistant cyathostomin isolates are now to be found in Europe (Geurden et al., 2014;Sallé et al., 2017), America , and Oceania (Scott et al., 2015). It is therefore required to alleviate drug selection pressure put on cyathostomin populations. Targeted-selective treatment can help reduce drug usage as horses differ in their intrinsic potential to resist infection. Indeed, it has been estimated that 80% of the total worm burden is produced by 20% of horses (Wood et al., 2012) and that 21% of the inter-individual variation has a heritable component (Kornaś et al., 2015). However, the factors underpinning this phenotypic contrast (Lester et al., 2013;Sallé et al., 2015) still remain unclear.
Equine strongyles are in close contact with a large community of microorganisms in the host intestines, estimated to reach a concentration of 10 9 microorganisms per gram of ingesta in the cecum alone (Mackie and Wilkins, 1988), spanning 108 bacterial genera (Steelman et al., 2012;Venable et al., 2017) and at least seven phyla (Costa et al., 2012(Costa et al., , 2015Shepherd et al., 2012;Weese et al., 2015;. Bacterial populations differ greatly throughout the various compartments of the equine gastrointestinal tract (e.g., duodenum, jejunum, ileum and colon) due to differences in the gut pH, available energy sources, epithelial architecture of each region, oxygen levels and physiological roles (Costa et al., 2015;Ericsson et al., 2016). The gut microbiota promotes digestion and nutrient absorption for host energy production and synthesizes folate (Sugahara et al., 2015), vitamin K 2 (Marley et al., 1986) and short chain fatty acids (SCFA) such as acetate, butyrate and propionate (Nedjadi et al., 2014;Ericsson et al., 2016). The gut microbiota also protects the host from pathogens and stimulates and matures the immune system and epithelial cells (reviewed by Nicholson et al., 2012). Although bacteria represent the major portion of the hindgut microbiota in horses (Dougal et al., 2013;Fernandes et al., 2014), the quantity of anaerobic fungi, protozoa and archaea relative to total bacteria found in horse feces is about 0.26, 3.7, and 0.0015, respectively (Dougal et al., 2012).
The physical presence of helminths in the intestinal lumen can alter the gut microbiota activity and composition (Midha et al., 2017;Peachey et al., 2017). These perturbations have been studied in various host-nematode relationships including rodent models (Reynolds et al., 2014;Fricke et al., 2015;McKenney et al., 2015;Su et al., 2017), ruminants (Li et al., 2011(Li et al., , 2016El-Ashram and Suo, 2017), pigs Wu et al., 2012), and cats (Duarte et al., 2016). However, it remains unresolved whether nematode infection has a beneficial (Lee et al., 2014) or a detrimental (Houlden et al., 2015) impact on the gut microbiota diversity, richness and functions. The mechanisms responsible for these gut microbiota alterations are also unclear and could arise indirectly due to the immune response helminths trigger in their host (Reynolds et al., 2014(Reynolds et al., , 2015Fricke et al., 2015;Cattadori et al., 2016;Zaiss and Harris, 2016), such as regulatory T cell stimulation and lymphoid tissue modifications and changes in the intestinal barrier (Boyett and Hsieh, 2014;Giacomin et al., 2016). On the other hand, the secretion of putative anti-bacterial compounds (Mcmurdie and Holmes, 2012;Holm et al., 2015) or modifications in the intestinal environment that supports helminth survival might directly modify the gut microbiota (D'Elia et al., 2009;Midha et al., 2017).
The putative interactions between equine strongyle infection, the gut microbiota and host physiology are unknown in horses. To address these questions, grazing ponies with an extreme phenotypic resistance or susceptibility toward natural strongyle infection were monitored over a 5-month period. We aimed to provide insights into the host response and the gut microbiological composition associated with strongyle natural infection in order to guide the development of new microbiotabased control strategies.

Selection of Animals
Twenty female Welsh ponies (10 resistant, R, and 10 susceptible, S; 5 ± 1.3 years old) were selected from an experimental herd of 107 ponies (INRA, UEPAO, Nouzilly, France) regularly monitored for fecal egg counts (FEC) during the grazing season since 2011. Selection was based on the FEC records database consisting of 753 observations (n = 146, 278, 214, and 115 for spring, summer, autumn, and winter measures) of the whole herd. Every individual pony had been recorded at least three times over a minimum of 2 years (seven observations per pony on average across the whole herd). FEC data were log-transformed to correct for over-dispersion and fitted a linear mixed model accounting for environmental fixed effects (month of sampling, year of sampling, time since last treatment, age at sampling). The individual was considered as a random variable to account for the intrinsic pony potential against strongyle infection. Estimated individual effects were centered and reduced to express each individual potential as a deviation from the mean on the logarithmic scale. Based on these estimates, the 20 individuals with most extreme estimates were chosen to create a resistant R group (mean individual effect with −1.18 phenotypic deviation from the mean and average age of 5.6 years) and a susceptible S group (mean individual effect with +1.45 phenotypic deviation from the mean and average age of 4.7 years). The median egg counts /g feces across the past 5 years were 800 for S (the mean was 897, the first quartile was 350, and the third quartile was 1,375), whereas the median egg counts/g feces were 0 for R (the mean was 24.6, the first quartile was 0, and the third quartile was 150; Figures S1A,B).
As it has been previously established that gut microbiota profiles might be shaped by host genetics (Lozupone et al., 2012;Goodrich et al., 2014), this information was considered to understand the individual variance underlying microbiota composition. The "kinship2" package in R was used to create the genetic relationship matrix. Both pedigree tree and correlation structure matrix are depicted in Figures S1C,D, respectively.
All the procedures were conducted according to the guidelines for the care and use of experimental animals established by the French Ministry of Teaching and Research and the regional Val de Loire Ethics Committee (CEEA VdL, n • 19). The protocol was registered under the number 2015021210238289_v4 in the experimental installations with the permit number: C371753. All the protocols were conducted in accordance with EEC regulation (n • 2010/63/UE) governing the care and use of laboratory animals and effective in France since the 1st of January 2013.
Grazing started by the end of a 3-month moxidectin residual period (mid-June 2015). Ponies grazed from mid-June to the end of October 2015 at the Nouzilly experimental station (France). The experimental pasture (7.44 ha) consisted of Festuca arundinacea, Phleum prateonse, Poa abbreviate, Holcus lanatus, and Dactylis glomerata. During all phases of the experimental period, ponies were provided access to water ad libitum.
Moxidectin is known to have moderate efficacy against encysted cyathostomin larvae (Xiao et al., 1994;Reinemeyer et al., 2015). Therefore, egg excretion can occur at the end of the residual period, as a result of encysted larvae completing their development into adults. To eliminate this residual excretion and to avoid any interference with strongyle infection at pasture, a treatment targeting the only luminal immature and adult stages (Strongid R paste, Zoetis, Paris, France; single oral dose of 1.36 mg pyrantel base per kg of body weight) was implemented at day 30.
Fecal samples were collected from the rectum. Fecal aliquots for microbiota analysis were immediately snap-frozen in liquid nitrogen and stored at −80 • C until DNA extraction, whereas fecal aliquots to measure the fecal egg counts were immediately sent to the laboratory.
The pH in the feces was determined after 10% fecal suspension (wt/vol) in saline solution.
FIGURE 1 | Experimental design and sampling. A set of 20 female ponies [10 susceptible (S) and 10 resistant (R) to strongylosis] were selected based on their fecal egg counts history during previous pasture seasons and were kept inside during the winter. In the spring, they were treated with moxidectin, to ensure that they were totally free from gastrointestinal nematodes (even from putative encysted larvae) and were kept indoors for 3 months. Thereafter, once the effect of moxidectin treatment was no longer detected, the ponies were moved to a 7.44 ha pasture to start the study. At day 30 of the study, a pyrantel treatment was administered to all the animals in order to eliminate any residual infections that could interfere in the protocol. A longitudinal monitoring of the parasitism level in each animal was performed through five time points from June to October. At each time point, fecal samples were collected from all ponies on 0, 24, 43, 92, and 132 days after the beginning of the grazing season to carry out fecal egg counts, pH measurements, and microbiota profiling. Blood samples were taken at the same time points to analyze biochemical and hematological parameters. Body weight (BW) was also recorded at the same time points.
Blood samples were taken from each pony and collected in EDTA-K3-coated tubes (5 mL) to determine hematological parameters and heparin tubes (10 mL) to determine biochemical parameters. After clotting, the heparin tubes were centrifuged at 4,000 rpm during 15 min and the harvested plasma was stored at −20 • C until analysis. Additionally, blood collected in EDTA-K3tubes was used to measure the different blood cells.
For each individual, body weight was measured monthly and average daily weight gain was also calculated.
None of the ponies received antibiotic therapy during the sampling period and diarrhea was not detected in any ponies.

Fecal Egg Counts
Fecal egg counts (eggs per gram of wet feces) was measured as a proxy for patent strongyle infection. FEC was carried out using a modified McMaster technique (Raynaud, 1970) on 5 g of feces diluted in 70 mL of NaCl solution with a density of 1.2 (sensitivity of 50 eggs/g). A Wilcoxon rank-sum test with Benjamini-Hochberg multiple test correction was used to determine whether there was a significant difference between groups across the experiment. An adjusted p < 0.05 was considered significant.

Blood Hematological and Biochemical Assays
For blood hematological assays, blood was stirred for 15 min at room temperature to facilitate oxygenation. Different blood cells were analyzed, including leucocytes (lymphocytes, monocytes, neutrophils, basophils and eosinophils), erythrocytes and different blood parameters such as hematocrit, mean corpuscular volume and the thrombocytes. The total blood cells were counted with a MS9-5 Hematology Counter R (digital automatic hematology analyzer, Melet Schloesing Laboratories, France).
The serum biochemical parameters (albumin, cholesterol, globin, glucose, phosphatase alkaline, total proteins and urea) were measured with a colorimetric method using Select-6V rings with the M-Scan II Biochemical analyzer (Melet Schloesing Laboratories, France).
Mixed-effects analysis of the variance (ANOVA) or Wilcoxon rank-sum tests were conducted for continuous variables fitting a normal or non-normal distribution respectively to delineate whether there was a significant difference between the average values of phenotype traits for the different groups, using a significance level of p < 0.05. Blood cell counts were corrected for the mild dehydration by fitting the hematocrit as a co-variable in the model.

Weather Data
Daily precipitation and temperatures were recorded at a meteorological station located 14 km from the experimental area.

Pasture Contamination Analysis
Pasture contamination, expressed as the number of infective larvae (L3) per kg of dry herbage, was measured throughout the grazing season. At every occasion, four plucked grass were picked up in 100 random locations across pasture (Gruner and Raynaud, 1980). A sub-sampling of 600 g of collected grass was mixed with 20 mL of neutral pH soap diluted in 5 L of tap water and left at room temperature for 24 h. Grass was then separated using a gross sieve (0.5 cm mesh) and discarded. The 5 L water was then passed through both a 125 µm mesh sieve and a 20 µm mesh sieve to retain smaller particles, including strongyle larvae. Recovered particles were diluted in 40 mL of tap water and subsequently dispatched into 4 × 10 mL glass tubes. A flotation technique was used to isolate larvae from other abiotic contaminants as follows. Tubes were centrifuged at 2,500 rpm for 5 min, before discarding supernatant and adding a dense solution (NaCl, density 1.18-1.2) and a cover slip onto glass tube. Gentle centrifuge was subsequently performed at 1,500 rpm for 8 min thus allowing larval material adherence onto the cover slip that was further examined under an optical microscope. This last step was performed three times.

Microorganisms DNA Extraction From Feces Samples
Total DNA was extracted from aliquots of frozen fecal samples (200 mg; 100 samples at different time points from 20 ponies), using E.Z.N.A. R Stool DNA Kit (Omega Bio-Tek, Norcross, Georgia, USA). The DNA extraction protocol was carried out according to the manufacturer's instructions (Omega-Bio-Tek, Norcross, Georgia, USA).

V3-V4 16S rRNA Gene Amplification
The V3-V4 hyper-variable regions of the 16S rDNA gene were amplified with two rounds of PCR using the forward primer (5 ′ -CTTTCCCTACACGACGCTCTTCCGATCTACGGRAG GCAGCAG-3 ′ ) and the reverse primer (5 ′ -GGAGTTCAGAC GTGTGCTCTTCCGATCTTACCAGGGTATCTAATCCT-3 ′ ) modified in order to include Illumina adapters and barcode sequences which allow for directional sequencing. The targeted region resulted in an amplicon of a region of ∼459 bp that when sequenced with paired-end reads had at least ∼50 bp of overlapping sequence in the middle. The first round of amplification was performed in triplicate in a total volume of 50 µL containing 10 ng of DNA, 2.5 units of a DNA-free Taq DNA Polymerase and 10X Taq DNA polymerase buffer (MTP Taq DNA Polymerase, Sigma). Subsequently, 10 mM of dNTP mixture (Euromedex, Souffelweyersheim, France), 20 mM of each primer (Sigma, Lezennes, France) and Nuclease-free water (Ambion, Thermo Fisher Scientific, Waltham, USA) were added. Ultrapure Taq DNA polymerase, ultrapure reagents, and plastic were selected in order to be DNA-free. The thermal cycle consisted of an initial denaturation step (1 min at 94 • C), followed by 30 cycles of denaturation (1 min at 94 • C), annealing (1 min at 65 • C) and 1 min of extension at 72 • C. The final extension step was performed for 10 min at 72 • C. Amplicons were then purified using magnetic beads (Clean PCR system, CleanNA, Alphen an den Rijn, The Netherlands) as follows: beads/PCR reactional volume ratio of 0.8X and final elution volume of 32 µL using Elution Buffer EB (Qiagen). The concentrations of the purified amplicons were checked using a NanoDrop 8000 spectrophotometer (Thermo Fisher Scientific, Waltham, USA).
Sample multiplexing was performed using 6 bp unique indexes, which were added during the second PCR step at the same time as the second part of the P5/P7 adapters used for the sequencing step on the Illumina MiSeq flow cells with the forward primer (5 ′ -AATGATACGGCGACCA CCGAGATCTACACTCTTTCCCTACACGAC-3 ′ ) and reverse primer (5 ′ -CAAGCAGAAGACGGCATACGAGATNNNNNN GTGACTGGAGTTCAGACGTGT-3 ′ ).
This second PCR step was performed using 10 ng of purified amplicons from the first PCR and adding 2.5 units of a DNA-free Taq DNA Polymerase and 10X MTP TaqDNA polymerase buffer (Sigma). The buffer was complemented with 10 mM of dNTP mixture (Euromedex), 20 mM of each primer (Eurogentec, HPLC grade) and Nuclease-free water (Ambion, Life Technologies) up to a final volume of 50 µL. The PCR reaction was carried out as follows: an initial denaturation step (94 • C for 1 min), 12 cycles of amplification (94 • C for 1 min, 65 • C for 1 min and 72 • C for 1 min) and a final extension step at 72 • C for 10 min. Amplicons were purified as described for the first PCR round. The concentration of the purified amplicons was measured using Nanodrop 8,000 spectrophotometer (Thermo Scientific) and the quality of a set of amplicons (12 samples per sequencing run) was checked using DNA 7,500 chips onto a Bioanalyzer 2,100 (Agilent Technologies, Santa Clara, CA, USA). All libraries were pooled at equimolar concentration in order to generate equivalent number of raw reads with each library. The final pool had a diluted concentration of 5 nM to 20 nM and was used for sequencing. Amplicon libraries were mixed with 15% PhiX control according to the Illumina's protocol. Details on sequencing, PhiX control and FastQ files generation are specified elsewhere (Lluch et al., 2015). For this study, one sequencing run was performed using MiSeq 500 cycle reagent kit v2 (2 × 250 output; Illumina, USA).

Sequencing Data Preprocessing
Sequences were processed using the version 1.9.0 of the Quantitative Insights Into Microbial Ecology (QIIME) pipeline (Caporaso et al., 2010;Rideout et al., 2014) and by choosing the open-reference operational taxonomic units (OTU) calling approach (Rideout et al., 2014).
First, forward and reverse paired-end sequence reads were collapsed into a single continuous sequence according to the "fastq-join" option of the "join_paired_ends.py" command in QIIME. The fastq-join function allowed a maximum difference within overlap region of 8%, a minimum overlap setting of 6 bp and a maximum overlap setting of 60 bp. The reads that did not overlap (∼20% of the total) were removed from the analysis. The retained sequences were then quality filtered. De-multiplexing, primer removal and quality filtering processes were performed using the "split_libraries"_fastq.py command in QIIME. We applied a default base call Phred threshold of 20, allowing maximum three low-quality base calls before truncating a read, including only reads with >75% consecutive high-quality base calls, and excluding reads with ambiguous (N) base calls (Navas-Molina et al., 2013).
Subsequently, the sequences were clustered into OTUs against the GreenGenes database (release 2013-08: gg_13_8_otus) (DeSantis et al., 2006) by using the uclust (Edgar, 2010) method at a 97% similarity cutoff. The filtering of chimeric OTUs was performed by using Usearch version 6.1 (Edgar et al., 2011) against the GreenGenes reference alignment (DeSantis et al., 2006). A phylogenic tree was generated from the filtered alignment using FastTree (Price et al., 2010). Singletons were discarded from the dataset to minimize the effect of spurious, low abundance sequences using the "filter_otus_from_otu_table.py" script in QIIME. To confirm the annotation, the resulting OTU representative sequences were then searched against the Ribosomal Database Project naïve Bayesian classifier (RDP 10 database, version 6; Cole et al., 2009) database, using the online program SEQMATCH (http://rdp.cme.msu.edu/seqmatch/seqmatch_intro.jsp). Finally, consensus taxonomy was provided for each OTU based on the taxonomic assignment of individual reads using GreenGenes and RDP databases. Using OTU abundance and the corresponding taxonomic classifications, feature abundance matrices were calculated at different taxonomic levels, representing OTUs and taxa abundance per sample. The "Phyloseq" (Mcmurdie and Holmes, 2012) and "Vegan" (Dixon, 2003) R package were used for the detailed downstream analysis on abundance matrix.
In the end, a total of 8,010,052 paired-end 250 bp reads were obtained, 6,428,315 of which were retained as high-quality sequences ( Table S1). On average, a total of 58,015 sequences per sample were obtained in the study, with a mean length of 441 ± 15 bp. These sequences were clustered into 15,784 OTUs using the open reference-based OTU-picking process ( Table S2). Among them, 12,069 were classified taxonomically down to the genus level (Table S2). OTU counts per sample and OTU taxonomical assignments are available in Table S2. The unclassified taxa were filtered out from the analysis because the main goal of the current study was to identify specific taxa related to host susceptibility to strongyle infection.
The α-diversity indexes (observed species richness, Chao, 1984 andShannon, 1997) were calculated using the "Phyloseq" R package (Mcmurdie and Holmes, 2012). Shannon's diversity index is a composite measure of richness (number of OTUs present) and evenness (relative abundance of OTUs). The nonparametric Kruskal-Wallis test was used to compare α-diversity indexes between groups and between the different time points. If the Kruskal-Wallis test was significant between the different time points, the post-hoc Dunn test analysis, which was performed with the dunnTest function in the "FSA" package, was applied to determine which levels of the independent variable differ from each other level.
Relative abundance normalization was applied, which divides raw counts from a particular sample by the total number of reads in each sample.
To estimate β-diversity, un-weighted and weighted UniFrac distances were calculated from the OTU abundance tables, and used in principal coordinates analysis (PCoA), correspondence analysis (CA), and non-parametric multidimensional scaling (NMDS) with the "Phyloseq" R package. The Permutational Multivariate Analysis of Variance (PERMANOVA) on unweighted and weighted UniFrac distance matrices were applied through the adonis2 function with 9,999 permutations from "Vegan" R package to test for groups effect. Moreover, Bray-Curtis similarity coefficients were calculated from genera count matrix and plotted in a NMDS and CA graph to show the similarity among samples using the "Vegan" R package. Then after, the PERMANOVA test was also applied to test for groups' and time points' effect.
In addition to multivariate analysis, we used the analysis of similarities (ANOSIM) to test for intragroup dispersion. ANOSIM is a permutation-based test where the null hypothesis states that within-group distances are not significantly smaller than between-group distances. The test statistic (R) can range from 1 to −1, with a value of 1 indicating that all samples within groups are more similar to each other than to any other samples from different groups. R is ≈0 when the null hypothesis is true, that distances within and between groups are the same on average. Because multiple comparison corrections for ANOSIM were not available, the number of permutations being used on those calculations was increased to 9,999.
The Wilcoxon rank-sum test with Benjamini-Hochberg multiple test correction was used to determine the differentially abundant OTUs, phyla, families, and genera between groups. An adjusted p ≤ 0.25 was considered significant. This threshold was employed in previous microbiome studies because allows compensation for the large number of microbial taxa and multiple comparison adjustment (Lim et al., 2017).
This targeted locus study project has been deposited at DDBJ/EMBL/GenBank under the accession KBTQ01000000. The version described in this paper is the first version, KBTQ01000000. The bioproject described in this paper belongs to the BioProject PRJNA413884. The corresponding BioSamples accession numbers were SAMN07773451 to SAMN07773550.

Functional Metagenomic Predictions
The functional prediction for the 16S rRNA marker gene sequences was done using the phylogenetic investigation of communities by reconstruction of unobserved states.
(PICRUSt) (Langille et al., 2013). After excluding the unknown OTUs from the GreenGenes reference database and normalizing by 16S rRNA gene copy number, functional metagenomes for each sample were predicted from the Kyoto Encyclopedia of Genes and Genomes (KEGG) catalog and collapsed to a specified KEGG level. We used Wilcoxon rank-sum test with Benjamini-Hochberg multiple test correction to evaluate predicted pathway-level enrichments between groups. An adjusted p < 0.05 was considered as significant.

Network Inference at the Genus Level
Networks at the genus level were inferred between groups at different time points. In order to prevent the compositional effects bias typical of the classical correlations methods, we calculated the correlations among genera using the partial correlation and information theory (PCIT) approach, which identifies significant co-occurrence patterns through a datadriven methodology based on partial correlation and information theory as implemented in the PCIT algorithm (Reverter and Chan, 2008). Further details are depicted in Ramayo-Caldas et al. (2016). The genera with <0.1% mean relative abundances were excluded to acquire the results for the taxa that met the statistical conditions for correlation estimations. Nodes in the network represent the genera and edges that connect these nodes represent correlations between genera. Based on correlation coefficient and p-values for correlation, we constructed cooccurrence networks. The cutoff of p-values was 0.05. The cutoff of correlation coefficients was determined as r ≥ |0.35|. Network properties were calculated with the NetworkAnalyzer plugin in Cytoscape. We used the "iGraph" package in R to visualize the network. Strong and significant correlation between nodes (r ≥ |0.60|) were represented with larger edge width in the network.
Amplified fragments of the target genes were used and diluted 10-fold in series to produce seven standards, ranging from 2.25 × 10 7 to 2.25 × 10 13 copies per µg of DNA for bacteria and protozoa and ranging from 3.70 × 10 6 to 3.70 × 10 12 copies per µg of DNA for anaerobic fungi. Each reaction contained, in a final volume of 20 µL, 10 µL of Sybergreen Mix (Power SYBR Green PCR Master Mix, ThermoFisher, Ullkirch-Graffenstaden, France), 0.6 µM of each primer to final concentration of 300 mM, and 2 µL of standard or DNA template at 0.5 ng/µL. The primer concentration of anaerobic fungi was 200 mM and 150 mM for protozoa. The DNA template was 0.5 ng/µL. In all cases, the thermal protocol for qPCR amplification and detection included an initial step of denaturation of 10 min (95 • C), followed by 40 amplification cycles (15 s at 95 • C; 60 s at 60 • C). After each run, melting curves between 60 and 95 • C were evaluated to confirm the absence of nonspecific signals. For each sample and each gene, qPCR runs were performed in triplicate. The standard curve obtained from the reference genomic fragment was used to calculate the number of copies of bacteria, protozoa or anaerobic fungi in feces. Taking into account the molecular mass of nucleotides and fragment length, we calculated the copy number as follows: Copy number per nanogram = (NL * A * 10 −9 )/ (n * mw), where NL is the Avogadro constant (6.02 × 10 23 molecules per mol), A was the molecular weight of DNA molecules (ng), n is the length of the amplicon in base pairs or nucleotides, and mw is the molecular weight per bp or nucleotide.
Wilcoxon rank-sum tests were calculated for all possible group combinations. A p < 0.05 was considered significant.

RESULTS
The effects of natural strongyle infection on gut microbiota composition and host phenotypic variables were determined in ten resistant and ten susceptible grazing ponies over a 5-month grazing season (Figure 1). Microbiological, parasitological, hematological and biochemical measures were performed at five time points (Figure S2), hereafter referred to as days after the onset of the grazing season or grazing days (gd).

A Mixture of Abiotic and Biotic Environmental Stressors During the 43-Day Transitioning Period Induced Shifts in Immunological and Microbiological Profiles
Measured FEC demonstrated a residual egg excretion in one susceptible individual at day 0 and in two susceptible ponies after 24 gd (Figure 2A). This was due to the known imperfect efficacy of moxidectin against encysted stages of strongyle. Because we were interested in studying the effects of strongyle exposure on the gut microbiota, a pyrantel treatment was administered after 30 gd to reset luminal parasite stages to zero in every pony. This treatment had a short-lived effect and resulted in zero FECs at 43 gd (Figure 2A). Therefore, the 0-43 gd period was considered as a transitioning period resulting in mild parasite exposure and changes in diet composition, management, and environmental conditions. Indeed, ponies shifted from an indoor lifestyle (associated with a hay-concentrate based diet) to an outdoor pasture-based diet at day 0. Additionally, a heat wave took place over the first 43 gd resulting in almost no rainfall (0.2 mm) and high temperatures peaking at 37.7 • C ( Figure S3A). This situation ultimately led to grass senescence and lowered pasture quality ( Figure S3B). Mild dehydration also occurred in ponies as supported by a constant rise in measured hematocrits from 32.33 to 39.63% during the first 43 gd (Figure 2B). To account for this, blood cell counts were corrected by hematocrit through time. Recorded hematological data showed that ponies were neither anemic nor thrombocytopenic throughout the experiment (Figure S4).
Concomitantly to these challenging conditions, increased levels of some white blood cell populations were observed, namely eosinophils (2.75-fold increase, p = 5.98 × 10 −13 , Figure 2C), and monocytes (1.64-fold increase, p = 2.94 × 10 −7 ; Figure 2D) during the first 43 gd. Notably, R ponies had significantly higher levels of these immunological cells relative to S ponies (Figures 2C,D), while S ponies presented higher levels of lymphocytes ( Figure 2E). Grass samples analysis did not provide any evidence of infective larvae until 43 gd.
16S rRNA gene sequencing was used to profile the fecal microbiota of S and R ponies through time. In spite of contrasted susceptibility to parasite infection, measures of α-diversity (Observed species, Chao1 richness, and Shannon diversity index) were not significantly different between the two groups (Wilcoxon rank-sum test, p > 0.05, Figure 3A). Furthermore, the overall community structure showed no statistically significance difference in un-weighted (presence/absence) Unifrac analysis (PERMANOVA, p > 0.05) or abundance-weighted analysis (PERMANOVA, p > 0.05) during the 43 days transitioning period at the OTUs level (data not shown). Nevertheless, a large variation on the β-diversity was observed at 43 gd (Figure 3), with half of the individuals (4 S and 6 R ponies) clustering separately from all the other samples (ANOSIM test, R = 0.289, p < 0.001). This stratification matched with the age of animals. The gut microbiota species composition and spatial organization at day 43 in the oldest animals (4 S and 6 R, 6 ± 0.81 years old) were similar to that of day 0 and 24, whereas the gut microbiota of the youngest animals at 43 gd (4 R and 6 S, 4.3 ± 1.15 years old) was comparable to that of day 92 and 132 gd.
Despite the lack of significant alterations in the overall gut microbial community structure between S and R ponies during the 43 gd, we next evaluated the association between the abundance of specific gut microbial taxa and the susceptibility to parasite infection. No significant differences were identified FIGURE 3 | Estimation of the α-diversity indexes and β-diversity in susceptible and resistant animals across the experiment. (A) Estimation of the α-diversity indexes in susceptible (S) and resistant (R) across the experiment. The box color indicates the time point analyzed: (dark blue = 0 days, cyan = 24 days, pink = 43 days, green = 92 days, and yellow = 132 days); (B) Principal Coordinate analysis of weighted Unifrac distances to compare fecal communities at the level of OTUs that differ between S and R animals across the experiment. Both PC axes 1 and 2 were plotted. Together they explained 60.1% of whole variation; (C) Correspondence analyses of weighted Unifrac distances to compare fecal communities at the level of OTUs that differ between S and R animals across natural parasite infection. Both CA axes 1 and 2 were plotted; (D) Genus-level network representation between ponies across the experiment linked within a specified Jaccard distance of 0.85. Two samples were considered "connected" if the distance between them was <0.85. In all cases, the relative position of points was optimized for the visual display of network properties. The point's shape indicates the susceptibility to strongylosis (triangle: S; round: R), the node color indicates the time point analyzed: (dark blue = 0 days, cyan = 24 days, pink = 43 days, green = 92 days, and yellow = 132 days).

The Predicted Levels of Resistance Matched Observed Fecal Egg Counts Between Resistant and Susceptible Ponies and Shifts in the Gut Microbiological Composition
By the end of the transition period, e.g., after 43 gd, ponies were considered adapted to the diet and their new environmental conditions. Following the patent infection, parasite egg excretion was significantly higher (adjusted p < 0.05) in the S group after 92 (235 eggs in S and 20 eggs in R ponies on average) and 132 gd (340 eggs in S and 135 eggs in R ponies on average; Figure 2A).
Ponies' body weights ( Figure S6A) and average daily weight gains ( Figure S6B) did not show significant differences between groups after natural parasite infection, and none of them displayed clinical symptoms like lethargy or diarrhea. However, strongyle exposure induced an increase of circulating monocyte levels, which were higher (p < 0.05) in R in comparison to S ponies ( Figure 2D). But the opposite trend was found for circulating lymphocytes, which continued to be significantly enriched in the white blood cells population in S ponies during parasite infection (Figure 2E). Among the white blood cell population, R ponies also presented higher levels of eosinophils at 132 gd (p < 0.05, Figure 2C). Serum biochemical analyses revealed that the enzyme alkaline phosphatase continued to be significantly higher (p < 0.05) in R than S during the patent parasite infection, whereas the levels of albumin, cholesterol, globins, glucose, total proteins and urea remained similar between both groups ( Figure S7).
Strongyle mediated alterations in microbiota diversity and structure were investigated between the two groups of ponies. As for the transition period, both groups of ponies displayed equivalent microbial species richness and α-diversity indexes (p > 0.05, Figure 3A). However, a substantial increase in observed species richness and Chao1 indexes occurred at 92 and 132 gd relative to the other time points (p < 0.05; Dunn test; Figure S8). The weighted UniFrac distance followed by PCoA on the abundance table of OTUs ( Figure 3B) showed no distinct clustering between samples from the S and the R group, which was indicative of, if at all, minor differences in microbiota composition between the two groups of ponies during parasite infection. Nonetheless, the gut microbiota structure at 92 and 132 gd differed compared to the other time points (Figure 3B; PERMANOVA, p < 0.05), reflecting the joint effect of strongyle infection, lower temperatures and increased pasture quality and abundance ( Figure S3). Similarly, the correspondence analysis based on weighted Unifrac distance ( Figure 3C) and the Jaccard network ( Figure 3D) analyses suggested that the overall gut microbiota composition was largely similar between S and R ponies at each time point but differed from the 92 gd onward. In the same line, we did not detect significant differences between the two groups using Bray-Curtis similarity coefficients at the genus level. The stress values of the NMDS plots were high (0.207), and the cumulative inertia of the two first axes of the CA was low (13.22%) indicating that the changes were very small ( Figure S9).
Even if the overall α-and β-diversity did not differ significantly between groups following the infection, changes in relative abundance of certain genera concomitantly arose with strongyle egg excretion at 92 gd. Genera of the unclassified Ruminococcaceae and Eubacteriaceae families as well as very low abundant genera (e.g., Paludibacter, Campylobacter, Bacillus, Pseudomonas, Clostridium III, Acetivibrio) increased in S relative to R ponies (adjusted p ≤ 0.25; Figure 4). Moreover, the relative abundance of Acetivibrio, and Clostridium III highly correlated with FEC (Pearson correlation coefficient ρ > 0.60). This genera enrichment occurred simultaneously with a reduction of dominant genera such as Ruminococcus, Clostridium XIVa but also members of the Lachnospiraceae family (adjusted p ≤ 0.25, Figure 4). However, none of the differences observed at 92 gd were sustained until 132 gd. The complete list of differentially expressed genera, p-values and adjusted p-values is presented in the Table S3.
The shifts in the gut bacterial composition between S and R at 92 gd were predicted to cause functional modifications, as inferred from PICRUSt (Table S4). Functional predictions showed an enrichment trend in S ponies microbiota (adjusted p < 0.10) for the mineral absorption, protein digestion and absorption, as well as of some of the pathways related to cell motility (e.g., bacterial chemotaxis, bacterial motility proteins, flagellar assembly), lipid metabolism (sphingolipid metabolism), peroxisome, and signal transduction (phosphatidylinositol signaling system) among others compared to the R group at 92 gd ( Table S4).
The gut microbiota co-occurrence networks were marginally different between S and R ponies at 92 gd ( Figure S10). The topological properties were calculated to describe the complex pattern of inter-relationships among nodes, and to distinguish differences in taxa correlations between these two groups of ponies ( Table S5). The structural properties of the S network were slightly greater than the R network, indicating more connections and closer relationships of microbial taxa in the S group. Notably, S co-occurrence network displayed higher levels of betweenness centrality (which measures the number of shortest paths going through a given node) and higher degree levels (which describes the number of neighbors) relative to R network. Nodes with the highest degree and betweenness centrality values were identified as key genera in the cooccurrence networks. Key genera in S ponies were related to the Clostridiales order (e.g., Clostridium IV, Roseburia, Nakamurella) or Spirochaetales (e.g., Treponema), whereas key genera in the R network included members of Clostridiales order (e.g., Clostridium XIV, Roseburia) and Bacteroidales (e.g., Alloprevotella).
Interestingly, natural parasite infection increased the concentration of protozoa in S ponies at 92 gd ( Figure 5B), but did not alter bacterial concentrations ( Figure 5C) and fecal pH between the pony groups within each time point (Figure 5D). The concentration of anaerobic fungi was higher in the S group compared to R group at day 92 ( Figure 5A).

DISCUSSION
While alternative strategies are needed for a more sustainable control of horse strongyle infection, the factors contributing to the over-dispersed distribution of these parasites in their hosts remain poorly characterized (Wood et al., 2012;Kornaś et al., 2015;Debeffe et al., 2016). In their preferred niche, strongyles are surrounded by gut microbiota, and reciprocal interactions between them are expected. Our study aimed to identify the consequences of parasite infection on the gut microbiota and host physiology under natural conditions and to seek a microbiological signature of strongyle infection in resistant and susceptible ponies.
For the first time, we have characterized how gut microbiological communities and host biochemical and hematological parameters varied between resistant and susceptible ponies under natural parasite infection conditions. Eosinophilia and high levels of monocytes were observed in resistant animals, while minor changes were found in the gut microbiological composition. Therefore, it seems that horses with FEC within the range of values reported herein, i.e., 0 to 800 eggs/g, do not exhibit major perturbations of their gut microbiota structure and composition. While the relationship between horse welfare and the severity of strongyle infection remains unclear, this finding may provide extra support for the current FEC cut-offs (200 or 500 eggs per gram) used to decide if treatment is needed or not (Nielsen et al., 2010).
Under our experimental setting, the first 43 days of the trial were heavily marked by both a heat wave that caused reduced pasture yield and quality, and a transition from an indoor to an outdoor management system. The combination of these factors induced a mixture of different physiological stresses for the ponies over these 43 days. Under these challenging climatic and nutritional conditions, R ponies displayed higher levels of eosinophils and monocytes in comparison to the S ponies. Eosinophilia has been reported in experimentally challenged horses (Murphy and Love, 1997). However, no larvae were recovered from pasture samples in our study suggesting mild, if any, contamination, in line with the drought conditions that are detrimental to their survival (Nielsen et al., 2007). In addition, the mild residual egg excretion observed after 24 grazing days occurred in only two susceptible ponies, which cannot explain the increased level of eosinophils in the R group. Therefore, this differential profile in immune cell populations may result from the stressful environmental conditions as previously reported (Collier et al., 2008).
The constituent phyla, families, and genera within the gut microbiota of R and S ponies during this first 43 days of the experiment were congruent with other studies performed on horses (Costa et al., 2012(Costa et al., , 2015Shepherd et al., 2012;Steelman et al., 2012;Dougal et al., 2013;Weese et al., 2015;Venable et al., 2017).
Differences in the microbiota structure and composition during the transition period were noted not between S and R individuals but rather between 4 and 6 years old individuals, suggesting that the gut microbiota of individuals at different ages who have varying effective immune response (Lyons et al., 1999) react to abiotic stressors differently. Extreme heat temperatures and consequent nutrient scarcity in the pasture might have been sufficient to cause a transient alteration in cortisol release patterns in horses (Aurich et al., 2015) and consequently shift the gut microbiota composition (reviewed by Clark and Mach, 2016). Differences in grazing behavior between individuals might also have played a role in shaping the gut microbiota. This has been observed in other animals such as sheep in which some individuals within a herd favor less N-rich swards of grass in order to reduce the rate of ingestion of parasitic larvae (Hutchings et al., 2000).
Notably, our results showed a difference in anaerobic fungal and protozoan concentrations between S and R ponies during the transition period. As these differences were found from the beginning of the experiment, we could not establish a relationship with any of the biotic or abiotic stressors occurring during this transition period. Nevertheless, these observations suggest that a potential intrinsic resistance to strongyle infection might be associated with reduced anaerobic fungal concentrations and increased protozoan concentrations. A larger scale survey linking together estimated breeding values for strongyle infection and microbiological structure could help resolve this hypothesis. In ruminants, Hsu et al. (1991) reported a reduction in ruminal fungal zoospores in individuals with higher protozoan concentrations possibly as a result of protozoal predation and competition for nutrients. Clearly, more studies are needed to fully understand the interactions between anaerobic fungi, protozoa and the susceptibility to parasite infection. However, indisputably, individuals with different susceptibility to parasite infection present different anaerobic fungal and protozoan concentrations under worm-free conditions.
The most interesting findings brought forward by this work were obtained during the natural strongyle infection from day 43 to 132 of the experiment (Figure 6A). Congruent with our initial hypothesis, the observed FEC matched the predicted resistance levels throughout the trial, hence supporting the high reproducibility of FEC (Debeffe et al., 2016;Scheuerle et al., 2016) and the feasibility to select for more resistant individuals (Kornaś et al., 2015). Eosinophil and monocyte counts remained higher in R compared to S ponies, particularly at 132 gd, whereas lymphocytes counts continued higher in S individuals. Eosinophils are the primary effectors of a Th 2 cell response during parasite infection (Weller and Spencer, 2017). Previous studies in horses infected with cyathostomins showed a Th 2 polarization during infection (Davidson et al., 2005) and an enrichment of eosinophils at the site of infection (Thamsborg et al., 1998;Collobert-Laugier et al., 2002). Interestingly, genetically resistant breeds of sheep also displayed higher levels of eosinophils at the site of infection (Terefe et al., 2009). However, more detailed immunological characterization under controlled conditions would be needed as our data are likely obscured by confounding variation in environmental factors (e.g., extreme heat temperatures and decrease of pasture quality before the parasite infection).
Despite the limited parasite infection level monitored throughout the experiment, FEC differences between S and R ponies at 92 gd were slightly reflected in the gut bacterial composition and predicted function. However, the estimated α-and β-diversity analysis suggested no differences in the overall microbiota structure between S and R individuals. These might be related to the fact that the resolution of taxonomical classification of sequences based on a limited segment of the 16S rRNA gene, such as the V3-V4 region, is relatively low (Jovel et al., 2016). One of the approaches to increase the resolution of taxonomical classification would be assemble individuals' reads into larger fragments, technically known as contigs, which are more amenable for taxonomic classifications (Jovel et al., 2016).
Strongyle natural infection in S ponies coincided with an increase in pathobionts, such as Pseudomonas, Campylobacter, and Bacillus and a decrease in commensal genera such as Clostridium XIVa, Ruminococcus, and members of the unclassified Lachnospiraceae family (Figure 6B). Firmicutes belonging to the families Ruminococcaceae (also referred as clostridial cluster IV) and Lachnospiraceae (also referred as clostridial cluster XIVa) comprise mostly of the butyrateproducing bacteria in the human gut (Geirnaert et al., 2017). Due to butyrate's anti-inflammatory properties, it could be suggested that the higher level of helminth infection in susceptible ponies altered the abundance of butyrate-producing bacteria, which therefore modulated inflammation in the gut (Li et al., 2016). Additionally, the reduced abundance of Clostridium XIV in S ponies could have had functional importance such as immune protection against the overgrowth of pathobionts like Pseudomonas and Campylobacter. Furthermore, it is notable that a significant number of predicted bacterial functional pathways in S ponies reflected immunological mechanisms, including pathogen sensing, changes in lipids, and activation of intracellular signal transduction pathways that are critical for immune system regulation and maintaining energy homeostasis (Vassart and Costagliola, 2011). Although the mechanisms of the interactions between these specific gut genera and the inferred functional attributes require further elucidation (e.g., shotgun metagenomics), our findings suggest that specific modulations of the gut microbiota might be an effective strategy for managing parasite infections in horse.
The shifts in the bacterial structure and composition at 92 gd due to the presence of parasites in the gut lumen were associated with changes in protozoan concentrations between groups. Protozoan concentrations increased in S ponies at 92 gd, which could have had an important physiological consequence, especially in host carbohydrate metabolism (Dougal et al., 2012), or have an unclear clinical significance. It could be that strongyle infection makes the large intestine a more suitable environment for commensal protozoa. It has been observed in chimpanzees that a positive correlation between intestinal parasites and commensal protozoa arose from complex interactions mediated by the host immune response and competition among taxa (McLennan et al., 2017). In horses, Güiris et al. (2010) have shown that helminths cohabit with commensal protozoa in the gastrointestinal tract in natural conditions. In line with the seemingly symbiotic relationship between H. polygyrus and bacteria from the family Lactobacillaceae in mice (Reynolds et al., 2015), it could be speculated that strongyle and protozoa may contribute to each other's success in the ecological niche of the equine intestines. Variations in protozoan concentrations might also reflect host adaptations for nutrient accessibility and acquisition during infection. In fact, protozoa have an important role in pectin degradation although have a limited role in cellulose digestion (Julliand et al., 1999). The key contribution of Treponema, a non-pathogenic carbohydrate metabolizing bacteria (Han et al., 2011), to the co-occurrence network of S ponies at 92 gd also suggests that compensatory mechanisms were induced to degrade fiber and supply the host with micro and macronutrients needed to face the infection. It is well-known that micronutrient deficiencies as well as numerous other nondietary mechanisms, such as inflammation, disruption of the gut barrier integrity and even viral, fungal or parasite overgrowth might dictate the microbial-microbial as well as microbialenvironmental interactions within the gut (Mach and Clark, 2017).
In conclusion, we have showed that predicted resistant ponies did excrete fewer eggs than their susceptible counterparts, further advocating for FEC measurement in the field to guide anthelminthic usage. This differential susceptibility status was also associated with divergent immunological responses and slight differences in microbiological composition under worm-free conditions that could be markers of resistance status but needs further validation.
Resistance was associated with high levels of circulating eosinophils and monocytes and lower levels of lymphocytes; however, the modifications of gut microbiota composition were modest. During parasite infection, susceptible ponies presented a decrease in known butyrate producing bacteria, such as Ruminococcus, Clostridium XIVa and members of the Lachnospiraceae family, which may have led to a disruption of mucosal homeostasis, intestinal inflammation and dysbiosis. In line with this hypothesis, an increase in pathobionts such as Pseudomonas and Campylobacter were observed in susceptible ponies as well as changes in several predicted immunological pathways. Moreover, these bacterial modifications in S ponies at day 92 occurred simultaneously with an increase in protozoan concentrations. Strongyles and these microbial taxa may hence contribute to each other's success or S ponies might favor the increase of these carbohydrate-degrading microorganisms to enhance the supply of micro and macronutrients needed to face the strongyle infection.
Our results therefore suggest that susceptibility to strongyle infection occurs in the presence of moderate gut microbiological factors that affect individual risk. This investigation should be followed up by experimental work in order to establish the causative reasons for variations in the gut microbiome observed in this study, as well as SCFA levels assessed due to their role in gut health.

AUTHOR CONTRIBUTIONS
GS, AB, and NM designed the experiment. AC, GS, and NM drafted the main manuscript text. NM designed and carried out the bioinformatics and biostatistical analyses, prepared all the figures and provided critical feedback on content. VB performed the RT-qPCR analyses. FR was in charge of pony maintenance and care throughout the experiment and managed sampling. AM analyzed the chemical composition of the diet. JC and CK performed fecal egg counts. MR performed the blood analysis. All authors reviewed the manuscript and approved the final version.

FUNDING
The production of the data sets used in the study was funded by grants from the Fonds Éperon, the Institut Français du Cheval et de l'Equitation, and the Association du Cheval Arabe. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

ACKNOWLEDGMENTS
We are grateful to Jean Marie Yvon, Yvan Gaude, Thierry Blard, Thierry Gascogne, Philippe Barriere, Francois Stieau, and Adelaide Touchard for the animal sampling and management, and to Marine Guego, Marine Beinat, and Julie Rivière for participating to the sample collection during the project. We also thank Diane Esquerré for preparing the libraries and performed the MiSeq sequencing. Lastly, we are grateful to the INRA MIGALE bioinformatics platform (http://migale.jouy.inra.fr) for providing computational resources.

SUPPLEMENTARY MATERIAL
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphys. 2018.00272/full#supplementary-material (D) Heatmap of the kinship coefficient matrix, which assess the genetic resemblance between ponies. Each entry in the matrix is the kinship coefficient between two subjects. Animals are arranged in the order of their genetic relatedness. Genetically similar animals are near each other. Note that the diagonal elements did not have values above unity, showing no consanguineous mating in the families. The treatment (S or R) is delineated below to the animal name. In the heatmap, red = high values of genetic relatedness, white = low values of genetic relatedness. Figure S2 | Overview of the data analysis in the study. Effect of parasite susceptibility on gut microbiota composition, gut-related parameters and host variables across time.
Step 1: Measurement of the gut microbiota composition between susceptible (S) and resistant (R) ponies across time. This step involved the analysis of the α-diversity and βdiversity between S and R ponies across time, as well as the analysis to assess gut bacterial genera whose relative abundances changed between groups across time, the prediction of the corresponding KEGG pathways, and the inference of the cooccurrence network.
Step 2: Measurement of the fecal parasite eggs counts and the gut related parameters between S and R ponies across time. The gut parameters included the anaerobic fungal, bacterial and protozoan concentrations, as well as the pH in the feces.
Step 3: Measurement of the host variables between S and R ponies across time. The host parameters included the body weight, the daily average gain, and the hematological and biochemical parameters in blood. (B) Chemical composition (crude protein (CP), neutral and acid detergent fiber (NDF and ADF), crude fiber (CF), crude ash, and acid detergent lignin (ADL) of hay at day 0 and pasture at day 24, 43, 92, and 132. From day 1 to day 43, the lack of rainfall resulted in significant soil moisture deficits and reduced growth rates of pasture. This period coincided with the late-flowering stage of the pasture, when stems and leaves are being depleted of nutrients and herbage maturation and lignification increases. After day 43, the environmental conditions were eminently favorable to start a second pasture cycle, with a quick herbage growth, high protein content and lower fiber content. We observed that herbage protein peaked at day 92, after the intense fall rains and the increased senescence of green material. Figure S4 | Hemogram data in susceptible and resistant animals across time. The evaluation of the hemogram involved the determination of the hematocrit, total white blood cell counts, total erythrocyte count, erythrocyte indices and platelet counts. The determination of the total white blood cells (A) was performed between S (purple boxes) and R (green boxes) animals across time. The quantification of different type of leukocytes: neutrophils (B) and basophils (C) were described between S and R animals across time. The values corresponding to red blood cell distribution width (D), microcytic (E) and macrocytic (F) platelets, as well as microcytic red blood cells (G), and macrocytic red blood cells (H) were plotted. In all cases, boxes show median and interquartile range, and whiskers indicate 5th to 95th percentile. * p < 0.05 for comparison between S and R ponies in each time point.   , globins (C), glucose (D), alkaline phosphatase (E), ratio albumin/goblins (F), total proteins (G) and urea (H) in susceptible (S, purple boxes) and resistant animals (R green boxes) were plotted. Boxes show median and interquartile range, and whiskers indicate 5th to 95th percentile. * p < 0.05 for comparison between S and R ponies in each time point. Figure S8 | Estimation of α-diversity indexes (observed species richness, Chao1 and Shannon) from the OTU table between time points. Estimation of the α-diversity indexes across the experiment. The box color indicates the time point analyzed: (dark blue = 0 days, cyan = 24 days, pink = 43 days, green = 92 days, and yellow = 132 days). A Dunn test was performed as post-hoc analysis to determine which levels of the independent variable differ from each other level. For all variables with the same letter, the difference between the means is not statistically significant (p > 0.05). If two variables have different letters, they are significantly different (p < 0.05). Figure S9 | Estimation of β-diversity indexes from genera table in susceptible and resistant animals across the experiment. (A) Non-metric multidimensional scaling (NMDS) of Bray Curtis distances to compare fecal communities at the level of genera that differ between susceptible (S) and resistant (R) animals across the experiment. Microbial communities were not significantly different among R and S treatments represented by convex hulls. Distance between objects on the plot represents relative dissimilarity. Stress >0.2 indicates that differences should be interpreted with caution; (B) Correspondence analyses of Bray-Curtis distances to compare fecal communities at the level of genera that differ between S and R animals across the experiment. Both CA axes 1 and 2 were plotted. Together they explained 13.22% of whole inertia. The convex hulls color indicates the treatment group: S = purple and R = green.
Figure S10 | Co-occurrence genera network at 92 days after the entry to the pasture for susceptible and resistant ponies. The correlations among genera were calculated using the PCIT method, which identifies significant co-occurrence patterns. The size of the node is proportional to genera abundance. Node fill color corresponds to phylum taxonomic classification. Edges color represent positive (red) and negative (blue) connections, the edge thickness is equivalent to the correlation values. Only genera with a relative abundance >0.10 were included.
Table S1 | Summary of study samples and fecal bacterial 16S rRNA gene amplicon sequence datasets.  S3 | Differences in relative abundance of genera between susceptible (S) and resistant (R) animals across time.
Table S4 | Differences in predicted KEGG pathway abundance between susceptible (S) and resistant (R) animals across time. Predicted pathways values across time are represented as %, where normalized counts from a particular pathway is divides by the total number of counts in each sample. Table S5 | Topological properties of co-occurring networks obtained from susceptible and resistant animals at day 92 of the experiment.