Molecular Detection and Genotyping of Enteric Protists in Asymptomatic Schoolchildren and Their Legal Guardians in Madrid, Spain

Asymptomatic carriage of diarrhoea-causing enteric protist parasites in the general population is poorly understood, particularly in mediumto high-income countries. This molecular epidemiological survey investigates the presence, molecular diversity, and household transmission of Giardia duodenalis, Cryptosporidium spp., Blastocystis sp., and Enterocystozoon bieneusi in schoolchildren aged 2–13 years (n = 74) and their legal guardians (n = 6) in Madrid, Spain. Enteroparasite detection and genotyping was conducted in stool samples by molecular (PCR and Sanger sequencing) methods. Potential associations linked to infections were investigated through epidemiological questionnaires. Giardia duodenalis was the most prevalent enteric parasite found (14%, 95% CI: 7.1–23), followed by Blastocystis sp. (10%, 95% CI: 6.2–22) and Cryptosporidium spp. (3.8%, 95% CI: 0.78–11). None of the participants tested positive for E. bieneusi. Sequence analyses revealed the presence of G. duodenalis assemblage B, sub-assemblage BIV in a single child. The three Cryptosporidium isolates obtained were assigned to C. hominis, two of them belonging to the gp60 subtype IbA10G2. Four Blastocystis subtypes were identified including ST2 (38%, 3/8), ST3 (25%, 2/8), ST4 (25%, 2/8), and ST8 (12%, 1/8). All G. duodenalis and Cryptosporidium isolates were detected in children only. Blastocystis ST3 and ST4 were circulating in members of the same household. Blastocystis carriage rates increased with the age of the participants. Presence of diarrhoea-causing enteric protists was common in apparently


Introduction
Several enteric protistan species contribute to the burden of diarrhoeal illness globally. Those considered of higher public health relevance include the protozoa Cryptosporidium spp., Entamoeba histolytica, Giardia duodenalis, and Cyclospora cayetanensis, and the microsporidia Enterocytozoon bieneusi [1][2][3][4]. Infections by these well-recognized pathogens affect primarily young children and immunocompromised individuals in poor-resource settings characterised by insufficient access to safe drinking water and sanitation [5]. Other species such as the stramenopile Blastocystis sp. and the trichomonad Dientamoeba fragilis have a less-defined pathogenic potential [6,7]. In contrast with the relative abundance of data existing in endemic areas, a comparatively much lower amount of information is currently available on the epidemiology of diarrhoea-causing protistan pathogens in medium-Parasitologia 2021, 1 84 to high-income countries. In these settings, enteric protist infections are erroneously perceived as less relevant threats to public health because of better health standards, lower prevalence rates, less severe symptoms, and rarely fatal consequences [8]. This is in spite that some of these pathogens are frequent cause of food-and waterborne outbreaks of gastrointestinal disease globally [9,10], resulting in significant economic burden associated to medical and treatment costs [11]. Additionally, in developed countries, G. duodenalis and Cryptosporidium spp. infections have been associated with growth retardation and cognitive impairment in children [12] and with long-term complications including perceived food intolerance, irritable bowel syndrome, and chronic fatigue in the general population [13].
Giardia duodenalis, Cryptosporidium spp., Blastocystis sp., and E. bieneusi exhibit extensive intra-species genetic diversity, allowing the identification of several genotypes/subtypes with variable host specificities and ranges. Some of them can infect mammal species other than humans and are therefore zoonotic. Giardia duodenalis is currently regarded as a multispecies complex comprising eight (A to H) distinct assemblages, of which assemblages A and B cause the bulk of human infections. Infection cases by assemblages C-F have been sporadically reported mainly in children and immunocompromised individuals [14]. Cryptosporidium encompasses at least 45 valid species, with C. hominis and C. parvum accounting for more than 90% of documented human cases of cryptosporidiosis [15]. At least 22 subtypes (ST) have been identified within Blastocystis sp., with ST1-4 being responsible for most of the human cases reported globally; infections by ST5-9 and ST12 are only sporadically documented [16]. Finally, several hundred E. bieneusi genotypes have been defined and clustered into 11 phylogenetic groups, with Group 1 and Group 2 including those of zoonotic potential [17].
In the Autonomous Region of Madrid (ARM, Central Spain), the occurrence and genetic diversity of G. duodenalis and Cryptosporidium spp. was initially assessed in children attending day care centres in Majadahonda [18], and in symptomatic individuals of all ages seeking medical attention at two major hospitals in Madrid city [19,20]. A subsequent large-scale molecular epidemiological survey investigated these pathogens, in addition to Blastocystis sp. and E. bieneusi, in asymptomatic schoolchildren in Leganés [21,22]. These studies provided important insights in the understanding of the epidemiology of enteric protists in Spain including (i) G. duodenalis and Blastocystis (but not Cryptosporidium spp. or E. bieneusi) were a common finding in the stools of apparently healthy children, (ii) transmission was mainly of anthropic nature at local foci (e.g., schools), and iii) the occurrence of gastrointestinal symptoms (mainly diarrhoea) was not linked to a given pathogen species/genotype. The present study attempts to complement and expand previous knowledge in a different schoolchildren population in Madrid city. To assess the potential occurrence of within-household transmission events, legal guardians of participating schoolchildren were also investigated.

