Amplicon-based next-generation sequencing of eukaryotic nuclear ribosomal genes (metabarcoding) for the detection of single-celled parasites in human faecal samples.

Comprehensive detection and differentiation of intestinal protists mostly rely on DNA-based methods. Here, we evaluated next-generation sequencing of eukaryotic nuclear ribosomal genes (metabarcoding) for the detection and differentiation of intestinal eukaryotic protists in the stool of healthy Tunisian individuals. Thirty-six faecal DNA samples previously evaluated by microscopy and ameboid species-specific PCRs were tested. The hypervariable regions V3-V4 and V3-V5 of the 18S rRNA gene were amplified using three universal eukaryotic primer sets and sequenced using Illumina®MiSeq sequencing. In addition, real-time PCR assays were used to detect Dientamoeba fragilis, Giardia duodenalis, and Cryptosporidium spp. The metabarcoding assay detected Blastocystis (subtypes 1, 2, and 3) and archamoebid species and subtypes (Entamoeba dispar, Entamoeba hartmanni, Entamoeba coli RL1 and RL2, Endolimax nana, Iodamoeba bütschlii RL1) in 27 (75%) and 22 (61%) of the 36 stool samples, respectively. Meanwhile, the assay had limited sensitivity for flagellates as evidenced by the fact that no Giardia-specific reads were found in any of the five Giardia-positive samples included, and Dientamoeba-specific reads were observed only in 3/13 D. fragilis-positive samples. None of the samples were positive for Cryptosporidium by any of the methods. In conclusion, a large variety of intestinal eukaryotic protists were detected and differentiated at species and subtype level; however, limited sensitivity for common flagellates was observed.

Introduction

Until very recently, intestinal parasites as a whole were considered organisms consistently taking a toll on human health. Most studies on pathogenic intestinal parasitic protists have focussed mainly on Cryptosporidium spp., Giardia (G.) duodenalis, and/or Entamoeba (E.) histolytica. Meanwhile, some intestinal protists appear to be commensals in the sense that they are not causes of intestinal symptoms, or only very rarely so, and may be found even with a high prevalence in healthy populations (Krogsgaard et al., 2018; Lass et al., 2017; Scanlan et al., 2014). Moreover, the presence of commensal protists may depend on or possibly even drive the composition of the prokaryotic part of the intestinal microbiota (Andersen et al., 2015; Krogsgaard et al., 2018; Morton et al., 2015). Nevertheless, it remains a fact that large proportions of populations in especially third-world countries (Barbosa et al., 2018; Faria et al., 2017; Higa et al., 2017; Lee et al., 2000) are carriers of species of e.g., Entamoeba, Endolimax, and Iodamoeba, the genetic diversity of which is extensive (Jacob et al., 2016; Poulsen and Stensvold, 2016; Stensvold, Lebbad, & Clark, 2012; Stensvold et al., 2011b), which is one of the reasons why the public health significance of infection/colonisation by these genera remains obscure and requires further study.

Several conventional laboratory methods are available for the detection and differentiation of intestinal eukaryotic protists (Siala et al., 2015). Although considered the gold standard for many decades, microscopic examination is generally insensitive compared with DNA-based methods and require the assistance of highly experienced microscopists (Blessmann et al., 2002; De Canale et al., 2009; Julio et al., 2015; Morio et al., 2018; Sarafraz et al., 2013; Stark et al., 2005; Verweij et al., 2007). Conventional PCRs have limitations such as PCR primer bias and inherent difficulty in detection of parasites in mixed-infection samples. As an example, sequencing of 18S ribosomal DNA (rDNA) from PCR products amplified using genus-specific primers may yield mixed chromatograms on Sanger sequencing (Scanlan et al., 2015). Genus-specific PCR therefore sometimes requires a cloning step, which makes this procedure very expensive, time consuming, and therefore inefficient. Multiplex PCRs are limited to three or four targets, cross-reaction between targets may occur especially in the case of the presence of closely related species, and the sensitivity decreases as more targets are included (Chan et al., 2016). However, the use of panels of species-specific PCRs remains an expensive, laborious, and time-consuming approach and is therefore cost-ineffective.

