Finding Sarcocystis spp. on the Tioman Island: 28S rRNA gene next-generation sequencing reveals nine new Sarcocystis species

Sarcocystosis is a parasitic disease caused by members of the Sarcocystis genus. With more than 200 species, an intermediate–definitive host (diheteroxenous) life cycle (Fayer 2004), and occasionally being dihomoxenous (Bannert 1994; Koudela & Modrý 2000), the Sarcocystis genus, is ubiquitous in various environments and hosts. Numerous findings of sarcocystosis among wildlife and livestock have been reported (Latif et al. 1999; Dubey et al. 2015; Gjerde et al. 2018). In humans, it can take the form of intestinal or muscular sarcocystosis.

Between 2011 and 2014, a series of acute muscular sarcocystosis had occurred among tourists who travelled to the tropical islands of Tioman and Pangkor, Malaysia (CDC 2012; Esposito et al. 2012; Esposito et al. 2014; Italiano et al. 2014). With confirmation from just a handful of biopsy samples, S. nesbitti was found to be the only known species to cause human muscular sarcocystosis on these two islands (Tappe et al. 2013; Lau et al. 2014). Water and food are stipulated to be the transmission routes of sarcocystosis on the Malaysian islands (Dubey 2015). Early surveillance monitoring was carried out in Tioman, but finding Sarcocystis among the environmental samples was elusive (Husna Maizura et al. 2012). A subsequent study in Tioman that examined various water samples taken from Salang, Juara and Ayer Batang reported the detection of S. singaporensis, S. nesbitti, Sarcocystis sp. YLL-2013 and one unknown Sarcocystis species based on the 18S rRNA gene sequences (Shahari et al. 2016).

This study continues the quest of finding Sarcocystis in Tioman, using next-generation sequencing (NGS) to target a region in the 28S rRNA gene with a new primer set 28S R7F–28S R8 Deg R, and the V9 hypervariable region of the 18S rRNA gene. With the deep sequencing ability of NGS and identification of unique 28S rRNA gene sequences, the presence of more Sarcocystis species in the environmental samples from Tioman was revealed. The sequence of primary polymerase chain reaction (PCR) products from initial screening examinations, and cloning results when available are referred in the reporting of these NGS acquired sequences. In addition, the corresponding amplified regions (R7–R8 and V9) belonging to eight in-house unknown Sarcocystidae sequences with long nucleotides from the 18S to 28S rRNA gene were extracted for comparison to assign a full-length ribosomal DNA sequence to the matching NGS sequence.

Samplings were conducted in four main locations on the Tioman Island, L1 to L4, in October 2014 and August 2015. Nine and 23 grab water samples and 19 and six grab soil samples were collected, respectively, in the two consecutive years, totalling up to 32 grab water samples and 25 grab soil samples. Samples were acquired from rivers, estuaries, mangrove and bay. The samples were kept in sterile Nasco Whirl-Pak bags, and stored on ice in cool boxes or refrigerated until being processed. Samples acquired in 2014 were processed within two months, and within two weeks for the 2015 samples. DNA was extracted from each water sample in duplicates; each replicate with 100 ml of water content aspirated from the sampling bags, filtered through 0.45 μm cellulose nitrate membrane filters (Sartorius, Germany) using the Millipore vacuum manifold system, and followed by genomic DNA extraction with PowerWater DNA Isolation Kit (Mo Bio, USA). DNA was extracted from each soil sample in triplicates; the first and second replicates each with genomic DNA extracted from 0.25 g of soil sample using the DNeasy PowerSoil DNA Isolation Kit (Mo Bio, USA), and the third replicate from 10 g of soil samples using the PowerMax Soil DNA Isolation Kit (Mo Bio, USA). Extracted DNA was stored at −20 °C. A muscle sample from water monitor lizard provided by The Tioman Island Sarcocystosis Investigation Team (Esposito et al. 2014) served as the positive control in this study.