Occurrence of Protist Enteroparasites
A total of 80 stool samples and associated questionnaires were collected from asymptomatic schoolchildren (n = 74) and their legal guardians (n = 6). Of these, 44 (55%) were male and 34 (43%) were female. The sex of two participants (2.5%) was unknown. The age of the participating schoolchildren ranged from 2 to 13 years with a median age of 7 years. The age of the participating legal guardians ranged from 38 to 55 years with a median age of 44 years. The age of four participants was unknown. The number of participants by age group was as follow: 0-4 years, n = 14; 5-9 years, n = 46; 10-14 years, n = 10; >35 years, n = 6. Of the 80 participants, four (5.0%) were born in countries other than Spain including Bolivia (1), China (1), Cuba (1), and Ecuador (1). This variable was unknown for five (6.3%) participants. Only those PCR amplicons confirmed by Sanger sequencing were considered as positive for any of the enteric protist species investigated here. The full dataset used to determine the prevalence rates and molecular diversity of G. duodenalis, Cryptosporidium spp., Blastocystis sp., and E. bieneusi and the sociodemographic and clinical variables used to estimate risk infections in the present study are presented as a spreadsheet (Supplementary  Table S1).
The distribution of protist enteroparasites according to the sex and age group of the volunteering participants is shown in Figure 1. Giardia duodenalis (18% vs. 8.8%) and Blastocystis sp. (14% vs. 5.9%) were more prevalently found in males than in females, whereas the opposite was true for Cryptosporidium spp. (2.3% vs. 5.9%). Giardia duodenalis was detected in schoolchildren only, particularly in those aged 0-4 years. A similar trend was observed for Cryptosporidium spp. In contrast, Blastocystis sp. followed a marked age-related pattern, with individuals of older age bearing higher carriage rates.
participants. Only those PCR amplicons confirmed by Sanger sequencing were considered as positive for any of the enteric protist species investigated here. The full dataset used to determine the prevalence rates and molecular diversity of G. duodenalis, Cryptosporidium spp., Blastocystis sp., and E. bieneusi and the sociodemographic and clinical variables used to estimate risk infections in the present study are presented as a spreadsheet (Supplementary Table S1).
The distribution of protist enteroparasites according to the sex and age group of the volunteering participants is shown in Figure 1. Giardia duodenalis (18% vs. 8.8%) and Blastocystis sp. (14% vs. 5.9%) were more prevalently found in males than in females, whereas the opposite was true for Cryptosporidium spp. (2.3% vs. 5.9%). Giardia duodenalis was detected in schoolchildren only, particularly in those aged 0-4 years. A similar trend was observed for Cryptosporidium spp. In contrast, Blastocystis sp. followed a marked age-related pattern, with individuals of older age bearing higher carriage rates.

