Identifying group-specific primers for environmental Heterolobosa by high-throughput sequencing.

Diversity of Heterolobosea (Excavata) in environments is poorly understood despite their ecological occurrence and health-associated risk, partly because this group tends to be under-covered by most universal eukaryotic primers used for sequencing. To overcome the limits of the traditional morpho-taxonomy-based biomonitoring, we constructed a primer database listing existing and newly designed specific primer pairs that have been evaluated for Heterolobosea 18S rRNA sequencing. In silico taxonomy performance against the current SILVA SSU database allowed the selection of primer pairs that were next evaluated on reference culture amoebal strains. Two primer pairs were retained for monitoring the diversity of Heterolobosea in freshwater environments, using high-throughput sequencing. Results showed that one of the newly designed primer pairs allowed species-level identification of most heterolobosean sequences. Such primer pair could enable informative, cultivation-free assays for characterizing heterolobosean populations in various environments.

Introduction

Heterolobosea is a distinctive lineage of amoebae, amoeboflagellates and flagellates within the eukaryotic supergroup Excavata, encompassing about 150 species (Adl et al.,  2012; 2019; Pánek et al.,  2017). Almost all free‐living heteroloboseans are largely soil and freshwater sediment aerobes, although the group also includes anaerobes, obligate halophiles and many thermophiles, revealing adaptation to various extreme environments (for review, see Pánek et al.,  2017). Most heteroloboseans adopt a three‐phase asexual life cycle involving a trophozoite (amoeba) stage that can reversibly switch to either a flagellate or a cyst form. In exceptional cases, heteroloboseans can lack flagellate stage (e.g. Vahlkampfia spp.; Page, 1988) or trophozoite stage (e.g. Stephanogon spp.; Yubuki and Leander, 2008). The most studied heterolobosean members are from the Naegleria genus. As an example, N. gruberi has long been used as a cell biology model for studying flagellar apparatus development (Dingle and Fulton, 1966; Lee, 2010). N. fowleri is largely described as facultative human pathogen responsible for primary amoebic meningoencephalitis (Maciver et al.,  2020). While the disease is still very rare (a few hundred documented cases worldwide since 1968), the mortality is very high (>90%) (Gharpure et al.,  2020; Maciver et al.,  2020). N. australiensis was reported to be pathogenic for mice (De Jonckheere et al.,  1983), while Naegleria sp. was identified in brain infections of catfish species (Dykova and Lom, 2004). Strains of Tetramitus ovis, Vahlkampfia sp. and Paravahlkampfia sp. were isolated from keratitis patients, also suggesting pathogenicity (De Jonckheere and Brown, 2005; Visvesvara et al.,  2009). In addition, some heteroloboseans have been described as reservoir of various pathogenic bacteria such as Legionella pneumophila and Mycobacterium spp. (Newsome et al.,  1985; Ovrutsky et al.,  2013). Taken together, and taking into account potential effect of global warming on heteroloboseans, these studies support the importance of biomonitoring heteroloboseans in environment such as freshwaters used or exposed to various anthropogenic activities, in order to provide a better understanding of the Heterolobosea ecology, as well as the putative risks for human and animal health.

To this day, the diversity and ecology of Heterolobosea remain understudied and poorly described due to the limit of the traditional methods (Smirnov, 2003; Mrva and Garajová, 2018; Samba‐Louaka et al.,  2019). As for other amoebal lineages, the biomonitoring techniques usually require a cultivation step, followed by either a morphological identification or DNA sequencing for identification. Indeed, the compulsory cultivation step has limitations, such as low sensitivity, lack of specificity, high response time and possibly an underestimation of populations due to the presence of viable but non‐cultivable microorganisms.

High‐throughput sequencing (HTS) of PCR‐amplified loci, such as ribosomal RNA, has enabled the assessment of abundances, diversity and relatedness, of various microbial communities. While the HTS identification of microbial protists has significantly increased over the last decade (e.g. Oliverio et al.,  2020; Singer et al.,  2021), few studies are dedicated to the study of the Heterolobosea group. This is added to the fact that most universal eukaryotic primers do not provide a satisfactory coverage for this group, making it likely under‐represented in global surveys (Aucher et al.,  2020; Vaulot et al.,  2021).

