Comparison of two molecular barcodes for the study of equine strongylid communities with amplicon sequencing

Mock community design and DNA extraction

To compare both barcodes and to quantify biases in cyathostomin community prediction, DNA mixtures were made from single worms of 11 species at most. These mixtures are refered to as ‘mock communities’ as their compositions were known, although they do not represent typical worm communities as encountered in the field, where multiple individuals from each species would be available. This choice was grounded by (i) the limited availability of the worm material, (ii) the inability to distinguish between the role of differential worm contribution to the DNA pool and the contribution of within-species diversity to any biases in community diversity estimate after pooling worms together.

Mock communities were built from morphologically identified equine cyathostomin specimens from pooled faecal samples in Ukraine (Kuzmina, Dzeverin & Kharchenko, 2016) and Poland. For each species and community size, a single adult male was digested using proteinase K (Qiagen, Hilden, Germany) in lysis buffer, before DNA extraction using a phenol/chloroform protocol. DNA was precipitated overnight in ethanol and sodium acetate (5M) at −20 °C and washed twice with 70% ethanol. The resulting DNA pellet was resuspended in 30 µL of TE buffer (10 mM Tris-HCl, 0.1 mM EDTA, pH 8.0), and DNA was quantified using using the Qubit® double-stranded high-sensitivity assay kit (Life Technologies™, Carlsbad, CA, USA) with a minimum sensitivity of 0.1 ng/µL. Extracted DNA was stored at −20 °C.

To quantify the impact of the number of species in the community, mock communities of 11 and five species were built and subjected to amplicon sequencing using the ITS-2 rDNA and mitochondrial COI barcodes (Table 1). Within each community, species DNA was either added on an equimolar basis or at their respective concentrations to mimic differences occurring when mixing species with differential contributions (heterogeneous communities; Table 1). Of note, two Cyathostomum pateratum individuals were added in the five- and 11-species communities to assess the impact of inter-individual variation. To further measure the resolution ability of the nemabiome approach, two-species communities were made with Cyathostomum catinatum and C. pateratum. Both species were either in imbalanced ratios (three-fold difference in both directions) or equal DNA concentration.

Parasite material collection for comparison of the metabarcoding performances across sample type

Parasite material was collected from six Welsh ponies with patent strongylid infection (Table 2). Faecal matter (200 g) was recovered from the rectum and incubated with 30% vermiculite at 25 °C and 60% humidity for 12 days before third-stage larvae were recovered using a Baerman apparatus.

To estimate larval concentration, thirty 5-µL aliquots were sampled to count the larvae in each. Strongylid eggs were extracted from another 200 g of faeces. For this, faecal matter was placed onto a coarse sieve to remove large plant debris, before further filtering was made on finer sieves (150 µM and 20 µM mesh). Kaolin (Sigma K7375; Sigma-Aldrich, St. Louis, MO, USA) was then added (0.5% w/v) to the egg suspension to further pellet contaminating debris (5 min centrifugation at 2,000 rpm). The supernatant was discarded and the egg pellet was resuspended in a dense salt solution (NaCl, d = 1.18) and centrifuged slowly (1,200 rpm for 5 min), before this final egg suspension was placed on a 20 µM mesh sieve for the last wash. The eggs were then counted in 30 drops of 5 µL to estimate their concentration. Adult worms were collected from the same ponies at 18 and 21 h after a pyrantel embonate treatment (Strongid®; Zoetis, Malakoff, France; 6.6 mg/Kg body weight). Pyrantel efficacy was 87.5% (95% confidence interval [78.5–93.7%]) for this isolate as described elsewhere (Boisseau et al., 2022). DNA extraction was performed as for the mock community samples. All DNA concentrations were standardized to 5 µL/mL.

COI and ITS-2 primer design