Molecular Characterization of G. duodenalis Isolates
The 11 DNA isolates that yielded a positive result for G. duodenalis by real-time PCR (qPCR) generated cycle threshold (Ct) values ranging from 18 to 37 (median: 33). Of these, 55% (6/11) had qPCR Ct values >30. Genotyping data were only available for a single isolate (from a nine years-old male) with a Ct = 18 that was successfully amplified at the glutamate dehydrogenase (gdh), ß-giardin (bg), and triose phosphate isomerase (tpi) genes of the parasite. Sequence analyses identified this isolate as assemblage B and sub-assemblage BIV (Table 1).

Molecular Characterization of G. duodenalis Isolates
The 11 DNA isolates that yielded a positive result for G. duodenalis by real-time PCR (qPCR) generated cycle threshold (Ct) values ranging from 18 to 37 (median: 33). Of these, 55% (6/11) had qPCR Ct values >30. Genotyping data were only available for a single isolate (from a nine years-old male) with a Ct = 18 that was successfully amplified at the glutamate dehydrogenase (gdh), ß-giardin (bg), and triose phosphate isomerase (tpi) genes of the parasite. Sequence analyses identified this isolate as assemblage B and sub-assemblage BIV (Table 1).

Molecular Characterization of Cryptosporidium spp. Isolates
All three Cryptosporidium-positive isolates detected by PCR amplification of the small subunit ribosomal RNA (ssu rRNA) gene were identified as C. hominis, and their partial sequences had 100% identity with reference sequence AF108865. Two of them yielded amplicons of the expected size at the PCR targeting the 60 kDa glycoprotein (gp60) gene of the parasite. Sequence analyses revealed the presence of gp60 subtype IbA10G2 in these samples, in both cases identical to reference sequence AY262031. Table 1. Diversity, frequency, and molecular features of Giardia duodenalis sequences at the glutamate dehydrogenase (gdh), ß-giardin (bg), and triose phosphate isomerase (tpi) loci generated in the present study. GenBank accession numbers are provided.

Risk Association Analysis
The results of the univariable analysis of potential risk associations between sociodemographic, epidemiological, clinical, and behavioural factors and a higher likelihood of infection/carriage by G. duodenalis, Cryptosporidium spp., or Blastocystis sp. are shown in Table 3. For Cryptosporidium spp., the only risk association determined was contact with family members with diarrhoea (p = 0.01). Despite the marked age-related pattern observed for Blastocystis sp. in Figure 1c (with individuals of older age bearing higher prevalence rates of the protist), this distribution was not statistically significant (p = 0.1). The same result was obtained when the four original age groups (0-4, 5-9, 10-14, and >35 years) were considered (X-squared = 1.9365, df = 3, p-value = 0.6). No risk associations could be demonstrated for G. duodenalis and Blastocystis sp. infections. Because of the low number of associations found, multivariable analysis was not attempted.