Only some primer pairs, mostly targeting part of the 18S rRNA gene, claimed to be specific to Naegleria and Vahlkampfia genera or to the N. fowleri species, were reported in the literature in relation to the Heterolobosea group (Sheehan et al.,  2003; Qvarnstrom et al.,  2006; Schild et al.,  2007; Le Calvez et al.,  2012; Scheikl et al.,  2014; Régoudis and Pélandakis, 2016; Taravaud et al.,  2018). For HTS of amplified loci, the primer selection, amplicon size and the amplification efficiency across the whole Heterolobosea taxonomic range are critical to obtain an accurate representation of such communities (Sze and Schloss, 2019). Therefore, it seems necessary to compare these existing primers and to search for new primers which will make it possible to precisely characterize environmental communities of Heterolobosea, using amplicon sequencing targeting part of the 18S rRNA gene.

In this study, we tested several primer pairs targeting Heterolobosea 18S rRNA gene sequences, consisting of existing and newly designed specific primer pairs. These primers were assessed in silico against an 18S rRNA gene database and experimentally by PCR amplification from DNA isolated for amoeba culture collection. Finally, two primer pairs were retained for assessing diversity of Heterolobosea in freshwater environments using HTS, via the Illumina MiSeq technology.

Results and discussion

In silico evaluation of previously and newly described primer pairs targeting Heterolobosea.

Primer selection is critical to obtain an accurate representation of Heterolobosea communities in an environmental sample. Only four primer pairs were collected from the literature based on their reported ability to specifically amplify Heterolobosea spp. and Naegleria spp. rRNA 18S sequences, in particular (Table 1, Table S1). To increase our primer database, ten primer pairs were newly designed using PrimerDesign_M from a subset of rRNA 18S sequences from Heterolobosea class available in the current SILVA SSU 132 database (Table S2). These new primer pairs cover different variable regions of the Heterolobosea 18S rRNA gene sequence (Fig. S1). The performance of each previously and newly designed primer pairs was evaluated using TestPrime by running an in silico PCR on the SILVA database. Coverages expressed as a percentage, for taxonomic group in all the taxonomies detailed by SILVA, were reported in Table 1 The most appropriate primer pairs for our HTS application associate the highest coverage of Heterolobosea with the lowest amplification among other taxa and with an amplicon size ranging from 300 to 550 bp. From this in silico taxonomic coverage of the primer pairs against the SILVA SSU, Nae1/Nae2, Vahl730f_C/Test6R, Vahl730f_C/R1200, HF_592/HR_982, HF_1198/HR_1680 and HF_1202/HR_1680 were selected for experimental evaluations. For these six selected primer pairs, the coverage for the Eukaryota domain was low (< 3.4%) as well as the coverage for Amoebozoa and Opisthokonta classes (< 1.5% and < 5.9%, respectively). In addition to a suitable coverage of the Heterolobosea class (> 73%), the newly designed primer pairs also had the highest coverage for Naegleria genus representing 100% of the 21 Naegleria sequences in the SILVA SSU 132 database that were expected to be amplified (Table 1).

Table 1

In silico taxonomic coverage, expressed as a percentage, of the previously and newly described primer pairs against the SILVA SSU 132.