Mitigating the limitations of state-of-the-art molecular methods, next-generation sequencing (NGS) has made it possible rapidly to generate large amounts of data from parasitic species from a single individual, population, or environmental sample in a single run (Ekblom and Galindo, 2011; Hutchison III, 2007; Nowrousian, 2010). In addition, NGS technology offers several advantages over traditional Sanger sequencing, including high throughput, reduced cost, increased speed, high sensitivity (typically, detection of ~10 ng DNA is possible), and reduced labour to avoid cloning bias.

Most metabarcoding studies to date have focussed on the biodiversity of prokaryotic communities for practical reasons such as the availability of universal prokaryotic primers targeting 16S ribosomal DNA and detailed comprehensive databases that enable precise differentiation between genera and species (Clayton et al., 2018; Johnson et al., 2019; Lear et al., 2019). For eukaryotic organisms, the large variability in the copy numbers of 18S small subunit rRNA genes, extensive genetic diversity, lack of reference sequence data, and overall genome complexity have made NGS-based studies of the eukaryotic component of gut microbiota much less straight-forward, which explains, at least in part, the paucity of studies focusing on the eukaryotic components of gut microbiota (Amaral-Zettler et al., 2009; Escobar-Zepeda et al., 2015). Attempts to assess eukaryotic diversity using NGS-based techniques targeting 18S rDNA are in the early stages for fungi (Auchtung et al., 2018), nematodes, and marine microbes (Hino et al., 2016; Porazinska et al., 2012), and the diversity of intestinal eukaryotic protists colonizing the human gut has been elucidated only to a very limited extent.

Recently, we used our in-house metabarcoding assay for identification and differentiation of intestinal prokaryotic and eukaryotic communities, targeting hypervariable regions of nuclear small subunit rRNA genes in various types of matrices (Fredensborg et al., 2020; Krogsgaard et al., 2018; Lear et al., 2019; Stensvold et al., 2020). The objective of the present study was to demonstrate the applicability and usefulness of NGS-based detection and differentiation of intestinal eukaryotic protists in human samples from healthy Tunisian individuals.

Materials and methods

DNA samples

A total of 36 DNA samples used in a previous study and extracted from stool specimens evaluated by microscopy for intestinal eukaryotic protists (Chihi et al., 2019) were used for this study (Table 1). Briefly, faecal specimens used in this study had been preserved in S.T.A.R. buffer (Roche Diagnostics, Indianapolis, IN, USA) at −20 °C prior to DNA extraction. Total genomic DNA was extracted from the stool samples pre-treated by bead beating using the QIAamp DNA Stool MiniKit (Qiagen, Inc., Valencia, CA, USA). These DNA samples were stored at −20 °C until use. These samples were already tested by a panel of eight species-specific PCRs assays (SSA) identifying the intestinal archamoebids Entamoeba (E.) coli, Entamoeba polecki, Entamoeba hartmanni, Entamoeba histolytica, E. dispar, Endolimax (E.) nana, and Iodamoeba (I.) bütschlii RL1 and RL2 (Chihi et al., 2019) (Table 1).

Table 1

Identification and differentiation of intestinal eukaryotic protists by amplicon-based next-generation sequencing of 18 rDNA (metabarcoding) and real-time PCR for Dientamoeba fragilis and Giardia duodenalis in 36 faecal samples previously tested for intestinal protists by microscopy and Archamoeba species-specific PCRs (Chihi et al., 2019).

Sample ID⁎	Microscopy⁎	Species-specific PCRs (18S rDNA)			Metabarcoding
		Archamoebid species⁎	G. duodenalis	D. fragilis	Archamoebid species	Stramenopiles species	Flagellate species