Discussion
This molecular-based survey provides novel data on the occurrence and molecular diversity of enteric protist species in Spanish paediatric populations, contributing to improve our knowledge on their epidemiology and transmission. Obtained results complement and expand those generated by our research group in the ARM. The main findings of the study are the confirmation that the protist enteroparasites G. duodenalis, Cryptosporidium spp., and Blastocystis sp. are a common finding in the stools of apparently healthy schoolchildren, and to a lesser extent, in their legal guardians. Evidence supporting the occurrence of within-household transmission was also provided for Blastocystis sp.
Giardia duodenalis was the most prevalent protist enteroparasite found in the present study (14%). Compared with previous PCR-based data available in the ARM, that infection rate was slightly lower than those reported in young children attending day care centres in Majadahonda (16%, 15/90) and in a similar schoolchildren population in Leganés (17%, 263/1,512) [18,22]. In the rest of the country, paediatric infections by G. duodenalis have been identified in the range of 3-36% by light and immunofluorescence microscopy, ELISA for the detection of copro-antigens or rapid diagnostic tests [23,24]. Regarding genotyping data, sequence analyses revealed sub-assemblage BIV in the only G. duodenalis isolate successfully genotyped in the present study. In the ARM, BIV was also the only G. duodenalis sub-assemblage detected in young children attending day care centres [18], and the most prevalent genetic variant of the parasite (79%) in a large school population [22]. The fact that BIV was also the predominant sub-assemblage (62%) in outpatients with clinical manifestations in this very same Spanish region strongly suggests that the genotype of G. duodenalis does not determine per se the outcome of the infection and the occurrence of symptoms [19].
The observed overall prevalence of Cryptosporidium infection (4%) was higher than those (1-3%) reported by PCR-based methods in similar studies conducted in the ARM [18,22]. Infections rates ranging from 1% to 10% have been previously documented in asymptomatic paediatric populations by non-molecular methods in other Spanish regions [24][25][26].
Not surprisingly, C. hominis was the only Cryptosporidium species found in the present study, affecting children only. Indeed, C. hominis has been reported as the most common (65-88%) Cryptosporidium species infecting symptomatic Spanish populations [12,27,28], including those (82%) surveyed in the ARM [20]. In addition, C. hominis caused 71% of the Cryptosporidium infections detected in asymptomatic schoolchildren in the ARM [22]. As in the case of the present survey, most of the above-mentioned C. hominis infections were attributed to the gp60 subtype IbA10G2. Taken together, all these molecular data indicate that the frequency of C. hominis infection is similar in all the Spanish populations investigated to date irrespectively of their clinical status.
In this survey 10% of the investigated children and adults carried Blastocystis sp., a prevalence rate lower than that (13%) identified also by PCR in an equivalent (but much larger) asymptomatic schoolchildren population in the ARM [22]. Both prevalence rates were in the low range of those (8-23%) previously reported by conventional microscopy in other Spanish populations including asymptomatic preschool-and schoolchildren [29], symptomatic outpatients [30], and HIV-positive children [31]. Our sequence analysis confirmed the presence of four Blastocystis STs, being ST2 the most prevalent (38%), followed by ST3 (25%), ST4 (25%), and ST8 (12%). Obtained results were remarkably similar to those (ST2: 36%, ST1: 23%, ST3: 22%, ST4: 19%, ST8: 1%) documented previously in schoolchildren in the ARM [22], with the exception that in the present survey, ST1 was not detected. Interestingly, ST2 (62%) was also the most prevalent Blastocystis ST identified in a community survey conducted in northern Spain [32]. An important contribution of this study is the demonstration that the distribution of Blastocystis sp. followed a clear age-related pattern, with individuals of older age bearing the highest carriage rates. This finding may be the direct consequence of the persistent enteric colonization of the protist, suggesting that Blastocystis sp. should be primarily regarded as a non-pathogenic commensal that forms part of the gut microbiota of the host. Similar age-related patterns for Blastocystis sp. carriage have been observed in healthy children in Spain [21], in orphan children and their caregivers in Thailand [33], and in symptomatic and asymptomatic schoolchildren in Mozambique [34]. In addition, the finding of ST3 and ST4 in a father and his daughter sharing the same household highlights the relevance of asymptomatic carriers in the spreading of faecally-orally microorganisms including Blastocystis sp.
Confirming previous results in apparently healthy schoolchildren populations [22], the microsporidia E. bieneusi was not detected in the surveyed paediatric and adult population. In Spain, E. bieneusi infections have been identified mainly in immunocompromised individuals including HIV/AIDS patients [31,35], transplant recipients [36], and elderly people [37], indicative of the opportunistic nature of this pathogen.
In this study, no animal-adapted species/genetic variants of G. duodenalis and Cryptosporidium spp. were identified, suggesting that the main source of infections by these pathogens should be anthropic in nature, most likely through direct human-to-human contact. However, the finding of Blastocystis ST8 in a 9-year-old male with no obvious risk factors for parasite infections was interesting. This very same Blastocystis ST was previously detected in a 10-year-old female in the ARM [22]. ST8 carriage has been documented primarily in captive and wild non-human primates, but also in swine, several avian species, and rodents [38]. Additionally, zoonotic transmission of ST8 has been documented between captive primates and their handlers in a zoological garden in the United Kingdom [39]. In our study, the source of infection by ST8 was unknown, but a zoonotic transmission event was highly suspected.
The results obtained in the present study and some of the conclusions reached may be biased by some design and methodological limitations. The most important was related to the reduced sample size (particularly those from legal guardians) that have affected the statistical power of the analyses conducted. Low G. duodenalis burdens (evidenced by the elevated qPCR Ct values observed) compromised the success rate of the PCR protocols used for genotyping purposes, as the amplified markers (gdh, bg, and tpi) were all single-copy genes with limited diagnostic sensitivity. In addition, our qPCR protocol lacked amplification internal control. Therefore, it is possible that an unknown number of true Giardia-positive samples were undetected due to the inhibition of the amplification reactions.