Primer pairs	Amplicon length (bp)	Euk	Amo	Dis	Ac	Tub	Vr	Exc	Het	Ng	Opi	Nu
Previously described
Nae1/Nae2	400							1.5	40.0	100.0
Nae3‐For/Nae3‐Rev	164							0.1	2.0	6.0
NaeglF192/NaeglR344	153							0.8	24.3	60.0
Vahl_560_F/ Vahl_730_R	170							2.0	70.0	66.0
Newly described
Vahl730f_C / Test6R a	330		0.1	0.3				2.0	73.0	100.0
Vahl730f_C / R1200 a	470	0.1	0.1	0.3				2.0	73.0	100.0
HF_308 / HR_648	340	0.6	3.4	3.1				1.2	35.1	80.0
HF_592 / HR_982	390	0.8	1.1	0.6		4.2		5.5	75.0	94.7	0.4
HF_413 / HR_1114	701	11.0	12.3	3.2	0.7	36.3		1.6	36.6	83.3	4.1
HF_628 / HR_1114	486	10.8	15.8	7.5	1.4	37.9		10.3	40.8	100.0	4.3
HF_628 / HR_982 a	354	12.7	12.0	3.2		1.9		5.0	79.6	100.0	17.6
HF_1245 / HR_1680 a	435	3.4	1.5	0.6				3.4	85.7	100.0	6.0
HF_1198 / HR_1680 a	482	3.4	1.5	0.6				3.4	85.7	100.0	5.9
HF_1202 / HR_1680 a	478	0.1						3.3	85.7	100.0

Aca, Acanthamoeba; Amo, Amoebozoa; Dis, Discosea; Euk, Eukaryota; Exc, Excavata; Het, Heterolobosea; Ng, Naegleria; Nu, Nucleria; Opi, Opisthokonta; Tub, Tubulinea; Vr, Vermamoeba.

The colour code indicates the coverage of primers, ranging from blue (low coverage) to red (full coverage).

Primers selected for PCR conducted on DNA from various amoeba isolates.

Because microbial eukaryotes are highly diverse, many studies concluded that no universal primer pairs can be designed to target all taxa with the same efficiency, while avoiding a resort to excessively degenerated primers. In a recent study, it was underlined that Excavata, including the Heterolobosea class, constitute a clade that is most often discriminated, in comparison with other groups (Vaulot et al.,  2021). An alternative approach to these universal primers is the application of group‐specific primers, allowing an optimal compromise between coverage and specificity of the primers. Nevertheless, our in silico analysis highlighted that the primer pairs described so far in the literature were not very efficient to target the Heterolobosea sequences present in the SILVA database since the coverages range from 24 to 70%. In addition, the Vahl_560_F/Vahl_730_R primer pair, giving the best results (70% of coverage) produced an amplicon which size (4pprox.. 170 base pairs; Le Calvez et al.,  2012), is too short and not discriminant enough for HTS applications and taxonomic identification. In fact, these primers were mostly designed to perform qPCR, requiring short amplicons. On the contrary, several primer pairs newly designed in this study cover more than 70% of the Heterolobosea sequences present in the SILVA database and were compatible in length with HTS. Moreover, they seem quite specific because they cover less than 5% of the other tested taxa.

PCR evaluation of the in silico‐selected primer pairs on isolates of amoebae

In order to confirm the in silico specificity of the six selected primer pairs, PCR experiments were conducted on DNA from various amoeba isolates, belonging to the Discosea, Tubulinea, Heterolobosea and Nucletmycea classes (Table S3). The results indicate that the Nae1/Nae2 primer pair failed to amplify DNA from all the Heterolobosea isolates like those from the Vahlkampfia and Willaertia genera (Fig. 1). Again, this primer pair was previously designed to amplify Naegleria and it succeeded to do so. Three other primer pairs were able to amplify DNA from all the Heterolobosea isolates but produced non‐specific amplification from isolates that do not belong to Heterolobosea. Thus, only two primer pairs (i) HF_1198/HR_1680 and ii) HF_1202/HR_1680) gave amplification of a single PCR product at the expected size, without amplifying any of the isolates from non‐target amoebal classes. Although we are aware that this evaluation does not guarantee that those primers will not amplify any DNA other that Heterolobosea in the environment, they seem specific enough to be tested on environmental samples. The amplicon size produced by using those two primer pairs is of 5pprox. 480 bases pairs (bp). This sequence size is likely long enough to be taxonomically resolutive, while allowing a sufficient overlap between forward and reverse sequences when using Illumina MiSeq 2x300 bp chemistry, thus ensuring to reconstruct the complete amplicon.