We aimed to define a 450-bp amplicon within the 650-bp fragment of the cytochrome c oxidase subunit I (COI) locus previously described (Bredtmann et al., 2019; Louro et al., 2021). This would leave at most 50 bp overlap, thereby allowing sequencing error correction and better amplicon resolution (Edgar & Flyvbjerg, 2015). For this, 18 mitochondrial sequences of 17 strongylid species with complete mitochondrial genomes available at that time (October 9, 2020) were retrieved from GenBank using the PrimerMiner package v4.0.5 (Elbrecht & Leese, 2017b). They encompassed 11 species of Cyathostominae (AP017681: Cylicostephanus goldi, GQ888712: Cylicocyclus insigne; NC_032299: Cylicocyclus nassatus; NC_035003: Cyathostomum catinatum; NC_035004: Cylicostephanus minutus; NC_035005: Poteriostomum imparidentatum; NC_038070: Cyathostomum pateratum; NC_039643: Cylicocyclus radiatus; NC_042141: Cylicodontophorus bicoronatus; NC_042234: Coronocyclus labiatus; NC_043849: Cylicocyclus auriculatus; NC_046711: Cylicocyclus ashworthi) and six species of Strongylinae (AP017698 and Q888717: Strongylus vulgaris; NC_026729: Triodontophorus brevicauda; NC_026868: Strongylus equinus; NC_031516: Triodontophorus serratus; NC_031517: Triodontophorus nipponicus). These sequences were aligned with Muscle v.3.8.21 (Edgar, 2004). Subsequently, this alignment was used to quantify sequence heterozygosity for 450-bp sliding windows using a custom python script (File S1). The consensus sequence of the region with the highest diversity, i.e., best discriminant across species, was isolated to design primers with the Primer3 blast web-based interface (Untergasser et al., 2012). Parameters were chosen to have an amplicon product of 400–450 bp, primers of 20 bp with melting temperatures of 60 ± 1 °C. Primer sequences were subsequently degenerated to account for identified SNPs in the mitochondrial sequence alignment, yielding a 24-bp long forward (5′-RGCHAARCCNGGDYTRTTRYTDGG-3′) and 25-bp long reverse (5′-GYTCYAAHGAAATHGAHCTHCTHCG-3′) primers. For the ITS-2 barcode, we relied on previously described NC1 (5′-ACGTCTGGTTCAGGGTTGTT-3′) and NC2 (5′-TTAGTTTCTTTTCCTCCGCT-3′) primers applied on strongylid species (Gasser et al., 1993) and used in previous phylogenetic (Hung et al., 2000) and metabarcoding experiments (Poissant et al., 2021). In both cases, a random single, double or triple nucleotide was added to the 5′ primer end to promote sequence complexity and avoid signal saturation. A 28-bp Illumina overhang was added for the forward and reverse sequences respectively, for subsequent ligation with Illumina adapters.

Library preparation and sequencing

For library preparation, PCR reactions were carried out in 80 µL with 16 µL HF buffer 5X, 1.6 µL dNTPs (10 mM), 4 µL primer mix containing forward and reverse primers, 0.8 µL Phusion High-Fidelity DNA Polymerase (2 U/µL; Thermo Scientific, Waltham, MA, USA), and 2 µL of a 5 ng/µL genomic DNA solution quantified using a Qubit® double-stranded high sensitivity assay kit (Life Technologies™, Carlsbad, CA, USA) with a minimum sensitivity of 0.1 ng/µL. A nested PCR approach was considered to compensate for the diversity found across species over the COI region and the suboptimal amplification efficacy for the considered barcode. The COI region was first amplified using the following conditions (Duscher, Harl & Fuehrer, 2015): 95 °C for 3 min, then 30 cycles at 95 °C for 30 s, 50 °C for 30 s, 72 °C for 30 s then a final extension of 2 min at 72 °C. The resulting PCR products were then diluted at 1/70^th for the second round using the degenerate primers targeting the chosen 500 bp amplicon. The PCR parameters for this second round were 95 °C for 3 min, followed by a pre-amplification with five cycles of 98 °C for 15 s, 45 °C for 30 s, 72 °C for 30 s, followed by 35 cycles of 98 °C for 15 s, 55 °C for 30 s, 72 °C for 30 s then a final extension of 2 min at 72 °C. In that case, the five-cycle pre-amplification at lower temperature were applied to compensate for suboptimal priming of the degenerate primers (Kwok et al., 1994). For comparison between the two barcodes, the number of amplification cycles was kept identical to the original application on equine strongylids (Hung et al., 2000): 95 °C for 3 min for the first denaturation, then 30 cycles at 98 °C for 15 s, 60 °C for 15 s, 72 °C for 15 s, followed by a final extension of 72 °C for 2 min.

For each sample, 20 µL were examined on 1% agarose gel to check for the presence of a PCR amplification band at the expected product size (or absence thereof for negative controls). PCR products were purified with magnetic beads (0.8X; AMPure XP, Beckman Coulter, Brea, CA, USA) following the manufacturer’s recommended protocol. A homemade six-bp index was added to the reverse primer during a second PCR with 12 cycles using forward primer (-AATGATACGGCGACCACCGAGATCTACACTCTTTCCCTACACGAC-) and reverse primer (-CAAGCAGAAGACGGCATACGAGAT-index-GTGACTGGAGTTCAGACGTGT-) for single multiplexing. The amplicons were purified using 1X Ampure XP beads (Beckmann Coulter, Brea, CA, USA) and the concentration was checked using a NanoDrop 8000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA). The quality of a set of amplicons was controlled using a Fragment Analyser (Agilent Technologies, Santa Clara, CA, USA).