S1	Entamoeba coli	E. coli	–	–	E. coli ST1	Blastocystis ST3, ST1	–
S2	E. coli	E. coli + E. dispar	–	–	E. coli ST2 + E. dispar	Blastocystis ST3, ST2	–
S3	E. coli	E. coli + Entamoeba hartmanni	+ (CT = 40.79)	–	E. coli ST2 + E. hartmanni	Blastocystis ST3, ST2, ST1	–
S4	E. coli	E. coli	–	+ (CT = 37.95)	E. coli ST1	Blastocystis ST1	–
S5	E. coli	E. coli	–	+ (CT = 35.13)	E. coli ST1	Blastocystis ST1	–
S6	E. hartmanni	E. hartmanni	–	–	E. hartmanni	–	–
S7	E. hartmanni	E. hartmanni	–	–	E. hartmanni	Blastocystis ST3, ST2, ST1	–
S8	E. hartmanni	E. hartmanni	–	+ (CT = 38.55)	E. hartmanni	Blastocystis ST2	–
S9	Endolimax nana	E. nana	+ (CT = 39.9)	–	E. nana	Blastocystis ST3	–
S10	E. nana	E. nana	–	+ (CT = 32.94)	E. nana	Blastocystis ST3	–
S11	E. nana	E. nana	–	+ (CT = 39.01)	E. nana	–	–
S12	E. nana	E. nana + E. hartmanni	–	–	E. nana	–	–
S13	E. nana	E. nana	–	–	E. nana	Blastocystis ST3, ST1	–
S14	E. nana	E. nana	–	+ (CT = 32.24)	E. nana	Blastocystis ST3	–
S15	Iodamoeba	Iodamoeba RL1 + E. dispar + E. nana	–	–	Iodamoeba RL1 + E. nana + E. dispar	Blastocystis ST2	–
S16	Entamoeba histolytica/dispar	E. histolytica	–	–	E. histolytica	–	–
S17	E. histolytica/dispar	E. dispar	–	–	E. dispar	Blastocystis ST2	–
S18	E. coli + E. nana	E. coli + E. hartmanni + E. nana	–	–	E. coli ST1 + E. hartmanni + E. nana	Blastocystis ST1	–
S19	Blastocystis	E. nana	–	–	E. nana	Blastocystis ST3	–
S20	Blastocystis	E. nana + E. hartmanni	–	–	E. hartmanni	Blastocystis ST1	–
S21	Blastocystis	–	–	–	–	Blastocystis ST1	–
S22	Blastocystis	–	–	–	–	Blastocystis ST3	–
S23	Blastocystis	E. hartmanni	–	–	E. hartmanni	Blastocystis ST3, ST1	–
S24	Blastocystis	–	–	–	–	Blastocystis ST3	–
S25	Giardia duodenalis	–	+ (CT = 32.6)	+ (CT = 29.45)	–	–	D. fragilis
S26	G. duodenalis	–	+ (CT = 32.03)	–	–	–	–
S27	G. duodenalis	E. coli + E. nana	+ (CT = 34.39)	+ (CT = 31.71)	E. coli ST1	Blastocystis ST3	–
S28	D. fragilis	–	–	+ (CT = 24.56)	–	–	D. fragilis
S29	D. fragilis	–	–	+ (CT = 34.55)	–	Blastocystis ST1	–
S30	D. fragilis	–	–	+ (CT = 24.82)	–	–	D. fragilis
S31	Negative	–	–	+ (CT = 31.73)	–	Blastocystis ST1	–
S32	Negative	–	–	–	–	Blastocystis ST3	–
S33	Negative	–	–	–	–	Blastocystis ST3, ST1	–
S34	Negative	–	–	+ (CT = 40.39)	–	Blastocystis ST3	–
S35	Negative	–	–	–	–	–	–
S36	Negative	–	–	–	–	Blastocystis ST1	–

CT = Real-time PCR cycle threshold value for positive [+] samples.

– = Negative.

SSA = species-specific assay.

Reported in Chihi et al., 2019 (Chihi et al., 2019).

Metabarcoding: Amplicon-based sequencing of eukaryotic ribosomal genes

Library preparation

Two PCRs were performed to construct a sequencing library using three different universal eukaryotic primers sets targeting nuclear 18S rRNA gene as previously described (Fredensborg et al., 2020; Krogsgaard et al., 2018; Ring et al., 2017; Stensvold et al., 2020). Briefly, the G3 and G6 primers target a ~ 300 bp-long DNA sequence spanning the hyper-variable regions V3-V4 of the 18S rRNA gene: G3F1: 5′-GCCAGCAGC-CGCGGTAATTC-3′ / G3R1: 5′-ACATTCTTGGCAAATGCTTTCGCAG-3′ and G6F1:

5′-TGGAGGGC-AAGTCTGGTGCC-3′ / G6R1: 5′-ACGGTATCTGATCGTCTTCG-ATCCC-3′. The G4 primers target approximately 500 bp in the hypervariable regions V3-V5 of the 18S rRNA gene: G4F3: 5′-CAGCCGCGG-TAATTCCAGCTC-3′ / G4R3:

5′-GGTGGTGCCCTTCCGTCAAT-3′.

Two microliters of genomic DNA from each of the 36 samples were used in PCR1 in a final volume of 25 μL using the RED Extract-N-Amp PCR Ready Mix (Sigma-Aldrich, St Louis, MO, USA) and 0.4 μM of each primer. The amplification profile consisted of initial denaturation for 3 min at 95 °C, followed by 20 cycles of 1 min at 95 °C, 1 min at 60 °C, 30 s at 72 °C, and a final elongation step for 4 min at 72 °C. The PCR1 products were prepared for sequencing by a second PCR (PCR2), which involved attaching an adaptor A, an index barcode i5, and a forward sequencing primer site to the 5′-end of the amplicons and an adaptor B, an index barcode i7, and a reverse sequencing primer site to the 3′-end of the amplicons. The amplification reaction of PCR2 was performed using the same PCR conditions as described above. Each PCR amplicon possessed specific barcode sequences, allowing the multiplexing of all samples in the pooled library. The PCR2 products were quantified and normalized using the Quant-IT™ dsDNA High Sensitive Assay Kit (Thermo Fisher Scientific) and mixed in equimolar amounts between samples to generate the pooled amplicon library (PAL) for further assays.

Deep sequencing using the Illumina®MiSeq platform

The PAL was purified with magnetic AMPure XP beads from Agencourt (Beckman Coulter) in a two-step process to remove undesirable DNA amplicons. First, small amplicons with a nucleotide (nt) count below 300 were removed by a 10-μL PAL-to-24-μL AMPure beads ratio, following the manufacturer's protocol. The magnetic beads sediment at the bottom of the tubes and the supernatant containing the undesirable amplicons were removed. Next, the magnetic beads were washed with 70% ethanol for 30 s and eluted in a volume of 40 μL of Tris-EDTA buffer. Secondly, amplicons larger than 1 kbp were removed by a 10-μL PAL-to-16-μL AMPure beads ratio, as previously described. The resulting AMPure bead-purified PAL was diluted with a 0.001 N NaOH concentration in order to obtain a final concentration of 11.5 pM suitable for Illumina®MiSeq sequencing. The sequencing was carried out by the Illumina®MiSeq desktop sequencer (Illumina Inc., San Diego, CA 29122, USA). The library was sequenced in a 2 × 250 nt paired-end setup with the 500-cycle MiSeq Reagent Kit V2 (Illumina Inc., San Diego, CA 29122, USA) following the manufacturer's recommendations.

Data analysis