Fig. 1

Specificity of selected primer pairs for detection of the Heterolobosea. Agarose gel electrophoreses (2%) of PCR products using 6 different primer pairs for the environmental isolates of amoebae mentioned in Table S3. Ac: Acanthamoeba, Ec: Echinamoeba, Vr: Vermamoeba, Ng: Naegleria, Tt: Tetramitus, Vh: Vahlkampfia, Wi: Willaertia, Nu: Nuclearia, C: control PCR reaction with no template DNA. The numbers on the left correspond to the molecular weight (bp) of the BenchTop 100 bp DNA ladder (L). Results from PCR experiment conducted at 56°C.

Both HF_1198/HR_1680 and HF_1202/HR_1680 primer pairs target the same 18S rRNA region and differ by a shift of 4 bp on the forward primer. To investigate whether this commonly amplified region could provide sufficient taxonomic resolution for genus and/or species‐level identification, a distance matrix based on theoretical amplicons yielded by HF_1202/HR_1680 was computed (Fig. 2A). This analysis highlighted that most heterolobosean genera could be clearly identified and separated from each other based on amplicon sequence. Furthermore, species‐level identification could be achieved for the majority of Naegleria and Tetramitus species (Fig. 2B and C). Of note, the analysed amplicon allows to formally identify and discriminate Naegleria fowleri sequences, from other Naegleria species; a divergence of 4 bp was identified with the closest relative species, N. lovaniensis (Fig. 2B).

Fig. 2

Sequence dissimilarity based on amplicons yielded by primer pair HF_1202/HR_1680, allows for a clear identification of heterolobosean genera (A). Each genus noted on the right side is colour‐coded. B) Dendrogram based on DNA distance matrix, depicting sequence divergence among the amplicon yielded by primer pair HF_1202/HR_1680 for Naegleria (B) and Tetramitus (C) genera. Scale bar represents the percentage of nucleotidic divergence, with each dotted line representing 1% divergence. [Colour figure can be viewed at wileyonlinelibrary.com]

Overall, the set of analyses conducted in silico and in vitro suggest that two pairs of primer newly designed in this study provided an amplicon enabling a sensitive identification, down to genus and sometimes species level of heterolobosean representatives. By focusing on Naegleria genus, we were able to differentiate most species by variations of one or several nucleotidic positions within the amplicon. However, no difference could be observed among sequences from N. gruberi, N. clarki and N. italica. Interestingly, N. fowleri sequences could be clearly distinguished from their closest relative that is N. lovaniensis, by a 4‐bp difference. Thus, the primer pair newly designed in this study would be proven useful for characterizing the presence of pathogenic Naegleria in environmental samples. By focusing on Tetramitus, the second most diversified genus among Heterolobosea, we showed that amplicon sequences also allowed for species‐level identification of most species. T. thorntoni could not be distinguished from T. jugosus. Overall, our analyses suggest that fine identification, likely at the species level, of most heterolobosean can be theoretically achieved using the primer pairs HF_1198/HR_1680 and HF_1202/HR_1680.

Experimental evaluation of two newly described primer pairs on environmental samples