One domestic ceramic water filter was acquired from a household during the 2014 sampling trip. The device filtered natural surface water and groundwater collected for household usage. It was soaked overnight in a 1 L glass measuring cylinder containing 400 ml of sterilised artificial groundwater diluent. The diluent contained 0.284 g MgSO₄·7H₂O, 0.032 g NaCl, and 0.266 g CaCl₂·2H₂O per litre of ultrapure water (Charles et al. 2009). The content in the cylinder was transferred into a sterile stainless steel container. The debris (filtered material) on the water filter was scrapped into the soaked content using a sterile stainless steel spatula. The content was then divided into eight 50 ml Corning centrifuge tubes and sonicated for 20 min to dislodge microorganisms that possibly had clumped up among the silt, followed by ultracentrifugation of the tubes at 13,000 g. The supernatant of two tubes was pooled together for genomic DNA extraction using the PowerWater DNA Isolation Kit (Mo Bio, USA). Soil pellet from each two centrifuge tubes was pooled into one PowerBead tube of the PowerMax Soil DNA Isolation Kit (Mo Bio, USA) using the PowerBead Solution, followed by genomic DNA extraction.

The extracted DNA of all the technical replicates from their respective samples were screened for the presence of Sarcocystidae through PCR, using a newly designed primer pair (Primer3Plus, http://www.bioinformatics.nl/cgi-bin/primer3plus/primer3plus.cgi) consisting of forward primer 28S R7 F (5′-TCGCGTTCATCGGGATTTGA-3′) and reverse primer 28S R8 Deg R (5′-TTTCAAGACGGGTCGGTTRA-3′). The expected product size is about 480 bp, based on Sarcocystidae sequences retrieved from the GenBank database. The PCR contained 25 μl of Q5 High Fidelity 2X Master Mix (New England Biolabs), 2.5 μl (10 μM) of each forward and reverse primer, 3 μl of genomic DNA, and 17 μl of Nuclease-Free Water (Qiagen, USA). PCR was performed with the following conditions: 98 °C for 3 min of initial denaturation, 40 cycles of amplification at 98 °C for 10 s, 65 °C for 45 s, 72 °C for 45 s; final elongation at 72 °C for 3 min and hold at 4 °C. The PCR products were screened with FlashGel DNA Cassette (1.2% agarose, Lonza, Switzerland). Presumptively positive samples (showing bands of about 500 bp in size) were electrophoresised on 1.5% agarose gel (HydraGene Instant Agarose Tablet, USA). The respective bands were excised and sent for Sanger sequencing (Genomics, Taiwan) to confirm the presence of Sarcocystidae. An environmental sample with at least one technical replicate confirmed as harbouring Sarcocystidae by Sanger sequencing was considered positive for Sarcocystidae.

PCR using primers 28S R7F–28S R8 Deg R that contained the Illumina NGS platform adapter, 5′-TCGTCGGCAGCGTCAGATGTGTATAAGAGACAGTTTCAAGACGGGTCGGTTRA-3′ and 5′-GTCTCGTGGGCTCGGAGATGTGTATAAGAGACAGTCGCGTTCATCGGGATTTGA-3′ were carried out. About 36 μl of each of the PCR products was electrophoresised. Bands of the expected size were excised, purified with Monarch DNA Gel Extraction Kit (New England Biolabs) and sequenced to confirm that Sarcocystidae has been amplified. The purified PCR products were then subjected to NGS library preparation using the Nextera XT Index Kit (Illumina, USA) and HiFi HotStart ReadyMix PCR Kit (Kapa Biosystems, USA), undergoing a reduced PCR cycle of eight. The NGS libraries were then quantified with the Kapa Library Quantification Kit (Kapa Biosystems, USA) using the Applied Biosystems 7500 Real-Time PCR Systems, normalised to 1.4 nM working concentration, and loaded into the Illumina Miseq System following the manufacturer's instruction, using the Miseq Reagent Nano Kit, v2 (300 cycles) and the PhiX Control Kit v3 (Illumina, USA). The reverse primer 28S R8 Deg R was assigned to the Read 1 adapter to enable single-end reading from the antisense direction, as public database information reveals a more diverse region from this direction among Sarcocystidae. For quality assurance, four replicates were also sent for paired-end 500 cycles NGS (BioBasic, Canada).

Ten candidates were sent for paired-end NGS of the V9 hypervariable region in the 18S rRNA gene (1st Base, Singapore), using the 1380F and 1510R universal primers. Seven of these 10 candidates had also been used for 28S rRNA gene NGS, either as the same technical replicate or different replicates of the same sample.

PCR product cloning was performed using either the pGC Blue Cloning & Amplification Kits (Lucigen, USA) with the E. cloni 10G Supreme Electrocompetent Cells, or NEB PCR Cloning Kit with NEB 10-Beta chemically competent cells (New England Biolabs). Samples R-2015, J-2015, Sarco, H-2015 and I-2015 were cloned through technical replicates R2-2015, J2-2015, Sarco, H2-2015 and I2-2015, and the products were sent for Sanger sequencing (Genomics, Taiwan). The first three technical replicates have also undergone 28S NGS, and replicate H1-2015 for 18S V9 NGS.

There were eight in-house unknown Sarcocystidae sequences, namely Combi, WFMax1B, WF2B, C1.9, C2.1, C2.5, 50.20Jun and J1.11. These sequences were acquired through the cloning of PCR products that amplified genomic DNA from the 18S to the 28S rRNA gene (submitted), using the same samples from this study. Their nucleotide fragments delineated by primers 28S R7F–28S R8 Deg R and 1380F–1510R were extracted for matching with the NGS reads.

The NGS data output was analysed using the USEARCH commands available from https://www.drive5.com/usearch/manual/cmds_all.html (accessed between July 2017 and April 2018). The decompressed fastq files were filtered with the criteria of maximum expected error rate (maxee) of 1.0, truncated into a suitable length of nucleotides and saved as fasta files. Paired-end fastq files were merged first before the similar quality filtration step. All quality filtered fasta files belonging to the 28S rRNA gene were pooled into a single fasta file to produce unique reads from the total libraries, followed by generation of Zotus (Zero-radius Operational Taxonomic Unit) or Otus (Operational Taxonomic Unit), and mapping of the respective raw fastq files to the Zotus or Otus. Similar steps were performed for the paired-end reads of 18S V9 region. The command lines for NGS data processing from fastq files to alpha and beta analysis based on Zotus generation are presented in Supplementary Material S1 (available with the online version of this paper). Otus were produced using the cluster_otus command (not shown).

Sequence display, alignment, analysis and editing were carried out using Mega 7 (Kumar et al. 2016), Clustal Omega (https://www.ebi.ac.uk/Tools/msa/clustalo/), NCBI Tree Viewer (https://www.ncbi.nlm.nih.gov/projects/treeview/) and iTOL (Letunic & Bork 2011). Four alpha diversity metrics were reported, namely Richness, Chao1, Simpson and Berger-Parker. The Richness metrics of the amplified 28S and 18S regions are portrayed as rarefaction curves with Excel. Bray-Curtis that is based on abundance was used as the main beta diversity analysis metric, supplemented by the Jaccard metric that is based on presence/absence. The phylogenetic tree was produced from aligned sequences with the Neighbour Joining method and 100 bootstrap replicates.

Sarcocystidae prevalence among the environmental samples acquired between 2014 and 2015 showed that 65.6% of the 32 grab water samples (5/9 in 2014 and 16/23 in 2015) and 28% of the 25 grab soil samples (6/19 in 2014 and 1/6 in 2015) contained Sarcocystidae. Note that the water filter did not count towards these prevalence statistics. Due to resource limitation, not all Sarcocystidae-positive technical replicates were subjected to NGS. Eighteen samples were chosen for the 28S rRNA gene NGS library preparation, producing 25 libraries, among which were three pairs of sample replicates, J1-2015 and J2-2015; S1-2015 and S2-2015; Sarco_S10 and Sarco.mar; and one group of related replicates that were processed from the water filter, namely WF-1A, WF-1B, WF-2A, Dmax.mar and WF-Sup. The reads and diversity metrics of these libraries are presented in Table 1. Reads from libraries D7.mar, Dmax.mar, Sarco.mar and D2.mar (paired-end sequencing) were trimmed down to 480 bp, whereas reads from the other 21 libraries (single-end sequencing) were trimmed down to 290 bp. Libraries Sarco_S10 (300 cycles) and Sarco.mar (500 cycles) served as the positive controls. D2.mar served as a ‘negative control’, as the Sanger sequencing result of the PCR product (with Illumina adaptors) was ineligible.

85.8% of the 28S rRNA gene NGS reads achieved a Phred quality score of Q30 and above. The number of raw reads pooled from all the 25 libraries was 447,328, which produced a total of 283,933 reads after applying the filtration criteria of maxee 1.0. This subsequently yielded 39,854 unique reads; among them 30,903 reads (77.5%) were singletons. Sixty-two Zotus were produced from the unique reads (Supplementary Material S2, available with the online version of this paper), of which 15 chimeras were removed. A total of 367,095 raw reads were successfully mapped to the Zotus. At least 84.2% of raw reads per sample could be mapped to the Zotus. In comparison, from the 39,854 unique reads, 33 Otus were produced, whereby 65 chimeras were removed (results not shown). The rarefaction curves of the 28S rRNA gene NGS were presented as log₁₀ trend lines of the richness index, as shown in Figure 1. From the figure, it seems that the 1.4 nM of NGS libraries concentration could still be increased under this experimental workflow, but is nonetheless sufficient.

The use of Zotu (UNOISE algorithm) was considered more appropriate in this study as compared to Otu (UPARSE algorithm) because the heterogeneity of sequences in the amplified 28S region is unknown. Hence, it is better to apply the ‘all correct biological sequences are assigned’ principle of the UNOISE algorithm, rather than the ‘no two sequences are >97% identical’ principle of the UPARSE algorithm (Edgar & Flyvbjerg 2015; Edgar 2016). The 97% similarity threshold could be insufficient to differentiate Sarcocystis species (Edgar 2018). On the other hand, a sequence that is produced by both the UNOISE and UPARSE algorithm may be more credible as a correct biological sequence.

The highest observed diversity among the libraries, as portrayed by the ‘Richness’ index, was reported by the sample replicate pair of J1-2015 and J2-2015, and also T2-2015, which were water samples acquired from two different locations, L1 and L3 (Table 1). However, the ‘Chao 1’ diversity index that also captures low-frequency reads in a library, estimated that P2-2014, C3-2014 and R2-2015, the first two being soil samples, would have the highest diversity. Based on the ‘Richness’ and ‘Chao1’ indices, water samples 2A-2014, 2A-2015, 8B-2014, V1-2015 and D7.mar were constantly reported as having low diversity. This implies the presence of dominant species in the samples, which is also reflected by their Simpson values of about 0.4 or higher.

Based on the Bray-Curtis beta analysis of the libraries as being depicted in Figure 2, water sample 8B-2014 that was acquired from location L1, which had common sightings of water monitor lizards, was closely related to the muscle sample of water monitor lizard (Sarco_S10 and Sarco.mar, positive controls). Two water samples, V1-2015 and D7.mar, were clustered under the same clade with four replicates of the water filter, WF-1A, WF-1B, WF-2A, and Dmax.mar. Interestingly, V1-2015, which was acquired from a mangrove area in L3, had the most similar content to WF-2A, surpassing the other three replicates of the water filter, which was taken from location L1 in 2014. Furthermore, the supernatant content of the water filter, WF-Sup, were related to 2A-2015, a soil sample taken from L2 in 2015. The Jaccard beta analysis of the libraries presented similar relationships as the Bray-Curtis analysis of Figure 2 (results not shown). These overall findings suggest that there is no specific clustering of Sarcocystidae based on different locations of the island, over the two subsequent years of 2014 and 2015. Nevertheless, sampling from the mangrove area and water filter appears to be useful, besides sampling different sections of rivers.

The potential identity of the 62 Zotus produced from the 25 libraries and their frequency of presence are shown in Supplementary Material S3 (available online). Forty-six of the Zotus made up five main ‘clusters’ of Zotus, showing various percentages of similarity to Sarcocystis sp. No. 5, S. singaporensis, S. speeri, S. zuoi, and S. zamani. Among them, Zotu 28 and Zotu 41 matched perfectly to S. zamani (GenBank KU244528.1) and S. singaporensis (GenBank KU341123.1), respectively. The identity of the remaining Zotus could be S. cymruensis, S. lari, S. halieti, Eimeria gruis and Goussia janae. The Zotus related to S. singaporensis, S. zamani and S. cymruensis showed at least 90% of similarities, with 100% query coverage, except for Zotu 47, 59, and 7 that had lower query coverage. Other Zotus have a poor similarity of below 90% to their related species. If the six Zotus that returned no matches in GenBank Megablast search and Zotu.20 were excluded, it is possible that this study has identified 53 unique Sarcocystidae species (Zotu.28 and Zotu.41 not included) through the NGS approach. Among these unique Zotu sequences, 36 had a similarity of less than the 97% threshold to their related species, and thus likely are different Sarcocystidae species.

There were a total of 67 Sarcocystis species with 28S rRNA gene information in the Genbank database as of July 2018. All these were initially used to produce a phylogeny tree together with the 62 Zotus. Species that did not show additional phylogenetical information were removed, leaving 17 representative species. The phylogenetic relationships of these 17 Sarcocystidae species, 62 Zotus and eight in-house unknown Sarcocystidae sequences (R7–R8 region of the 28S rRNA gene, Supplementary Material S2) are delineated in Figure 3. Three Zotus turned out to be identical: Zotu.27 to J1.11, Zotu.12 to Combi, and Zotu.58 to 50.20Jun.

The ‘cluster’ of Zotus related to Sarcocystis singaporensis (Supplementary Material S3) was phylogenetically divided into one monophyletic group consisting of Zotu 36, 56, 38, 60, 25, 10, 11, 4, 23, 8, 49, and 27, which also included J1.11 and S. singaporensis KU341122.1; and polyphyletic groups of 16 other Zotus (Zotu 17, 37, 41, 43, 44, 13, 59, 33, 40, 35, 51, 21; 47, 24, 34, 62) and S. singaporensis KU341123.1 (Figure 3). The ‘cluster’ of Zotus related to Sarcocystis sp. No. 5 (Supplementary Material S3) was grouped into a single monophyletic group consisting of Zotu 26, 55, 42, 50, 3, 45, 32, 1, 14, together with WF2B. Zotu 5, 58, 39 and 12 are in the same monophyletic group with 50.20Jun and Combi that also includes Sarcocystis species isolated from mammals and ducks. This group shares the same clade with another monophyletic group that comprises Zotu 31 and 61, which were related to Sarcocystis species of viper and cattle origins. Zotu 22 that has a poor query coverage with Megablast search appears together with three other Zotus that have no result in Megablast search, the Zotu 15, 48 and 54, in the monophyletic group that also comprises of Zotu 46, 53 and 19, C1.9, C2.1, C2.5, and the paraphyletic group of Zotu 6 and 18. These 12 taxons may be related to non-SarcocystisSarcocystidae species.

Eight Zotus have sequences identical to primary PCR products; four are identical to cloned sequences; three appeared both as primary PCR product and clones. These 15 Zotus with their presence verified by additional Sanger sequencings are in bold in Figure 3. The four sample replicates (D7.mar, Dmax.mar, Sarco.mar and D2.mar) subjected to paired-end (500 cycles) NGS were re-analysed as a data subset. This produced seven Zotus with 480 bp. Five of them matched to Zotu 12, 31, 42, 58 and 61 at 290 bp length. They were denoted as Zotu.n.480 in parenthesis next to the respective Zotu headings in Supplementary Material S2. The remaining two Zotus with 480 bp, Zotu5.480 and Zotu7.480, were reported as the last entries in the same reference.

Based on the frequency of presence, abundance, and identity uniqueness (Nguyen et al. 2015), nine Zotus are named, as shown in Table 2. Each of these Zotus has no higher than 98% similarity to a known species (Supplementary Material S3). The multiple sequence alignment of the named Zotus is reported in Supplementary Material S4 (available online). Sarcocystis varanus (Zotu 12) and Sarcocystis biawaki (Zotu 58) were a ‘doppelganger’ pair. Both were found as dominant sequences of the positive control (sample replicate S10 and Sarco.mar), which also appeared in samples from two different rivers of location L1, and in one of the water filter replicate (D.max.mar) (refer Supplementary Material S5, available online). The ‘doppelganger’ pair has only one nucleotide of difference, but their simultaneous presences in other samples and replicates suggest that their identities are best kept distinguished. The Otu algorithm tolerates a higher level of nucleotide difference, thus did not distinguish the pair (Table 2). Sarcocystis J1.11 (Zotu 27) and Sarcocystis rompini (Zotu 8) were clustered under the same monophyletic group (Figure 3), with a difference of three nucleotides. Seven of the named Zotus were detected in primary PCR and/or cloning, and hence supplemented with full-length sequence of 480 bp (Table 2, Supplementary Material S2). The remaining two, Sarcocystis malayani (Zotu 13) and Sarcocystis pahangi (Zotu 14), were at least produced by both the Zotu and Otu algorithm, and have only 290 bp of sequence length.

Table 3 presents the NGS sequencing results of the 18S rRNA gene V9 hypervariable region of 10 environmental samples. 71.6% of the reads achieved a Phred quality score of Q30 and above. The total number of merged reads from the 10 samples was 2,003,590. This produced 132,734 unique reads; among them 81,160 reads (61.1%) were singletons. 6,995 Zotus were yielded, with 15 chimeras removed. Values of the Richness metrics and the estimated Chao1 are very close, indicating saturation of detection, as also displayed by the rarefaction curves in Figure 4.

Due to the abundance of unique reads produced by the 18S V9 NGS, a nucleotide strand of 51 bp (CCCTGCCCTTTGTACACACCGCCCGTCGCTCCTACCGATTGAGTGTTCCGG) was used to narrow down the search for reads related to Sarcocystis. The length of this strand was selected based on the eight in-house unknown Sarcocystidae sequences. Among the 132,734 unique sequences, 15 sequences matched with this 51 bp ‘screening’ nucleotide strand, of which only two are Sarcocystis species (Supplementary Material S6, available online). They are Uniq.V9.21756, which is 99% similar to S. nesbitti (HF544323.1) with 98% query cover, and Uniq.V9.39999, that is 91% similar to Sarcocystis fulicae with 100% query cover. Uniq.V9.21756 reported 105 reads in WFMax1B and 3 reads in P3A-2014; Uniq.V9.39999 had 34 reads in WFMax1B. These 15 sequences in the V9 region (160 bp), and the sequences of eight in-house Sarcocystidae species in the similar region are reported in Supplementary Material S7 (available online).

The difference is clear between the numbers of Sarcocystis species detected based on sequences of the 28S rRNA gene as compared to the sequences of V9 region of the 18S rRNA gene. That is, 62 Zotus results in potentially 50 Sarcocystis species (not including non-Sarcocystis Sarcocystidae and Zotus without search result by Megablast, Supplementary Material S3), versus 6,995 Zotus leading to the identification of only two Sarcocystis species. This shows surpassing sensitivity and specificity of the 28S R7F–28S R8 Deg R primer pair in targeted Sarcocystis amplification compared to the 18S V9 region. The primer pair of 528F and 706R that amplifies the 18S rRNA gene V4 hypervariable region was also assessed in silico, and was found to have a mismatch to S. singaporensis, an important ‘cluster’ of Sarcocystis in this study, hence disregarded. Amplifying hypervariable regions in the 18S gene for Sarcocystis detection was not a priority over the concerns of specificity. During the early phase of method exploration with the acquired samples, amplification of the 18S gene using published primers (Fischer & Odening 1998; Yang et al. 2001) yielded sequence fragments that were similar among multiple Sarcocystis species, according to a public database (data not shown).

The 98% coverage and 99% similarity of Uniq.V9.21756 to S. nesbitti HF544323.1 happened because the reference microorganism fell short of four nucleotides compared to the query sequence and has one nucleotide difference over the 156 bp of comparison. Importantly, Uniq.V9.21756 was identical to the parallel V9 region of WF2B, one of the in-house unknown Sarcocystidae species with nucleotide sequence ranging from the 18S to the 28S rRNA gene. This means Uniq.V9.21756, very likely, is WF2B. The corresponding region of WF2B in the 28S rRNA gene (R7–R8 region) was closely related to Sarcocystis mantioni (Zotu 42) and Zotu 50 (Figure 3, Supplementary Material S2). S. mantioni is the most prevalent sequence found in this study. In other words, notwithstanding the single nucleotide differences between WF2B and S. mantioni in the 28S R7–R8 region, and between Uniq.V9.21756 and S. nesbitti (HF544323.1) in the 18S V9 region, the most prevalent Sarcocystis species from the Tioman Island, is one that is highly similar to Sarcocystis nesbitti, as determined by NGS of two different identification genes (18S and 28S). This is the best possible postulation, as HF544323.1 is the longest publicly available sequence for S. nesbitti for the 18S gene. Sequences of other genes for this species have not been reported.

The cloning of sample I-2015 (using technical replicate I2-2015), which was a seawater sample from L1, produced three different sequences. One of it, D5.3, was identical to the sequencing results of primary PCR product used for cloning (results not shown), meaning that it is most likely the dominant species in the sample. The sequence of D5.3, which appeared twice among the clones, is reported as the last entry in Supplementary Material S2. D5.3 is 91% similar to Sarcocysits zamani KU244528.1, at 49% of query coverage. However, subsequent attempts to amplify sample I-2015 for NGS, and medium-range PCR targeting a longer PCR product with the same sample, were all not successful. This regrettably rendered the interesting quest for how much Sarcocystis species could survive in seawater, largely unanswered. Nevertheless, the initially successful PCR and the subsequent cloning results provided sufficient credibility to the reported sequence of C5.3 (467 bp) for future reference. This also forms the basis of the likelihood that Sarcocystis may survive in seawater, hence the possibility of contracting sarcocystosis through seawater recreational activities.

The use of NGS for quantitation has been cautioned due to replicate variability, whereby average Otu overlaps of only 17.2% and 8.2% were reported, respectively, among groups of two and three technical replicates (Zhou et al. 2011). Similarly, Zotus variability (28S R7–R8 region) is evident among the five technical replicates of the water filter, namely WF-1B_S5, WF-1A_S13, WF-2A_S14, WF-Sup_S21, and D.max.mar (Table 1, Supplementary Material S5). If replicate WF-Sup_S21 (supernatant) that was processed differently was left out, then the remaining four replicates shared 15.9% of Zotus. In comparison, the duplicate pair of J1-2015 and J2-2015 shared 71.1% of Zotus, and the pair of S1-2015 and S2-2015 shared 61.3% of Zotus. Referring back to the replicates of the water filter, the Chao1 (a measure of alpha diversity) values of WF-1B, WF-1A, WF-2A, and D.max.mar were 32.1, 17.5, 17.0 and 6.0, respectively (Table 1). The Chao1 values of J1-2015, J2-2015, S1-2015 and S2-2015 were 33.1, 31.5, 18.0 and 15.0, respectively. These data seem to imply that in the context of technical replicates, a higher percentage of Zotu overlaps might be expected if the alpha diversity metrics (measures of species richness in a sampling site) of the replicates are close in values. Following this simple presumptive ‘rule’, sample WF-1A and WF-2A that had close Chao1 values were found to be sharing 63.9% of their Zotus, as compared to the overall 15.9% among the four replicates. To circumvent the issue of technical replicate variability, methods such as removing singletons followed by sequence abundance-weighted OTU overlaps analysis were proposed (Wen et al. 2017). The 28S rRNA gene NGS results of this study suggest that before the inevitable sequence abundance-based ‘correction’ methods, the alpha diversity metrics among the technical replicates, such as Chao1 values, should be noted to give a preliminary idea as to whether a closely similar Zotu or Otu content can be expected. It is likely that technical replicates with very different alpha diversity values would have less shared Zotu/Otu/unique sequences, and hence, require deeper sequencing, if possible.

Hino et al. (2016) have previously shown that the V9 region of the 18S gene was able to detect various parasites from rat faeces. The microbial diversity of faeces is less compared to environmental samples. If pathogens are detected, their concentrations would most likely be significant, since they are already manifested clinically. PCR amplification workflow tested with clinical samples such as faeces, muscles or blood, is often assumed to be also useful for other applications. However, with environmental samples, this study, for example, has shown that the opposite could be true of the 18S V9 region. Similarly, PCR amplification of microalgae (dinoflagellates and diatoms) was reportedly prone to ‘contamination’ by fungus (Guo et al. 2016). The microbial diversity of environmental samples should, therefore, not be underestimated. In other words, PCR workflow optimised with clinical samples should not be assumed to be also useful for environmental samples – it needs to be tested with environmental samples itself.

This study presents a new region in the 28S rRNA gene for Sarcocystidae identification, which is better than the V9 region of the 18S rRNA gene. The Zotu sequences of the respective samples are not exhaustively discussed herein. Other than the nine named Zotus, some of the other interesting Zotus are Zotu 15, 3, 26, 40, 2, 29, 61, 50, and 59. Their actual identities, as well as the identities of all other Zotus reported in this study, await to be discovered, either as species that have been known using other identification genes like the 18S or ITS (internal transcribed spacer) as such, or as new species that would then be recognised as first detected on the Tioman Island of Malaysia. It is demonstrated here that NGS is a very useful tool in assessing outbreaks whereby a clearer picture of the microbiological entities contributing to the issue is needed, in order to devise effective outbreak control and prevention steps. Information from this study fills the critical environmental perspective gap of Sarcocystis research against the backdrop of various Sarcocystis species reported from animal hosts.

The author would like to thank the Director General of Health Malaysia for permission to publish this article. Dr Lu Ping Tan and Dr Zuraifah Asrah Mohamad from the Institute for Medical Research Malaysia are greatly appreciated for their review of the manuscript. Acknowledgement goes to personnel from The Tioman Island Sarcocystosis Investigation Team, the Environmental Health Research Centre and Molecular Pathology Unit of the Institute for Medical Research Malaysia, the Disease Control Division of the Ministry of Health Malaysia, and the National Zoo of Malaysia for their field, sample and laboratory assistance. This work was funded by the Ministry of Health Malaysia, with project code NIH/IMR/15-011 and registration number NMRR-15-2005-27199.