The final libraries had a diluted concentration of 5 to 90 ng/µL, and 150 ng of each library was pooled for combined library production. Quantification was done by qPCR using the Kapa Library Quantification Kit (Roche, Basel, Switzerland) and loaded onto the Illumina V3 500 cycles MiSeq cartridge (2 × 250 output; Illumina, San Diego, CA, USA) according to the manufacturer instructions. The quality of the run was checked internally using 15% of PhiX control, and then each pair-end sequence was assigned to its sample with the help of the previously integrated index.

Analytical pipelines

Quantitative PCR (qPCR) assay for species-specific amplification

To quantify any biases in PCR amplification before sequencing, single-species DNAs used to make the mock communities were subjected to quantitative PCR reactions with the ITS-2 and COI-specific primers. Other worms collected following pyrantel treatment in an experimental pony herd were also used for qPCR to avoid relying on a single sample. Their DNA was used for sequencing of the ITS-2 and COI barcodes and belonged to the C. ashworthi, C. catinatum, C. goldi, C. labiatus, C. leptostomum, C. longibursatus, C. nassatus and C. pateratum species. One male and one female were available, but for C. catinatum and C. longibursatus (single male collected).

The DNA was diluted at 1:250 and 1 µL of DNA was added to each reaction. qPCRs were carried out on a Biorad CFX Connect Real-Time PCR Detection System following the iQ SYBR GREEN supermix® protocol (1708882; Biorad, Marnes-la-Coquette, France). Reactions were run in triplicate for each species with 40 amplification cycles: 95 °C for 3 min for the first denaturation, then 45 cycles of 98 °C for 15 s, 60 °C for 30 s, 72 °C for 40 s, followed by a melt curve (65 °C to 95 °C).

Statistical Analyses

Statistical analyses were run with R v.4.0.2 (R Core Team, 2016). Community compositional data were imported and handled with the phyloseq package v1-34.0 (McMurdie & Holmes, 2013). Abundance data (read count for the ITS-2 rDNA region, or scaled read depth for COI) were aggregated at the species level using the taxglom() function of the phyloseq package v1-34.0 and converted to relative abundances for further analysis. These data were used to compare each barcode and pipeline predictive ability in the first comparison. After the most appropriate bioinformatics parameters were identified for each barcode, ASVs with outlier representation, i.e., less than 100 reads per base pair for the COI barcode or 40 reads for the ITS-2 rDNA barcodes, were regarded as likely contaminants and further discarded.

Species richness, alpha-diversity and beta-diversity analyses were conducted with the vegan package v2.5-7 (Oksanen et al., 2017). PerMANOVA was implemented using the adonis() function of the vegan package v2.5-7 (Oksanen et al., 2017).

To monitor the predictive ability of each pipeline and barcode, the precision (the proportion of true positives among all positives called, i.e., true positives and false positives) and recall (the proportion of true positives among all true positives, i.e., true positives and false negatives) of the derived community species composition were computed and combined into the F1-score as:

$$F1=2\times {\displaystyle \frac{recall\times precision}{recall+precision}}$$

This score supports the ability of a method to correctly identify species presence while minimising the number of false-positive predictions, i.e., a trade-off between precision and recall (Hleap et al., 2021).

Alpha diversity was estimated using the Shannon and Simpson’s indices using the estimate_richness() function of the phyloseq package (McMurdie & Holmes, 2013). The difference between the expected mock community expected and observed alpha diversity values was further considered to compare conditions and barcodes. The divergence between the inferred community species composition and the true mock community composition was estimated from a between-community distance matrix based on species presence/absence (Jaccard distance) or species relative abundances (Bray-Curtis dissimilarity) using the vegdist() function of the vegan package (Oksanen et al., 2017). Compositional differences between sample types were visualised with a non-metric multidimensional scaling (NMDS) with two dimensions using Jaccard and Bray-Curtis dissimilarity.

For each of these variables (F1-score, alpha-diversity differences and divergence) and within each barcode, estimated values were regressed upon bioinformatics pipeline parameters (mxee, truncation length and band size parameters for the ITS-2 barcode data; k-mer size, window size, mismatch penalty and mapping quality for the COI barcode) to estimate the relative contribution of these parameters and determine the most appropriate analytical pipeline for each barcode independently. The model with the lowest Akaike Information Criterion value was first selected with the stepAIC() function of the R MASS package v.7.3-55 (Venables & Ripley, 2002) to retain the most relevant parameter combination (model with the lowest AIC). Parameter values were then chosen according to their least-square mean estimate. These analyses were applied to every available data within each barcode.