Based on in silico and in vitro analyses, the two newly designed primers HF_1198/HR_1680 and HF_1202/HR_1680 were further tested by using them for HTS experiments. For this, total DNA isolated from river water and sediment samples were used as templates for PCR and sequencing, using an Illumina MiSeq platform. With each primer pair, a negative template control was performed and added to the sequencing run. In total, 2′387’524 raw sequences were generated (Table 2, Table S4). After the various filtering steps, 363′887 high‐quality, eukaryotic sequences were retained for further analyses (mean of 45′486 ± 15′736 sequences per sample, excluding controls which comprised 0 or 6 sequences). Rarefaction analyses indicated that diversity contained within all samples was well covered and that the number of sequences generated for this study is largely sufficient (Fig. S2). On average, 241 ± 26 amplicons sequence variants (ASV) were identified per sample. Taxonomic identification of ASV highlighted that among those, 25 ± 11 ASVs were affiliated to the Heterolobosea group. These heterolobosean ASVs represented 26 to 66% of identified eukaryotic sequences in all environmental samples (Table 2; Fig. 3). Heterolobosean sequences thus consistently represented an important fraction of generated sequences, for all samples. By comparing both primer pairs, 29.6 ± 4.1% of total sequences were affiliated to Heterolobosea using the HF_1198/HR_1680 primer pair, compared to 41.6 ± 16.6% with the HF_1202/HR_1680 primer pair.

Table 2

Numeric and taxonomic features of sequences generated by high‐throughput sequencing, per sample.

Identifier	Origin	Sample Type	Primer pair	Eukaryotic sequences	ASV number	Heterolobosean sequences	Heterolobosean ASV	Heterolobosoan sequences (%)	Shannon H (Heterolobosea)	Overall mean genetic distance
11‐1198_S36	Vienne	Freshwater	HF_1198‐HR_1680	32 878	265	8626	14	26.24	1.60	0.62
11‐1202_S19			HF_1202‐HR_1680	25 810	216	9421	18	36.50	1.81	0.67
13‐1198_S7		Sediment	HF_1198‐HR_1680	66 964	225	18 457	40	27.56	1.99	0.59
13‐1202_S29			HF_1202‐HR_1680	64 830	236	20 171	44	31.11	2.13	0.60
22‐1198_S3	Rhône	Freshwater	HF_1198‐HR_1680	51 644	292	18 379	14	35.59	0.91	0.49
22‐1202_S25			HF_1202‐HR_1680	28 288	222	18 691	23	66.07	1.00	0.64
25‐1198_S13		Sediment	HF_1198‐HR_1680	49 474	231	14 354	24	29.01	3.04	0.32
25‐1202_S35			HF_1202‐HR_1680	43 999	244	14 257	25	32.40	3.05	0.59

Fig. 3

Taxonomic assignments of amplicons yielded by primer pairs HF_1198/HR_1680 and HF_1202/HR_1680, from 4 environmental samples. Relative abundances are shown for each sample, and sequences affiliated to the Heterolobosea clade are highlighted in red. [Colour figure can be viewed at wileyonlinelibrary.com]

The diversity within Heterolobosea, as identified using a standard taxonomic pipeline implemented in QIIME2 (vsearch, coupled with SILVA database, see methods), was further explored at the genus level (Fig. 4A). Using the conventional pipeline, three genera could be identified, namely Naegleria, Neovahlkampfia and Tetramitus. Naegleria and Tetramitus accounted for the majority of sequences in all samples, while Neovahlkampfia was only identified in one sample and only with the primer pair HF_1202/HR_1680. Also, not all sequences could be identified down to the genus level, hence the presence of unidentified Heterolobosea and/or Tetramitia in 6 out of 8 samples analysed. In most samples, Naegleria sequences represented at least 75% of all sequences. Only two samples presented a majority (around 70%) of Tetramitus sequences. The tendency is that Tetramitus was more represented in sediments that in water. Similar results regarding Heterolobosea relative abundances and diversity were obtained when using the Protist Reference database PR2 (v4.14, Fig. S3; Guillou et al.,  2013).

Fig. 4

Taxonomic assignments of heterolobosean amplicons yielded by primer pairs HF_1198/HR_1680 and HF_1202/HR_1680, from 4 environmental sample. The conventional pipeline using VSEARCH and SILVA allows genus‐level taxonomic identification (A), while the improved pipeline described in this study allow for species‐level identification of most amplicons (B). [Colour figure can be viewed at wileyonlinelibrary.com]