Study Area and Stool Sample Collection
A cross-sectional molecular epidemiological study of diarrhoea-causing enteric parasites including the protozoan G. duodenalis and Cryptosporidium spp., the stramenopile Blastocystis sp., and the microsporidia E. bieneusi was carried out in voluntary asymptomatic schoolchildren (2-13 years) attending a public school in Madrid (central Spain) in October 2017. An informative meeting was held for interested families and parents/legal guardians were also invited to participate; participating schoolchildren and parents/legal guardians were provided with sampling kits (sterile polystyrene plastic flask with spatula and a unique identification number) to obtain individual stool samples. Parents/legal guardians assisted in the collection of stool samples from consenting schoolchildren and brought the samples to school. Samples were collected by members of the research team at scheduled times, transported to the Spanish National Centre for Microbiology, and stored at 4 • C. All stool samples were processed within 3 days of collection.

DNA Extraction and Purification
Genomic DNA was isolated from about 200 mg of each faecal specimen using the QIAamp DNA Stool Mini Kit (Qiagen, Hilden, Germany) according to the manufacturer's instructions, except that samples mixed with InhibitEX buffer were incubated for 10 min at 95 • C. Extracted and purified DNA samples (200 µL) were kept at 4 • C until further molecular analysis.

Molecular Detection and Characterization of Giardia duodenalis
A qPCR protocol amplifying a 62 bp-fragment of the ssu rRNA gene of Giardia duodenalis was used as initial screening test for the presence of the parasite [40]. Amplification reactions (25 µL) consisted of 3 µL of template DNA, 0.5 µM of primers Gd-80F and Gd-127R, 0.4 µM of probe (Supplementary Table S2), and 12.5 µL TaqMan ® Gene Expression Master Mix (Applied Biosystems, Foster City, CA, USA). Detection of parasitic DNA was performed on a Corbett Rotor Gene™ 6000 real-time PCR system (Qiagen) using an amplification protocol consisting on an initial hold step of 2 min at 55 • C and 15 min at 95 • C followed by 45 cycles of 15 s at 95 • C and 1 min at 60 • C. Water (no-template) and genomic DNA (positive) controls were included in each PCR run.
Giardia duodenalis isolates with a qPCR-positive result were re-assessed by sequencebased multi-locus genotyping of the single-copy genes gdh, bg, and tpi of the parasite. A semi-nested PCR was used to amplify a 432-bp fragment of the gdh gene [41]. PCR reaction mixtures (25 µL) included 5 µL of template DNA and 0.5 µM of the primer pairs GDHeF/GDHiR in the primary reaction and GDHiF/GDHiR in the secondary reaction (Supplementary Table S2). Both amplification protocols consisted of an initial denaturation step at 95 • C for 3 min, followed by 35 cycles of 95 • C for 30 s, 55 • C for 30 s, and 72 • C for 1 min, with a final extension of 72 • C for 7 min. A nested PCR was used to amplify a 511 bp-fragment of the bg gene [42]. PCR reaction mixtures (25 µL) consisted of 3 µL of template DNA and 0.4 µM of the primers sets G7_F/G759_R in the primary reaction and G99_F/G609_R in the secondary reaction (Supplementary Table S2). The primary PCR reaction was carried out with the following amplification conditions: one step of 95 • C for 7 min, followed by 35 cycles of 95 • C for 30 s, 65 • C for 30 s, and 72 • C for 1 min with a final extension of 72 • C for 7 min. The conditions for the secondary PCR were identical to the primary PCR, except that the annealing temperature was 55 • C. Finally, a nested PCR was used to amplify a 530 bp-fragment of the tpi gene [43]. PCR reaction mixtures (50 µL) included 2-2.5 µL of template DNA and 0.2 µM of the primer pairs AL3543/AL3546 in the primary reaction and AL3544/AL3545 in the secondary reaction (Supplementary Table S2). Both amplification protocols consisted of an initial denaturation step at 94 • C for 5 min, followed by 35 cycles of 94 • C for 45 s, 50 • C for 45 s, and 72 • C for 1 min, with a final extension of 72 • C for 10 min.