These variables were subsequently regressed upon barcodes and the mock community complexity (two species vs five species and more) to estimate how the predictive performances were affected by these predictors. Correspondence between input DNA and recovered reads was estimated using Spearman’s correlation applied to the homogeneous and heterogeneous communities separately (Table 1).

To test for differences in alpha diversity across sample type (eggs, larvae or adult worms), the Shannon index was regressed upon the sample type and barcode using the lm() function.

To estimate amplification efficiencies, the threshold cycle (Ct) values were regressed upon the log10-transformed DNA concentration for each species and barcode. The PCR efficiency was subsequently derived as: $E_{i,j}={{10}^{-\frac{1}{\beta i,j_{}}}}_{}^{}$, where $E_{i,j}$ is the efficiency, and $\beta i,j_{}$stands for the regression slope of species i and barcode j.

The R script and the necessary files used for analysis are available under the INRAE data repository at https://doi.org/10.57745/MNYRFQ.

Impact of the community size on the ITS-2 and COI barcode performances

To ensure a fair comparison, the most appropriate bioinformatic processing was determined for each barcode according to their predictive performances of mock community composition (Supplemental Information and Tables S1–S3). This comparison relied on the same two mock communities of five cyathostomin species (Table 1). The combination of the minimap default values for the mismatch penalty (B = 4) and the window size (w = 10) parameters, with a k-mer size of 10 base pairs and no further filtering on mapping quality (MQ = 0) was deemed as the most appropriate pipeline for the COI barcode in this study. For ITS-2, stringent tolerance in the maximal expected number of errors and a truncation length of 200 bp were chosen for downstream analyses.

With these settings, the number of reads available for the considered mock communities ranged between 6,164 and 151,606 read pairs for the COI barcode, while these numbers ranged between 1,243 and 88,845 non-chimeric pairs for the ITS-2 rDNA barcode.

No significant difference was found in the F1-score between the ITS-2 and COI barcodes overall (F_1,19 = 0.26, P = 0.61; Fig. 1). However, higher predictability was obtained with the ITS-2 barcode for the more complex community as suggested by the significant interaction term between mock community size and barcode predictors (t₁ = 2.7, P = 0.02; Fig. 1). Regarding species detection, the COI barcode ability to recover the true species composition was suboptimal. This barcode recovered nine correct species at most in the most complex communities and systematically overlooked C. leptostomum and C. labratus. On the contrary, it also identified C. coronatus (in all of the four 11-species mock communities) and C. minutus (in one out of the four 11-species mock communities) despite these species absence. Their relative abundances remained however lower than 1% for C. coronatus and 0.02% for C. minutus and were associated with low mapping quality (Phred quality score <6). The ITS-2 rDNA barcode performed better with ten species detected overall, but C. goldi was systematically overlooked in the 11-species mock communities.

The ITS-2 barcode gave better representations of the most complex cyathostomin community composition (with five or 11 species; Fig. 1B). This superiority was observed for species relative abundances (reduction of 0.34 ± 0.14 in Bray-Curtis dissimilarity relative to COI, P = 0.03). Still, it was milder for species presence/absence (decrease of 0.32 ± 0.16 in Jaccard distance relative to COI, P = 0.07). The ITS-2 barcode also gave the closest estimates of the expected alpha diversity (F_1,12 = 16.6, P = 1.5 × 10⁻⁴ and F_1,12 = 15.9, P = 1.7 × 10⁻⁴ for Simpson and Shannon indices; Fig. 1C). These differences vanished, however, when considering the mock communities composed of two Cyathostomum species (Fig. 1, P > 0.2 for the six parameters considered).

Last, the fraction of unassigned reads decreased with the mock community size (Spearman’s correlation coefficient ρ = −0.64, P = 0.02, n = 13) for the ITS-2 rDNA barcode (from 18.6% to 2.51% for the two- and 11-species communities, Fig. 1D and Table S2). Because of the mapping procedure applying to the COI region, this fraction remained negligible (2 × 10⁻⁵% for the most complex community and null otherwise; Fig. 1D). In summary, none of the two barcodes offered a perfect fit for the expected community composition, but the ITS-2 rDNA barcode was more accurate for the compositional description of the communities.

PCR amplification bias for the COI and ITS-2 rDNA barcodes