Because in silico analyses of amplicon taxonomic resolution suggested the possibility to identify ASV at the species level (Fig. 2), we tentatively re‐classified all sequences using a manually curated database of 90 high‐quality Heterolobosea sequences (Table S5). Identifications were performed using the BLAST algorithm, and species‐level identification was achieved when coverage and identities provided unambiguous results. Identities threshold were defined in light of previously computed distance matrix, manually inspected and confirmed for each ASV. Such an approach considerably improved the taxonomic identification of ASV in all samples (Fig. 4B). In total, eleven species were formally identified in the samples. While no difference could be observed regarding trends of genera abundances, numerous ASV could be identified down to species level, especially among Naegleria and Tetramitus genera, not all ASV could reliably be assigned to defined species, mostly due to a lack of identity or coverage. Overall, the global profiles of heterolobosean communities identified in all samples were comparable from one primer pair to another. Only minor differences emerged, such as the detection of Vahlkampfia sp. Sequences in Rhône water and of Naegleria pagei in Rhône sediments only observed using HF_1202/HR_1680. These differences are in line with the observed number of ASV and with diversity indexes, which were always slightly higher, though not significantly (Wilcoxon rank sum test with continuity correction, W = 5, P = 0.4678) for samples generated using HF_1202/HR_1680 primer pair (Table 2).

The characterization of amoebae belonging to the Heterolobosea group heavily relies on cultivation methods. The development of a primer pair, enabling the molecular analysis of heterolobosean diversity, gives access to the uncultivable fraction within this taxonomic group. Designing priming specifically covering this group was also motivated by the fact that heteroloboseans are frequently under‐covered by other, more universal primers used for describing eukaryotic diversity (see primers referenced in Vaulot et al.,  2021). Our study thus aims to provide a new tool for better characterizing Heterolobosea diversity, without cultivation steps.

The presented approach for taxonomic identification relies on manual inspection of ASV identities compared with a curated set of reference sequences. The use of ASV is essential in completing such analyses; this approach enables a nucleotide‐level comparison with reference sequences (Callahan et al.,  2017). A clear benefit of this approach is the enhancement of species‐level identification of heterolobosean, as compared to other methods. Most notably, the set of primers tested in this study, coupled with the above‐mentioned taxonomic identification method, theoretically allows the identification of the numerous heterolobosean species, as well as the unambiguous detection of Naegleria fowleri from high‐throughput amplicons sequences. The presented work thus provides a method and a new set of primers to identify heterolobosean diversity with superior precision, as compared to other known methods.

Even though relative abundances were reproducible between the two primer pairs, such values should be interpreted with caution. Indeed, 18S rRNA gene copy numbers can vary a lot within protists cells, both because of interspecific and intraspecific variations (Lavrinienko et al.,  2021). Moreover, a peculiarity of Naegleria species lies in the fact that the rRNA cistron is exclusively located on an intranuclear plasmid, found in approx. 4000 copies per cell (Clark and Cross, 1987). Early research suggested this feature might be conserved within other heterolobosean such as Tetramitus and Vahlkampfia genera, though it has not been recently investigated (Clark and Cross, 1988). Thus, large variation in 18S rRNA gene copy number across heterolobosean species could be observed, prompting caution during interpretation of relative abundance data.

Conclusion

In conclusion, we designed a new PCR primer pair dedicated to Heterolobosea identification by high‐throughput sequencing. We designed and tested several primer pairs in silico, on pure culture and on environmental samples. The experiments led to the selection of one primer pair, HF_1202/HR_1680, showing the highest specificity and allowing to detect the highest diversity among Heterolobosea. Interestingly, the amplicon sequences were discriminant enough to distinguish most Heterolobosea at the species level. This primer pair allows, for the first time, a global description of the Heterolobosea community in environmental samples. In these samples, Naegleria (in majority), Tetramitus, Vahlkampfia, Neovahlkampfia and Learamoeba were detected. This primer pair should be further tested on samples coming from various environments in order to confirm its specificity and universality towards the Heterolobosea group.

Experimental procedures

Primer design and in silico primer evaluation