Molecular Detection and Characterization of Cryptosporidium spp.
The presence of Cryptosporidium spp. was assessed using a nested PCR to amplify a 587-bp fragment of the ssu rRNA gene of the parasite [44]. Amplification reactions (50 µL) included 3 µL of DNA sample and 0.3 µM of the primer pairs CR-P1/CR-P2 in the primary reaction and CR-P3/CPB-DIAGR in the secondary reaction (Supplementary Table S2). Both PCR reactions were carried out as follows: one step of 94 • C for 3 min, followed by 35 cycles of 94 • C for 40 s, 50 • C for 40 s, and 72 • C for 1 min, concluding with a final extension of 72 • C for 10 min. Sub-typing of the isolates identified as C. hominis was attempted at the gp60 gene. Briefly, a nested PCR was conducted to amplify an 870-bp fragment of the gp60 locus [45]. PCR reaction mixtures (50 µL) included 2-3 µL of template DNA and 0.3 µM of the primer pairs AL-3531/AL-3535 in the primary reaction and AL-3532/AL-3534 in the secondary reaction (Supplementary Table S2). The primary PCR reaction consisted of an initial denaturation step of 94 • C for 5 min, followed by 35 cycles of 94 • C for 45 s, 59 • C for 45 s, and 72 • C for 1 min with a final extension of 72 • C for 10 min. The conditions for the secondary PCR were identical to the primary PCR except that the annealing temperature was 50 • C.

Molecular Detection and Characterization of Blastocystis sp.
Identification of Blastocystis sp. was achieved by a direct PCR targeting the ssu rRNA gene of the parasite [46]. This protocol uses the pan-Blastocystis, barcode primers RD5 and BhRDr (Supplementary Table S2) to amplify a PCR product of 600 bp. Amplification reactions (25 µL) included 5 µL of template DNA and 0.5 µM of the primer set RD5/BhRDr. Amplification conditions consisted of one step of 95 • C for 3 min, followed by 30 cycles of 1 min each at 94, 59, and 72 • C, with an additional 2 min final extension at 72 • C.