The imperfect match between the mock community and the prediction from the metabarcoding approach might relate to biases in the first PCR amplification. To test this hypothesis, qPCRs were performed on each species DNA from the same single individual used for library preparation (Fig. 2, and Table S3). The average amplification efficiency was 67.7 ± 0.24% for the COI barcode. It was above 90% for the two Cyathostomum species (Table S3) but it fell below 70% for five species. Among these, C. calicatus and C. longibursatus showed the lowest values (39% and 13% respectively, Fig. 2 and Table S3).

Omitting the outlier values found for the Cylicostephanus members (C. goldi undetected and too high efficiency for C. calicatus; Fig. 2 and Table S3), the ITS-2 rDNA barcode yielded more consistent and higher amplification efficiency on average (92.2 ± 0.03%; t₁₀ = 3.42, P = 0.006) than COI.

To confirm that the results were not specific to the considered set of worms, independent qPCR were run for male and female worms (Fig. 2 and Supplemental 3). In that case, the ITS-2 rDNA barcode showed efficiencies above 93% for every species and sex. However, the efficiency for the COI barcode dropped to 38 ± 28% on average (Fig. 2 and Supplemental 3).

The considered life stage has limited effect on the inference of community diversity

Differential worm fecundity or larval development may also contribute to distort the correlation between relative abundances measured from different sample types for some species. To test this hypothesis, nemabiome metabarcoding data were generated from six horses for three life stages of strongylids (eggs, infective larvae and adult worms collected after pyrantel treatment). The total number of species present in this population remains unknown to date. The COI barcode retrieved 12 species, of which C. radiatus was not found with the ITS-2 barcode. This latter barcode allowed the retrieval of 13 species and four additional amplicon sequence variants that were assigned at the genus level only (Fig. 3). ITS-2 data also recovered C. leptostomum and Craterostomum acuticaudatum (Strongylinae) that were not found with the COI marker (Fig. 3). The fraction of unassigned reads was 2.01% on average (range between 0% and 4.28%) and did not differ across sample types (1.6% for larval samples, 1.95% for egg samples and 2.46% for adult samples).

The community structure was generally in good agreement across sample types for both the ITS-2 and COI barcodes, although two larval samples exhibited outlying behaviours with the ITS-2 barcode (Fig. 3A). As a result, the samples from these two ponies were not considered further for analyses of the diversity with the two barcodes. The Shannon index showed no significant variation across the considered life stages for both barcodes (P = 0.91, F_2,20 = 0.1) as the observed differences fell below the resolutive power of this experiment (difference of 1.8 detectable with a significance level of 5% and power of 80%). Shannon index estimates obtained with the ITS-2 barcode were higher (difference of 0.50 ± 0.19 relative to COI, Student’s t₁ = 2.56, P = 0.019).

In agreement with this observation, PerMANOVA showed that the sample type was not a significant driver of the beta-diversity, explaining 10.6% (P = 0.53, F_2,11 = 0.53) and 6.42% (P = 0.98, F_2,11 = 0.31) of the variance in species relative abundances for the COI and ITS-2 barcode respectively (Figs. 3B and 3C). The same applied when considering the Jaccard distance on species presence/absence, whereby the sample type explained 12.3% (P = 0.87, F_2,11 = 0.63) and 8.92% (P = 0.99, F_2,11 = 0.44) of the variance for the COI and ITS-2 barcode respectively (Fig. S1).

The consistency of species relative abundance between sample types varies across cyathostomin genera

Differential larval development and fecundity may affect the correlations between different sample types. Relying on the ITS-2 based data, i.e., the most reliable barcode in our setting, correlations between species relative abundances from the three samples types were estimated for each genus. The highest consistency between inferred relative abundances was found between the egg and the larval samples (Pearson’s r = 0.96, P < 10⁻⁴, n = 64).

On the contrary, the values obtained using adult samples were less well correlated to the two others (Pearson’s r = 0.77 and 0.75 between adult worms and eggs or infective larvae respectively, P < 10⁻⁴, n = 64). These correlations were high for Cylicocyclus spp (Pearson’s r between 0.97 and 0.99, n = 20), intermediate for Cyathostomum spp and poor for both Coronocyclus spp and Cylicostephanus spp (Table 3). In the lack of major PCR amplification biases found, this may owe to differential larval development for members of the Coronocyclus spp and Cylicostephanus spp. In both cases, significant correlation was found between the eggs and larvae (Pearson’s r = 0.98 and 0.89 for Coronocyclus spp and Cylicostephanus spp respectively, Table 3) but not between the adult worms and the other two stages (P > 0.14 in all cases, Table 3). While this may result from a lack of power for Coronocyclus spp (n = 8), the number of observations was similar between Cylicocyclus spp and Cylicostephanus spp.