Specific primer pairs target the 18S rRNA gene sequence relative to Heterolobosea, were collected from the literature. In addition, a reference 18S rRNA sequences database was created with sequences from 17 Heterolobosea, available from the current SILVA SSU 132 database (https://www.arb‐silva.de/documentation/release‐132/; Table S2). Only 18S rRNA sequences with length longer than 1400 nt were considered for alignment using Muscle (Edgar, 2004). From an optimal sequence alignment, PrimerDesign_M (Yoon and Leitner, 2015) was performed to determine sites of conservation across the Heterolobosea sequences and identify primer pairs. The Flex fragment overlap option was used, as well as the following settings: an expected read length of 280–470 nucleotides (nt), no barcode tag, a minimum primer length of 18 nt, a maximum primer length of 25 nt, a detection limit at 5%, a complexity limit at 32 and a maximum difference in Tm at 3°C. In silico taxonomic coverage of the previously and newly designed primer pairs against the SILVA SSU 132 was performed using TestPrime, allowing for a maximum of 2 mismatches (www.arb‐silva.de/search/testprime/, accessed in February 2019).

For evaluating the taxonomic resolution of amplicons yielded by selected primer pairs, a distance matrix was computed with sequence listed in Table S5. Sequences were aligned using MUSCLE (Edgar, 2004), and sequence dissimilarity was computed, using kimura parameter, in R environment and ape package (Paradis et al.,  2004; R Development Core Team, 2013). Heatmap was generated using gplot package (Warnes, 2008), while dendrogram was generated using vegan (Oksanen et al.,  2013).

DNA extraction, PCR amplification and sequencing

DNA was extracted from 13 purified isolates of amoebae deposited in our laboratory culture collection of environmental free‐living amoebae (FLA) using the Nucleospin Tissue kit according to manufacturer's instructions. Isolates were selected from the Discosea, Tubulinea, Heterolobosea and Nucletmycea classes (Table S3).

Water (n = 2) and sediment (n = 2) samples were collected between June and October 2018 from the Vienne and Rhône rivers, France. Water samples (1 l) were filtered through a 3 μm nitrocellulosic membrane (Millipore) within 24 h after sampling. Total DNA was extracted directly from membranes (i.e. without classic enrichment step) using the DNeasy PowerWater DNA and from 500 mg (wet weight) of sediment using the Dneasy PowerSoil Pro DNA Kit, with the QIAcube workstation according to Qiagen's instructions. The integrity of the extracted DNA was observed with 0.8% agarose gel electrophoresis. The amount of extracted DNA was quantified using a NanoDrop™ One spectrophotometer (Thermo Scientific™). DNA from FLA isolates and from environmental samples were diluted to a standard concentration of 5 ng μl−1, stored at −20°C until use for molecular applications.

All PCR reactions were carried in a total volume of 50 μl PCR mix with 0.5 μM of each forward and reverse primer (Table S1), 1X PCR mix buffer with 200 mM of each dNTP and 1 U of Q5 high fidelity DNA polymerase (NEB) and 25 ng of template DNA, using a vapo.protect™ Mastercycler (Eppendorf). DNA from FLA isolates and from environmental samples was amplified under the following conditions: an initial denaturation at 98°C for 1 min, followed by 30 cycles at 98°C for 10 s, 56°C for 30 s and 72°C for 30 s and a final extension performed at 72°C for 2 min. The amplified products were checked by 2% agarose gel electrophoresis. For sequencing, triplicate PCR products for each environmental sample were pooled, purified using the PCR clean‐up Kit with the QIAcube connect workstation according to the manufacturer's instructions (Qiagen) and quantified using a Qubit 2 fluorimeter (Invitrogen). Libraries were prepared using the NEBNext® Ultra II DNA library Prep Kit for Illumina (New England Biolabs). Sequencing was performed on the Illumina MiSeq sequencing platform with the V3 Illumina chemistry (2 × 300 cycles) at the ICM Institute (Paris, France; https://icm‐institute.org/fr/) according to standard protocols.

Sequence data processing

Raw fastq files were recovered and processed using QIIME2 software package. First, data were imported (‘qiime tools import’ command) into the software package and overall quality data were visualized (‘qiime demux summarize’ command). Quality filtering, denoising and merging of paired reads were performed using DADA2 with standard parameters, as implemented in QIIME2 (‘qiime dada2 denoise‐paired’ command; Callahan et al.,  2016). Reads were trimmed on their 5′ ends in order to remove primer sequences and truncated on their 3′ ends based on quality data (median quality < Q30). Thus, forward reads were truncated at 280 nt and reverse reads at 250 nt. The amplicon sequence variants (ASV) resulting from this step were then aligned using mafft (Katoh, 2002), and a midpoint‐rooted phylogeny was computed using FastTree (Price et al.,  2010). Alpha diversity parameters (Chao1, Faith phylogenetic diversity, Shannon indices, Observed OTUs, Goods coverage) were computed (‘qiime diversity alpha‐rarefaction’ command) and visualized using QIIME 2 View online interface (Table 2, Fig. S2). Estimates of the mean evolutionary diversity indices for the entire population were calculated based on ASV sequence alignments, for each sample, using MEGA 11 (Tamura et al.,  2021).

Taxonomic classification of ASV was performed using vsearch and the SILVA release 132 database (‘qiime feature‐classifier classify‐consensus‐vsearch’ command; Pruesse et al.,  2007; Rognes et al.,  2016). ASV DNA sequences, as well as a table grouping ASV counts per sample and full taxonomic paths, were then exported (‘qiime tools export’ command). In further steps, only ASVs affiliated to the ‘Discicristata’ taxonomic group were considered. In order to retain only heterolobosean sequences, sequences from Euglenozoa were filtered out, as this group constitutes a distinctive class from the Heterolobosea within the Discicristata eukaryotic clade. To improve taxonomic classification down to genus and or species level, a collection of 18S rRNA sequences from Heterolobosea representative was constituted (Table S5), covering the 18S rRNA region amplified in this study. Heterolobosean ASVs were queried against this custom database using the BLASTn algorithm (Altschul et al.,  1990). Output was manually curated and ASVs were assigned to defined genera and species when unambiguous results were obtained, using threshold defined by DNA distance matrices (Fig. 3). Results were visualized using ggplot2 in R environment (Wickham, 2009).

Raw sequences were deposited on the Sequence Read Archive and are freely available under project number PRJNA725645.

Conflict of interest

None declared.

Fig. S1 Positions the multiple primers designed and used in this study, placed on an alignment of heterolobosean 18S rRNA gene sequences, generated using PrimerDesignM. The x‐axis shows nucleotidic position, and y‐axes show entropy (light brown) and complexity (green). Primers are displayed as arrowheads, coloured according to the combinations used. *: F_628 was also used in combination with R_982.

Fig. S2: Rarefaction curves obtained for each sample, using Chao1, Goods coverage, Faith Phylogenetic diversity and the number of observed OTUs. The sampling depth was normalized at 20′000 sequences.

Fig. S3. Taxonomic assignments of amplicons yielded by primer pairs HF_1198/HR_1680 and HF_1202/HR_1680, from 4 environmental samples, using VSEARCH and the Protist Reference Database PR2 v4.14 (Guillou et al., 2013). Taxonomic assignments are shown at level 4 (A), showing abundances of heterolobosean sequences, and at level 9, showing genus and species‐level identification (B).

Table S1. Primer pairs to identify Heterolobosea communities by high‐throughput sequencing

Table S2. 18S rRNA sequences of Heterolobosea from SILVA SSU 132 database used for the design of new primer pairs

Table S3. Isolates of amoebae to evaluate the specificity of the primer pairs

Table S4. Impact of the various quality control and filtering steps on the high‐throughput sequencing dataset.

Table S5. 18S rRNA sequences of Heterolobosea from SILVA SSU 132 database used species identification.

Click here for additional data file.