Molecular Detection and Characterization of Enterocytozoon bieneusi
Enterocytozoon bieneusi DNA was detected using a nested PCR to amplify a 390 bpfragment of the entire internal transcribed spacer as well as portions of the flanking large and small subunits of the rRNA gene of the parasite [47]. Reaction mixtures (50 µL) contained one µL of DNA template and 0.2 µM of the primer pairs EBITS3/EBITS4 in the primary reaction and EBITS1/EBITS2.4 in the secondary reaction (Supplementary Table S2). After denaturation at 94 • C for 3 min, amplification conditions of the primary PCR included 35 cycles of amplification (denaturation at 94 • C for 30 s, annealing at 57 • C for 30 s, and elongation at 72 • C for 40 s), followed by a final extension at 72 • C for 10 min. Conditions for the secondary PCR were identical to the primary PCR except only 30 cycles were carried out with an annealing temperature of 55 • C.
All the direct, semi-nested, and nested PCR protocols described above were conducted on a 2720 thermal cycler (Applied Biosystems). Reaction mixes included 2.5 units of MyTAQ™ DNA polymerase (Bioline GmbH, Luckenwalde, Germany), and 5× MyTAQ™ Reaction Buffer containing 5 mM dNTPs and 15 mM MgCl 2 . Laboratory-confirmed positive and negative DNA isolates for each parasitic species investigated were routinely used as controls and included in each round of PCR. PCR amplicons were visualized on 2% D5 agarose gels (Conda, Madrid, Spain) stained with Pronasafe nucleic acid staining solution (Conda). Positive-PCR products were directly sequenced in both directions using the internal primer set described above. DNA sequencing was conducted by capillary electrophoresis using the BigDye ® Terminator chemistry (Applied Biosystems) on an ABI PRISM 3130 automated DNA sequencer.

Sequence Analyses
Raw sequencing data in both forward and reverse directions were viewed using the Chromas Lite version 2.1 (Technelysium Pty Ltd., South Brisbane, Australia) sequence analysis program (https://technelysium.com.au/wp/chromas/). The BLAST tool (http://blast.ncbi. nlm.nih.gov/Blast.cgi) was used to compare nucleotide sequences with sequences retrieved from the NCBI GenBank database. Generated DNA consensus sequences were aligned to appropriate reference sequences using the MEGA version 6 software [48] to identify Giardia species and assemblages/sub-assemblages and Cryptosporidium species. Blastocystis sequences were submitted at the Blastocystis 18S database (http://pubmlst.org/blastocystis/) for sub-type confirmation and allele identification. The sequences obtained in this study have been deposited in GenBank under accession numbers MW810324-MW810326 (G. duodenalis), MW789158 and MW810327 (Cryptosporidium spp.), and MW789153-MW789157 (Blastocystis sp.).

Epidemiological Questionnaire
A standardised questionnaire was provided as part of the sampling kit to be completed by the children's legal guardians. Questions included: (i) demographic characteristics, e.g., age, sex, country of birth, and number of siblings; (ii) behavioural habits, e.g., hand and fruit/vegetable washing and whether there have been any occurrence of diarrhoea in the participant, their family members, their schoolmates, and/or pets; and (iii) additional questions on potential risk factors, e.g., types of drinking water, whether they had swum in pools or natural waters in the 2 weeks prior to sample collection, had any contact with pets and any recent travel abroad. Completed questionnaires, signed informed consents, and individual stool samples were returned for collection by each participating schoolchild as described above.

Statistical Analyses
Data entry from the epidemiological questionnaires was completed using EpiData version 4.2.0 (EpiData Association, Odense, Denmark) and analysed using R software (R Foundation for Statistical Computing, Vienna, Austria). Prevalence of any G. duodenalis, Cryptosporidium spp., and Blastocystis sp. infections (single or multiple) in the study population was calculated. The chi-squared and/or Fisher's exact test were used to compare G. duodenalis, Cryptosporidium spp., and Blastocystis sp. infection rates by sex and categorical risk factors (including the other protists as a possible risk factor). A p value < 0.05 was considered evidence of statistical significance.