Raw data were analysed with the software ‘BION’ (http://box.com/bion [last accessed on 25 January 2022]). The pipeline accepts raw sequence data and includes steps for de-multiplexing, primer extraction, sampling, sequence- and quality-based trimming and filtering, de-replication, clustering and chimera-checking, reference data mapping, and taxonomic mapping and formatting. Sequences passing quality control were mapped to reference data, using the Silva database v. 12.8, with a custom-built eight tier-based eukaryotic taxonomy and a similarity cut-off of 85% k-mer similarity, with confirmatory sequences extraction and querying DNA sequences of interest using the National Center for Biotechnology Information (NCBI)’s Basic Local Alignment Search Tool (BLAST) (https://blast.ncbi.nlm.nih.gov/Blast.cgi - last accessed on 25 January 2022). The software allows non-overlapping paired end reads for analysis.

TaqMan-based real-time PCR for Dientamoeba (D.) fragilis, G. duodenalis, and Cryptosporidium spp.

Since previous studies using the metabarcoding assay had revealed limited sensitivity for flagellates and sporozoa (Stensvold et al., 2021; Stensvold et al., 2020), detection of Dientamoeba (D.) fragilis, G. duodenalis, and Cryptosporidium spp. was performed by amplification of the 18S rRNA gene (98 bp, 62 bp, and 138 bp, respectively) by TaqMan Real-Time PCR. Primers and probes used for detection of D. fragilis and G. duodenalis were described by Verweij and colleagues (Verweij et al., 2004; Verweij et al., 2007). Primers and probe for Cryptosporidium spp. were described by Berg et al. (Berg et al., 2021). For each real-time PCR reaction, positive, negative, and PCR inhibition controls were included as part of the standard diagnostic protocols in place at the Laboratory of Parasitology, Department of Bacteria, Parasites & Fungi, Statens Serum Institut, Copenhagen.

PCR reactions were performed in 50-μL volumes using Hotstar Taq master mix (Qiagen), 5 mM MgCl2, and 5 μL of the DNA sample. The amplification profile consisted of 10 min at 95 °C, followed by 50 cycles of 15 s at 95 °C and 60 s at 60 °C. PCR products were analysed using Sequence Detection Software v.2.3 (Thermo Fisher Scientific). A sample was considered positive, if a target-specific exponential curve was observed with a cycle threshold (Ct) value ≤42 read at a ΔRn value of 0.1. All samples were tested in duplicates.

Results

The three universal eukaryotic primer sets of the metabarcoding assay provided a total of 1,947,435 paired-end reads reflecting the eukaryotic communities existing across all 36 faecal samples. This total represented all the 18S rDNA sequences assigned to genus/species level by the automated taxonomic mapping enabled by BION.

About half of the sequences (1,058,989 paired-end reads [54%]) aligned to intestinal eukaryotic protists, with a mean of 29,390 sequences/sample (range, 61–81,536). The other eukaryotic sequences mainly corresponded to fungi (Ascomycota and Basidiomycota), plants (Phragmoplastophyta and Ochrophyta), and vertebrates (mainly host-associated DNA).

Among the intestinal eukaryotic protist-specific sequences, a remarkably high number (923,212 [87.2%]) could be attributed to Blastocystis, while 135,436 sequences (12.7%) were identified as belonging to archamoebid species; only very few sequences (341 [<0.1%]) reflected flagellate-specific DNA. Only two samples did not contain any reads representing intestinal eukaryotic protists; however, since DNA from the host (Vertebrata), fungi (mostly Ascomycota, Mucoromycotina), and plants (Phragmoplastophyta) was detected, we did not suspect that inhibition of the PCRs had been an issue in these samples.

The diversity of intestinal eukaryotic protists as detected and differentiated by metabarcoding assay is summarized in Table 1.

Blastocystis

The metabarcoding assay identified Blastocystis-specific reads in 75% of the samples (n = 27), which was a remarkably high proportion compared with the findings obtained by microscopy; only six samples (16.6%) had been identified as Blastocystis-positive by microscopy (Table 1). Among the six samples negative for intestinal parasitic protists by microscopy, only one sample was intestinal protist-negative by the metabarcoding assay; the remaining five samples were positive for Blastocystis. By querying consensus sequences using the Nucleotide BLAST at NCBI, only a single Blastocystis subtype was seen in 20 of 27 Blastocystis-positive samples (74%); meanwhile, seven Blastocystis-positive samples contained reads from more than one subtype (26%; Table 1). ST3 was the most common subtype in our sample set (n = 16, [59%]), followed by ST1 (n = 14, [51%]) and ST2 (n = 6, [22%]).

Archamoebids

Among the 135,436 sequences belonging to archamoebids, 93% were successfully classified to species level. Only 7% of the sequences were taxonomically assigned only to family level (Entamoebidae) or genus level (Entamoeba); further inspection of the nucleotide BLAST results revealed that these clusters represented E. hartmanni. Twenty-two (61.1%) samples out of 36 were positive for at least one Archamoebid species by the metabarcoding assay. Using the SSA as the gold standard, the sensitivity of the metabarcoding assay for amoebic species was 93.5%. One sample was positive for E. histolytica by both the metabarcoding assay and the SSA (Table 1), while E. dispar was identified in three samples, again with a 100% concordance between metabarcoding assay and the SSA. A total of seven samples were positive for E. coli, and further inspection of the nucleotide BLAST results enabled the identification of E. coli subtypes, with five samples being positive for E. coli ST1 and two samples being positive for ST2. Hence, 100% agreement was observed between metabarcoding assay and the SSA with regard to detection of E. coli, including subtype designation (ST1 vs. ST2). E. hartmanni was detected in seven samples by the metabarcoding assay, all of which had also tested positive by the SSA. However, the metabarcoding assay scored sample S12 negative for E. hartmanni; this sample was found positive by the SSA. Endolimax nana was identified in eight samples by the metabarcoding assay. All these samples were also positive by the SSA. However, two samples positive for E. nana by the SSA were negative by both metabarcoding and microscopy (Table 1). With regard to I. bütschlii, metabarcoding detected RL1 in one sample, hence confirming the result of the SSA (Table 1).

Flagellates

Across the sample set, five and 13 samples were positive for G. duodenalis and D. fragilis, respectively, by real-time PCR (Table 1). However, very few sequences were assigned to flagellate taxa by the metabarcoding assay; those that were, were identified as belonging to D. fragilis. The metabarcoding assay failed to detect any G. duodenalis-specific DNA sequences, even in the three samples that were G. duodenalis-positive by microscopy. All samples positive for G. duodenalis presented a Ct value >32 (Table 1). With regards to D. fragilis, the metabarcoding detected D. fragilis-specific DNA in the three samples that presented a Ct value <30, but failed to detect D. fragilis-specific DNA in the nine samples with Ct values >30 (Table 1).

Cryptosporidium

No samples were positive for Cryptosporidium spp. by microscopy, real-time PCR, or the metabarcoding assay.

Discussion

Advances in molecular technologies have improved detection and differentiation of micro-eukaryotic taxa and enabled further insight into their roles in health and disease. Moreover, they are currently facilitating One Health approaches to surveying, managing, and controlling infectious diseases. The ideal potential of using a metabarcoding assay is the ability to use a ‘one-fits-all’ approach to detecting and differentiating intestinal eukaryotic protists with as high specificity and sensitivity as possible. 18S rDNA remains one of the most informative targets for identification and differentiation of eukaryotic communities due to its conservation among species and genera and the existence of nine hypervariable regions, which allow genetic differentiation sometimes at subtype/genotype/ribosomal lineage level (e.g., Blastocystis and archamoebids) (Clark and Diamond, 1991; Hadziavdic et al., 2014; Stensvold, Lebbad, & Verweij, 2011; Verweij and Stensvold, 2014). However, one of the challenges associated with amplicon-based NGS of eukaryotic organisms is the lack of 18S rDNA sequences in publicly available databases for certain parasite species (Hino et al., 2016; Tanaka et al., 2014). For instance, until now, only one complete 18S rDNA sequence is available in the NCBI Genbank database for E. nana, although there is evidence of extensive genetic diversity within the genus of E. nana sensu latu (Poulsen and Stensvold, 2016); this is surprising, given the likelihood of E. nana colonizing or infecting maybe hundreds of millions of people based on what can be extrapolated from data published in parasite surveys (Espinosa Aranzales et al., 2018; Faria et al., 2017; Seguí et al., 2018). To this end, genomic studies of the species E. coli, E. hartmanni, and I. bütschlii are also strikingly limited (Stensvold, 2019; Stensvold et al., 2012b).

Despite the limited number of DNA sequences available in Genbank, our metabarcoding assay detected and differentiated intestinal eukaryotic protists across our sample panel to a much larger extent than microscopy. This is in contrast with a recent study, in which a broad-specificity 18S rDNA primer set failed to detect important protozoa such as Entamoeba and Plasmodium (Kounosu et al., 2019). Optimal detection and discrimination of the diversity of protists depends on conditions such as the design of universal primers, amplicon length, competing template, DNA extraction, degree of PCR inhibition, quality of reference databases, and the variable region(s) targeted, just to mention some of the factors involved.

Cysts of E. dispar are morphologically identical to those of E. histolytica. Distinction between the pathogenic and the non-pathogenic species therefore relies primarily on DNA-based methods. The metabarcoding assay was able to differentiate these two species despite targeting the 18S gene, which offers limited discriminatory power for these two species. Indeed, using the 18S V9 region, Zahedi et al. (Zahedi et al., 2019) were not able to differentiate between E. histolytica and E. dispar.

On the contrary, extensive genetic diversity has been observed within the species of E. coli; almost 13% across the SSU rRNA gene between subtypes ST1 and ST2 (Stensvold et al., 2011b). Our assay showed 100% concordance with results obtained by the SSA (Chihi et al., 2019) with regard to the detection and differentiation of E. coli subtypes (Table 1).

The metabarcoding assay was less sensitive (87%) than specific PCR for detection of E. hartmanni. One sample (S12) was positive for E. hartmanni by conventional PCR in our previous study but came out negative by the metabarcoding assay in the present study. However, although this sample was considered E. hartmanni-positive by SSA, this finding remains unconfirmed by sequencing, owing to the low band intensity of the obtained PCR product (Chihi et al., 2019).

Metabarcoding was also slightly less sensitive than specific PCR for detection of E. nana. In fact, two samples, S20 and S27, which tested positive for E. nana by the SSA with a low-intensity PCR product, came out negative for E. nana in the metabarcoding assay. This suggests that the metabarcoding assay could have difficulty detecting species from samples with small amount of species-specific DNA due to the strong competition of universal primers of the remaining eukaryote DNA present in the sample (Hino et al., 2016; Tanaka et al., 2014). Besides, samples for which less than three sequence reads were available were considered negative. It should be noted that while the PCR product for S27 obtained by conventional PCR could be successfully sequenced using the Sanger method, the PCR product for S20 could not.

Two highly different ribosomal lineages (RLs) (~31% genetic divergence across the SSU rRNA gene) have been described for I. bütschlii sensu latu (Stensvold et al., 2012b). One sample was scored positive for I. bütschlii RL1 by the SSA method and the metabarcoding method, indicating consensus in terms of detection and differentiation (Chihi et al., 2019) (Table 1).

The sensitivity of the metabarcoding assay for detecting and differentiating DNA from flagellates was low. Only three of 13 samples positive for D. fragilis by real-time PCR could be confirmed positive by the metabarcoding assay. Perhaps unsurprisingly, the Ct values for the real-time PCR-positive samples confirmed by the metabarcoding assay were lower than those observed for the samples that could not be confirmed positive by the assay. In addition, we noticed that the three samples positive by metabarcoding for D. fragilis were negative for Blastocystis and archamoebid species. Hence, the absence of competition from abundant DNA from organisms with high-copy numbers of the 18S rRNA gene might explain why the metabarcoding assay enabled detection of D. fragilis in these three samples (Bik et al., 2013; Bradley et al., 2016; Cooper et al., 2016; Medinger et al., 2010; Tanaka et al., 2014). Blastocystis and archamoebid species have high copy numbers of rRNA genes (almost 200 copies) (Bhattacharya et al., 1989; Rezaei Riabi et al., 2018; Stensvold and Clark, 2016a), while those of at least some flagellates have low copy numbers (Zhu et al., 2005) (e.g., 60 copies of rDNA genes for G. duodenalis).One study investigated the effect of rDNA copy number on the relative abundance estimates of all species in the sample (Medinger et al., 2010). It was suggested to use either targeted blocking primers to inhibit amplification of abundant rDNA or use group-specific PCR or limit the analysis to species with similar rDNA repeat numbers (Cooper et al., 2016; Medinger et al., 2010).

When it comes to the sensitivity of the metabarcoding assay for detecting D. fragilis we recently purported the theory that the amount of DNA available in the sample from trophozoites may be very limited due to e.g., degradation (Stensvold and Clark, 2020), and cysts of D. fragilis have only very rarely been described.

Certain hypervariable regions could be suitable for most but not all species. For example, a previous study showed limited detection efficacy for nematodes and cestodes by targeting the 18SV9 hypervariable region; nevertheless, G. duodenalis was detected (Tanaka et al., 2014). The hypervariable region V9 may be more suitable for the detection of G. duodenalis. Additionally, the exact relationship between these species and our primer affinities could result in the formation of secondary structures. In fact, some species may have structures that bind with more difficulty than others and that could cause a low binding affinity for our primers, despite a 100% primer sequence match with target species.

A strongly limited ability to detect flagellate protozoa by both metagenomics and metabarcoding were observed in previous studies, especially with regard to G. duodenalis, but also for Cryptosporidium (Parfrey et al., 2014; Ramayo-Caldas et al., 2020; Stensvold et al., 2021; Wylezich et al., 2019); the reason for this remains obscure. The authors suggested that the primers were not suitable for detection G. duodenalis and that the DNA extraction method was not appropriate for detection of Cryptosporidum (Parfrey et al., 2014; Ramayo-Caldas et al., 2020; Stensvold et al., 2021; Wylezich et al., 2019). Nevertheless, the metaborcoding method used in the present study has been shown previously to enable detection of Cryptosporidium (Stensvold et al., 2021), although with a somewhat reduced sensitivity.

Several studies point towards the possibility of Blastocystis being the most widespread micro-eukaryote in the human gastrointestinal tract (Andersen and Stensvold, 2016; Beghini et al., 2017; Seyer et al., 2017; Stensvold and Clark, 2016a, Stensvold and Clark, 2020), and among the eukaryotic sequences produced in this study, those reflecting Blastocystis were by far the most common. In the present study, hypervariable regions targeted by the metabarcoding assay proved sufficiently discriminative to identify different Blastocystis subtypes in stool samples. Among the many subtypes (e.g., ST1-ST9 and ST12) so far detected in humans, the first three (subtypes 1, 2, and 3) are common all over the world, while only in Europe, ST4 can be considered common in humans (Scanlan et al., 2014; Stensvold, Alfellani, & Clark, 2012; Stensvold and Clark, 2016a, Stensvold and Clark, 2016b). Globally, ST3 appears to be the most prevalent in humans (Stensvold et al., 2012a; Stensvold and Clark, 2016a, Stensvold and Clark, 2016b). A previous study in Tunisia showed that ST3 was the predominant subtype, followed by ST1 and ST2 (Ben Abda et al., 2017). Our data were in agreement with this observation. An advantage of 18S amplicon-based NGS is the ability to explore mixed Blastocystis subtype DNA template. A previous study demonstrated that NGS-based metabarcoding is 16 times more likely to detect multiple STs in a single human host than Sanger sequencing and thereby powerful in detecting low-abundance subtypes that are missed by Sanger sequencing (Maloney et al., 2019). In the present study, the metabarcoding assay identified three combinations of mixed subtypes (ST3 + ST1, ST3 + ST2, and ST3 + ST1 + ST2) in a total of seven Blastocystis-positive samples, representing 26% of all Blastocystis-positive samples. The presence of ST3 in all mixed-Blastocystis subtype examples indicates that ST3 was the most common subtype in the sampled population.

Our metabarcoding assay based on 18S rDNA sequencing not only estimates the overall diversity and relative abundance of eukaryotic microorganisms but also informs studies on intra-generic/species diversity serving to establish taxonomic relationships between organisms and with the potential also to inform studies on host specificity and geographical variation for parasites of interest. Therefore, this metabarcoding assay could be useful for cost-effective and standardized large-scale sampling for molecular epidemiological studies of various intestinal eukaryotic protists with expectedly high applicability to studies developing a One Health approach to parasite surveillance.

We conclude that metabarcoding made it possible rapidly and cost-effectively to detect and differentiate a large variety of intestinal eukaryotic protists present in stool samples in a single run. The main challenge remaining in the present version of the metabarcoding assay is the limited sensitivity in terms of detecting DNA from flagellate species.

Funding

LOA and CRS are part of the PARADISE consortium, supported by funding from the programme under Grant Agreement No 773830: One Health European Joint Programme.

Declaration of Competing Interest

The authors declare that they have no competing interests.