Ethical statement
The local animal care and use committee reviewed and approved the study protocol (CEEA Val de Loire; reference: APAFIS#17427-2018080211451233 v7). All protocols were conducted following EEC regulation (no 2010/63/UE) governing the care and use of laboratory animals, effective in France since the 1st of January 2013.
Animal cohort selection and maintenance
Forty female Welsh ponies (average age 4.77 ± 2.05 years) were selected from an experimental herd of 107 ponies [78] regularly monitored for FEC during the grazing season since 2013. The selection was based on the FEC records database consisting of 703 observations of the whole herd at the time of selection. Every individual pony had been recorded at least three times over a minimum of two years (seven observations per pony on average across the entire herd). FEC data were log-transformed to correct for overdispersion and fitted a linear mixed model accounting for environmental fixed effects (month of sampling, year of selection, time since last treatment, age at sampling) as previously described [19, 79]. The individual was considered a random variable to account for the intrinsic pony susceptibility against cyathostomin infection. Based on these estimates, the 40 individuals with the most extreme forecast were chosen to create two groups of high- (HIGH) or low-shedders (LOW). The median FEC was 515.28 ± 564.04 for HIGH and 41.17 ± 99.99 for the LOW groups across the past six years.
The 40 Welsh ponies grazed from mid-April to the first week of October 2019 in a group across 19 hectares. During that time, they were subjected to targeted pyrantel treatment (Strongid® oral paste, Zoetis, Paris, France, 1.36 mg pyrantel base per kg of body weight) based on FEC measurements in mid-July 2019. They had a variable anthelmintic history but received no larvicidal treatment directed at immature cyathostomins stages three months before the study. Neither had antibiotic or anti-inflammatory therapy in the previous two months (Fig S7). Intestinal disorders such as colic were recorded in one individual before the start of the study. The horse was treated with the antipyretic metamizole (Dipyralgine®, Med’Vet, France, 5 mL per 100 kg of body weight) but did not receive antibiotics (Fig S7).
In September 2019, all individuals were housed and maintained in groups of three under natural light conditions in a 14 m2 pen with slatted floors, which precluded further nematode infections. Animals were fed 5 kg of hay per day and 600 g of concentrate pellets (Tellus Thivat Nutrition Animale Propriétaire, Saint Germain de Salles, France) consisting of barley (150 g/kg), oat bran (162 g/kg), wheat straw (184.7 g/kg), oats (200 g/kg), alfalfa (121.7 g/kg), sugar beet pulp (50 g/kg), molasses (30 g/kg), salt (7.3 g/kg), carbonate Ca (5.5 g/kg), and a mineral and vitamin mix (2 g/kg), on an as-fed basis. The mineral and vitamin mix contained Ca (28.5%), P (1.6%), Na (5.6%), vitamin A (500,000 IU), vitamin D3 (125,000 IU), vitamin E (1,500 IU), cobalt carbonate (42 mg/kg), cupric sulfate (500 mg/kg), calcium iodate (10 mg/kg), iron sulfate (1 g/kg), manganese sulfate (5.8 g/kg), sodium selenite (16 mg/kg), and zinc sulfate (7.5 g/kg) on an as-fed basis. In all cases, the concentrate in the stalls was offered and controlled individually by a caretaker who monitored the animals daily to ensure they were not sick or injured. Water was available ad libitum. The study provided no food additives, prebiotics, or probiotics that could affect gut microbiota composition. The study described herein started after a 3-week acclimation period that was considered sufficient to account for changes in diet composition, management, and environmental conditions.
Study design and sampling
On day 0, the HIGH and LOW groups were randomly divided into two subgroups blocked for age. The subgroups HIGH-TRT and LOW-TRT received a dose of pyrantel (Strongid® oral paste, Zoetis, Paris, France, 1.36 mg pyrantel base per kg of body weight) to eliminate the adult worms in the gut lumen while leaving developing stages untouched. Subgroups HIGH-CTL and LOW-CTL were not treated. To monitor the immediate anthelmintic treatment effect on gut microbiota, faecal material was collected 15, 18, 21, 24, 48, and 72 hours after treatment, corresponding to faecal excretion of the drug [80]. Then, every pony was subjected to longitudinal monitoring of faecal strongyle egg excretion and faecal microbiota weekly, including days 7, 15, 21, 28, 35, and 42 post-treatment (Fig S1). This established interaction patterns within the bacterial community and between the parasite and bacterial communities. On day 0, samples were collected before the anthelmintic treatment. The body weight of the animals was recorded every fifteen days.
Faecal sampling, faecal pH measurement, and coprology
Fresh faecal samples were collected from the rectum at every time point. Faecal aliquots for microbiota analysis were immediately snap-frozen in liquid nitrogen and stored at −80°C until DNA extraction. In contrast, faecal aliquots were processed on-site to measure the cyathostomin eggs’ faecal excretion level, and initiate larval culture for community composition analysis. The pH in the faeces was determined after 10% faecal suspension (wt /vol) in saline solution (0.15 M/ml NaCl solution).
FECs were measured as a proxy for patent cyathostomin infection. FEC was carried out using a modified McMaster technique [81] on 5 g of faeces diluted in 70 mL of NaCl solution with a density of 1.2 (sensitivity of 50 eggs/g). Immediately after homogenising, 0.5 mL aliquots of the solution were added to both chambers of a McMaster slide. At ten magnifications, parasite eggs were counted within each slide chamber under a light microscope. The number of eggs on the slide was then added and multiplied by 15 to obtain the eggs per gram for each sample. Upon observation of cyathostomin eggs during FEC analysis, the remaining faecal aliquots were subjected to larval culture to profile the nemabiome.
None of the horses demonstrated any signs of colic when observed over 24 h following anthelmintic administration. Faecal consistency and feed intake also remained standard for these animals following treatment.
Cyathostomin larval culture
Faecal samples were weighed and mixed with 30% vermiculite for incubation in a culture chamber (25°C and 60% humidity) to allow cyathostomin eggs to hatch and develop into L3 larvae. The faeces were stirred and moistened as needed for 12 days before the larvae were recovered using the Baermann technique. Faeces of the treated horses were cultured before treatment and on a weekly basis after the expected faecal egg reappearance period from day 28 to day 42. They were matched with that of the untreated ponies.
Microorganisms and cyathostomin larval DNA extraction from faecal samples
For the microbiota profiling, total DNA extraction from the 520 samples was performed as previously described [19]. DNA was extracted from 200 mg of faecal material using the EZNA Stool DNA Kit (Omega Bio-Tek, Norcross, Georgia, USA). According to the manufacturer’s instructions, the DNA extraction protocol was carried out (Omega Bio-Tek, Norcross, Georgia, USA). The DNA was then quantified using a Qubit and a dsDNA HS assay kit (Thermo Fisher).
Cyathostomin larval samples were incubated for one day and overnight at 57°C in lysis buffer (48 µL NaCl (1M), 6 µL Tris-HCl (pH 8.5; 1M), 120 µL EDTA (pH 8.0; 0.5M), 426 µL H20 RNase free, and 30 µL proteinase K solution). After the proteinase K digestion, 5 µL of RNAse solution was added and incubated for 1 h at 37°C. Then, 550 µL of phenol: chloroform: isoamyl alcohol (25:24:1) was added before centrifuging for 15 min at 14,000G in phase-lock tubes. The supernatant was then collected in a 1.5 ml tube, and 200 µL of chloroform was added and centrifuged again for 15 min at 14 000 G. The supernatant was again collected in a new tube, and 0.1X volume of 3M sodium acetate, 3X 100% ethanol and 2µL of pellet paint were added. Each sample was incubated at -20°C overnight. After incubation, a 30 min centrifugation at 14 000 G was performed. The supernatant was discarded, and 500 µL of alcohol was added to wash the pellet and centrifuged again (step performed twice). The sample was then air-dried (about 2 h) before being suspended in 30 µL of Tri-HCl.
V3–V4 16S rRNA and cyathostomin ITS-2 barcodes amplification
The bacterial V3-V4 hyper-variable regions of the 16S rRNA gene were amplified with two rounds of PCR as previously reported [19]. We added four negative controls during PCR cycles, two at PCR1 and two at PCR2, to control any source of contamination.
For cyathostomins, the ITS-2 region was PCR amplified using the NC1 and NC2 primers [82] as described previously [83]. A negative control sample and five distinct mock communities were added to establish putative biases in cyathostomin presence or absence [83]. In every case, the concentrations of the purified amplicons were checked using a NanoDrop 8000 spectrophotometer (Thermo Fisher Scientific, Waltham, USA), and the quality of a set of amplicons was checked using DNA 7500 chips onto a Bioanalyzer 2100 (Agilent Technologies, Santa Clara, CA, USA). The negative control samples did not yield a band on the agarose gel, and the concentration of the purified amplicons was undetectable (< 1 ng/μL).
All libraries were pooled at an equimolar concentration to generate an 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. According to Illumina's protocol, amplicon libraries were mixed with 15% PhiX control. For this study, one sequencing run was performed using MiSeq 500 cycle reagent kit v3 (2 x 250 output; Illumina, USA).
ITS-2 barcode data preprocessing and inferences on community composition
Amplicon sequencing data were quality-trimmed using cutadapt [84] (v. 1.14) (with options --pair-filter any --no-indels -m 50 –max-n 1 -q 15 -n 2), and subsequently handled with the Divisive Amplicon Denoising Algorithm (DADA2) algorithm [85] implemented in the R (v. 4.0.2) software version. Errors were learned, and the reads denoised, tolerating a single error for each of the forward and reverse reads, with further filtering of reads shorter than 200 bp and setting the band size parameter to 32 as recommended for ITS-2 amplicons before additional chimaera removal. The amplicon sequence variant (ASV) taxonomy assignment was subsequently performed using the IDTAXA algorithm [86] as implemented in the DECIPHER (v. 2.18.1) R package, enabling a 50% bootstrap cutoff and using the public nematode ITS-2 database v1.3 (https://www.nemabiome.ca/its2-database.html, last accessed, February 3rd, 2022). The ASV count table was subsequently handled with the phyloseq [87] package (v. 1.34). A total of 345 ASVs were identified, of which 106 with less than 50 occurrences were deemed contaminant and removed, and 16 were not assigned to any genus or species. Nemabiome samples showing less than 30 counts (n = 11 and the negative control) were not considered further. ASVs were subsequently aggregated at the species level using the tax_glom() function of the phyloseq package and removing unassigned ASVs. This left 14 strongylid species and three genus-wise ASVs for 52 samples. As a result of pyrantel treatment, 16 samples were missing due to low or no parasite egg excretion in the treated groups. However, five and three additional samples of the untreated and treated groups also failed (details provided in Supplementary information; Fig S8-S10). To deal with this issue and to promote the study of interactions between bacterial genera and cyathostomin, species counts of the treated groups were set to 0 while that of the control group were assigned to species average abundance recorded in the remaining samples (n = 8 for day 35; n = 5 for day 42). In the end, nine and 10 individuals were available for the HIGH-CTL and HIGH-TRT groups.
Species assignment for Enterococcus genus
Species-level assignments were made by taking a consensus of BLAST results for each Enterococcus ASV. That is, sequences were BLAST-ed against nt, excluding uncultured/environmental accessions. The species designations of all BLAST hits sharing the top score were collated. If all top-hit species designations agreed, a species assignment was made. Accessions with no species designation were ignored.
Host blood transcriptomic
To characterise the host response to cyathostomin infection, parasite removal and parasite re-emergence, blood samples were taken (Tempus Blood RNA tubes, Thermo Fisher) from a subset of six individuals within each group at day 0, 24h and 42 days post-treatment to perform mRNA sequencing (Supplementary information). Total RNAs were then isolated using the Preserved Blood RNA Purification Kit I (Norgen Biotek Corp., Ontario, Canada), according to the manufacturer’s instructions. RNA purity and concentration were determined using a NanoDrop ND-1000 spectrophotometer (Thermo Fisher), and RNA integrity was assessed using a Bioanalyzer 2100 (Agilent Technologies, Santa Clara, CA, United States). Total RNAs were subjected to poly-A tail enrichment and library preparation following the Illumina Truseq Stranded mRNA recommendations and sequenced on a NovaSeq 6000 platform. Adapter sequences were removed using cutadapt (v. 3.4) [84], and reads were filtered with Trimmomatic (v. 0.36) [88] for a minimum length of 100 bp and minimum leading and trailing base quality of 20. Retained reads were subsequently processed with Salmon (v. 1.4) [89] for pseudo-mapping against the Equus caballus v3 reference transcriptome (Ensembl v. 103), correcting for GC (--gcBias), position (--posBias) and sequence (--seqBias) biases and enabling selective alignment (--validateMappings). RNAseq counts were regressed upon the variables of interest to establish their respective contribution to expression variance using the variancePartition package (v. 1.20) [90]. RNA counts were filtered to retain genes with log(cpm) > 1 and corrected to account for gene-wise dispersion trend across samples with the voomWithDreamWeights() function [90]. The variance was subsequently partitioned across the effect of time, treatment, the respective pair-wise interactions between these three factors, the time x condition interaction and inter-individual variation using the fitExtractVarPartModel. The same framework was applied to the subset of high shedders to isolate the genes most affected by treatment, time, and their interaction effects. To define the transcriptomic signature associated with each variable, we retained the top 5% genes with a higher contribution to the variance. For the genes of interest, linear regression of gene count upon the treatment and time effects and blocking for the random individual effect was applied with the lme4 package. For specific comparison between groups (between high and low-shedders on day 0, between treated and untreated high-shedders on day 1) and to avoid relying on a single analytical framework, gene counts were handled with the DESeq2 package [91]. Analyses were applied to genes with at least ten counts across n individuals, n referring to the least considered group size, e.g. 12 (infection effect at day 0) or six (any other comparisons). Gene expression levels were compared between high and low-shedders before treatment (n = 12 individuals in each group) to identify the impact of cyathostomin infection in the host. To retain the most discriminant genes, a π score [92] was computed as:
Gene set enrichment analysis was run with the g:profiler2 package [93], retaining enrichment with an adjusted P-value below 5%.
Faecal microbiota data preprocessing
The DADA was also implemented using the DADA2 plug-in for QIIME 2 (v. 2021.2) to perform quality filtering and chimaera removal and construct a feature table consisting of reads abundance per ASV by sample [85]. Taxonomic assignments were given to ASVs by importing both SILVA 16S rRNA full-length database (release 138; https://www.arb-silva.de/documentation/release-138/) and SILVA 16S rRNA region-specific (515F/806R) database (release 138, 99% identity; Silva 138 99% OTUs from 515F/806R region of sequences) to QIIME 2 and classifying representative ASVs using the naive Bayes classifier plug-in [94]. The full-length and region-specific annotation tables were combined, avoiding redundancy. Subsequently, for each ASV with an unidentified genus, homology-based identifications were performed using blast (v. 2.9.0) on the NCBI 16S ribosomal RNA reference database (released on May 25th, 2020) containing 20 845 sequences belonging to archaea and bacteria for classification in Greedy run mode allowing for maximum three mismatches with an E-value smaller than 1e-8.
The feature table, taxonomy, and phylogenetic tree were exported from QIIME 2 to the R statistical environment (v. 4.1.0) and combined into a phyloseq object using the phyloseq R package (v. 1.36.0) [87]. The ASV counts per sample and ASV taxonomic assignments are available in Table S2. Abundance data were aggregated at genus, family, order, class, and phyla using the tax_glom() function of the phyloseq R package.
Microbiota biodiversity and richness analysis: core genera, α- and β-diversity
The core genera, cross-sectional genera found at every time point and individuals, were estimated using a detection threshold of 0.1% and a prevalence threshold of 95% in the microbiome R package (v. 1.15.3; http://microbiome.github.io).
This microbiome R package produced the measures of evenness, dominance, divergences, and abundance. The gut α-diversity indices between groups and time points were compared using a generalised linear model (glmer() function of the lme4 package (v. 1.1-27.1)), where the group, day, faecal pH, faecal DNA concentration, age, and the box were included as fixed effects and individuals as a random effect.
Bray-Curtis dissimilarity and unweighted and weighted UniFrac distances were calculated using the phyloseq R package to estimate β-diversity. For this, samples were rarefied at 16 875 reads of depth (minimum sampling depth in our data) to allow an equal depth using the rarefy_even_depth() function of the phyloseq R package. The minimal sequencing depth was sufficient for accurately profiling bacterial composition, as predicted by calculating the rarefaction curve for observed richness, Chao1, and Shannon index (which accounts for both abundance and evenness; Supplementary information). The β-diversity was visualised using the non-metric dimensional scaling (NMDS) in the vegan R package (v. 2.5.7) [95] using the metaMDS() function. The stress value was calculated to determine the dimensions for each NMDS. The partitioning of β-diversity into turnover (βTURN) and nestedness (βNEST) components was performed with the betapart R package (v. 1.5.6) [96].
The PerMANOVA test (a non-parametric method of multivariate analysis of variance based on pairwise distances) implemented in the adonis2() function from the vegan (v. 2.5.7) [95] R package allowed testing of the global association between ecological community structure and groups across time. The strata option (strata = horse) was used to account for repeated sampling of individual horses. Specifically, we tested the effects of group and time corrected by age, faecal pH, faecal DNA concentration, and the box on the variation of total dissimilarity between faecal microbiota. The significance of the effect of group and time was assessed in an F-test based on the sequential sum of squares estimated from a 10 000 permutations procedure. The significance threshold was chosen at an adjusted P < 0.05.
Pairwise comparisons of mean Bray–Curtis distances to group centroids among horse faecal samples were calculated using the permutational analysis of multivariate dispersion, permdisp() function in vegan package.
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 than any other representatives 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 unavailable, the number of permutations used on those calculations increased to 9 999.
The effect size of the host and environmental variables on the ASV level community ordination was tested using the envfit() function in the vegan R package (10,000 permutations, Benjamini-Hochberg adjusted P < 0.05). The envfit() function performs multivariate analysis of variance (MANOVA) and linear correlations for categorical and continuous variables. The effect size and significance of each covariate were determined by comparing the difference in the centroids of each group relative to the total variation. The obtained r2 gives the proportion of variability (that is, the main dimensions of the ordination) that can be attributed to the explanatory variables. Complementary, we employed the variation partitioning analysis varpart() function in the R package vegan to account for direct and indirect effects of environmental factors, where dissimilarities in faecal microbiota composition were considered as a response and divergence in host and environmental factors.
Microbial homogeneity across groups and time
The homogeneity of dispersions of microbiota composition between groups and time points was tested through Whittaker's index using the multivariate analyses of the homogeneity of group dispersion, that is, the distance of individual groups from the centroid of their group. The betadisper() function of the Vegan R package was employed. Moreover, to assess whether the within- and between-group variability in beta diversity were compared for each group separately through an ANOVA test on the group dispersions using the aov() function in R followed by the post hoc Tukey-Kramer at 0.95 and the multiple-test correction of Benjamini-Hochberg (adj P < 0.05).
The contribution of each genus to the ecosystem temporal dynamics through groups
MaAsLin 2 R package (v. 1.8.0) [97], a multivariate statistical framework utilising generalised linear models, was used to analyse the effect of group and time points jointly with associated covariates on genera abundances, allowing formulation of both fixed effects and random effects (for within-subject correlations) in a single unified framework. Age was not considered a potential confounding factor in the MaAsLin analyses. Multiple time points per subject were considered by including the subject in the model as random effects. MaAsLin transformed the raw abundance of each genus with an arcsine square root transformation before testing for associations. The FDR threshold was set at adjusted P < 0.05 for all MaAsLin tests, with boosting and quality control steps enabled.
Intraclass correlation coefficient
The intraclass correlation coefficient (ICC) estimation uses the variance components from a one-way ANOVA (among-group variance and within-group variance; ICC = variancebetween/[variancebetween + variancewithin]). Here, the ICC was calculated based on the longitudinal genera abundance, defining each group as a distinct class, using a linear mixed-effects model (lme4; v. 1.1.28 R package) with day fitted as a fixed effect and horse as a random effect to account for inter-horse variation, as applied in the function ICC() of the sjmisc (v. 2.8.9) R package. Inter- and intra-subject variation was estimated using all individuals within each experimental group. Higher the ICC value, the higher the correlation between time points. Moreover, conditional (R2c) and marginal (R2m) coefficients of determination were calculated using the function r.squaredGLMM() in the R package lme4. The model’s marginal R2m represents the variance explained by fixed effects, whereas the conditional R2c accounts for the variance explained by fixed and random effects combined.
Quantification of the prediction of response
The receiver operating characteristics (ROC) and precision-recall analyses were performed using the R package pROC (v. 1.18.0) package to quantify the infection status predictive power of the gut microbiota. Recall and precision were essential metrics for metagenomic classification [98]. The ROCit R package (v. 2.1.1) was employed to plot the ROC curve, and from this, it calculated an optimal cut-off was determined by the Youden method.
Pyrantel treatment effect on cyathostomin community: differential abundance analysis and multi-spatial Convergent cross-mapping (CCM)
Pyrantel efficacy was estimated from the Faecal Egg Count Reduction 14 days after pyrantel treatment [99] using the eggCounts package [100]. Egg Reappearance Period (ERP) in ponies’ faeces was considered the day when observed FEC reduction fell below 90% of that measured at 14 days post-treatment [99].
To investigate the pyrantel treatment effect on the cyathostomin community, we first focused on the difference in community composition observed on day 0 and day 42 after treatment, using available data from five untreated and seven treated ponies. Alpha diversity coefficients (estimated with the estimate_richness() function of the phyloseq package) were compared between the two-time points within each group using a t-test. The group and day effects were estimated using a PerMANOVA on Jaccard distance on species presence/absence and Bray-Curtis dissimilarity on relative abundances (using the ordinate() function of the phyloseq package), using the adonis2() function of the vegan package (1 000 permutations and accounting for pony effect using the strata option). To isolate species with significant variation in abundance between the two timepoints, count data were modelled using the DESeq2 package for the same 12 samples. This analysis was restricted to 15 species found in at least four individuals and with more than 50 counts to avoid convergence issues.
To establish how pyrantel was affecting species interaction, we inferred causal relationships between measured FEC, faecal pH and 43 faecal microbiota genera (with more than 20 000 counts; 93.8% of total abundance) using the multi-spatial Convergent Cross-Mapping (CCM) procedure [101]. This extension of the classical CCM [36] leverages a bootstrapping approach to infer the causative effect of one interactor on the other from their respective replicated time series [101]. As such, it is well adapted to shorter time series (less than 30-time points). For this analysis, the time-series data consisted of the replicate observations (n = 20 horses in the treated or untreated groups) gathered on day 0 and the weekly collected data after the expected egg reappearance period (days 28, 35, and 42) as parasites were not expected to re-appear before 28 days after treatment. Therefore, the considered time lag was uneven between the first and second observations a month apart, and this analysis implicitly assumed a steady state between day 0 and day 28. This assumption was considered conservative as variations observed in the untreated ponies may reduce the ability to detect relationships between bacteria and parasite species. In the treated individuals, the transition between day 0 and day 28 is the only one being informative about the system states. The best embedding value E, e.g. the number of past time lags used to predict the next value, was determined after examining the predictability obtained considering two or three historical records. The significance of the interaction was tested with 1 000 bootstraps retaining correlation with a P-value below 5%. This procedure was run on the treated and untreated groups separately.
Stability of the gut microbiota community following disturbance (pyrantel treatment)
Dynamic stability is an aggregated measure of the stability of a species interaction matrix across time, whose value distinguishes between unstable communities (stability > 1) or more stable assemblage (stability < 1) [102]. To establish how pyrantel treatment rewired the interaction network of bacterial species, we estimated the dynamic stability of the gut microbiota collected from the 20 treated and 20 untreated individuals following a previously described approach [102]. We first isolated the interacting major bacterial genera (43 taxa represented by more than 20 000 counts and amounting to 93.8% of total sequences) using the multispatial CCM package [101], varying embedding from 2 to 7 and retaining the value with maximal predictability. After 1 000 bootstraps, species pairs showing significant predictability (correlation above the 95% confidence interval upper limit) and convergence (difference between correlation estimated with the smallest and largest library size above 0.1) were deemed to have significant interaction with each other (n = 151 distinct pairs between 43 genera) [101]. Species interaction strengths between these pairs were subsequently derived at each time point using the multivariate S-map method [103]. This method produces weighted multivariate linear regression models that predict the future interaction values between one species and its interactant using past observed interactions in the time series. S-maps hence model a given species abundance at time t + 1 as the sum of an intercept and partial derivatives (interaction strength) of that species abundance with its interactant abundances at the time t [101]. In this respect, the number of partial derivatives must equal the best embedding E determined from CCM for that species. As a result, time lags of that species were used as a replacement when the number of interactant species was lower than E, as described previously [102]. The community dynamics can thus be summarised in an interaction strength matrix (partial derivative matrix) whose dominant eigenvalue (absolute value of its real part) corresponds to the dynamic stability [102]. S-maps were estimated with the rEDM (v. 0.7.5) R package [104], and eigenvalues were estimated as described previously [102]. This workflow was applied separately to the gut microbiota of the treated (n = 20) and untreated ponies (n = 20), yielding the dynamic stability coefficients associated with each of the 40 bacterial communities. As empirical dynamic modelling relies on evenly spaced time series, this analysis was restricted to the weekly observations, ignoring the six-time points collected during the first three days post-treatment.
DATA AVAILABILITY
The gut metagenome 16S rRNA targeted locus data are available in the DDBJ/EMBL/GenBank under the BioProject PRJNA438436 and accessions from SRR18489300 to SRR18489819. The accession numbers of the BioSamples included here are SAMN26941952 to SAMN26942471. The nemabiome data used in this paper spanned BioSamples SAMN29052217 to SAMN29052295 of the SRA BioProject PRJNA849212. RNAseq data have been released as part of the BioProject PRJNA849212.
Data sets and products generated from the raw sequence data are available at the INRAE institutional data repository powered by Dataverse with DOI: https://doi.org/10.57745/AS7NTU. The validation set data are available under BioProject PRJEB39250, PRJNA433202, ERP118335, and at http://doi.org/10.17632/g7chkjrp8f.1.
All other data is available in the Supplementary Data.
CODE AVAILABILITY
The R-scripts used are available at the INRAE institutional data repository powered by Dataverse with DOI with public access: https://doi.org/10.57745/AS7NTU.