Schistosomes in the Persian Gulf: novel molecular data, host associations, and life-cycle elucidations

Avian schistosomes, comprise a diverse and widespread group of trematodes known for their surprising ability to switch into new hosts and habitats. Despite the considerable research attention on avian schistosomes as causatives of the human cercarial dermatitis, less it is known about the diversity, geographical range and host associations of the marine representatives. Our molecular analyses inferred from cox1 and 28S DNA sequence data revealed presence of two schistosome species, Ornithobilharzia canaliculata (Rudolphi, 1819) Odhner, 1912 and a putative new species of Austrobilharzia Johnston, 1917. Molecular elucidation of the life-cycle of O. canaliculata was achieved for the first time via matching novel and published sequence data from adult and larval stages. This is the first record of Ornithobilharzia from the Persian Gulf and globally the first record of this genus in a potamidid snail host. Our study provides: (i) new host and distribution records for major etiological agents of cercarial dermatitis and contributes important information on host-parasite relationships; (ii) highlights the importance of the molecular systematics in the assessment of schistosome diversity; and (iii) calls for further surveys to reach a better understanding of the schistosome diversity and patterns of relationships among them, host associations, transmission strategies and distribution coverage.


Results
Three out of the 1,745 examined P. cingulata were infected with avian schistosomes. The infected snails were collected at two distinct localities named Genaveh (n = 2; prevalence = 1%) and Jask (n = 1; prevalence = 0.4%) (see also Fig. 1 for sampling locations). Successful amplifications were achieved for 28S and cox1 for all three isolates. The yielded sequences were 1254-1285 bp (28S rDNA) and 344-730 bp long (cox1). The two isolates from Genaveh shared an identical 28S rDNA sequence with a published isolates for Ornithobilharzia canaliculata from the USA ex Larus delawarensis and L. occidentalis (AF167085, AY157248, KP734309), while the isolate from Jask differed by 2.3% (29 bp) from the former ones. A BLASTn search indicated that the latter isolate belonged to the genus Austrobilharzia. The novel isolate from Iran differed by 12-19 bp (0.9-1.8%) from the published representatives of the genus. The closest relative was an otherwise unidentified isolate from the same host species, P. cingulata, from off Kuwait (12 bp, 0.9% genetic difference).
Cox1 sequence divergence between our two isolates of O. canaliculata from Iran was 9 bp (2.6%). In contrast to the identical 28S sequences between the novel and published isolates for O. canaliculata, cox1 sequences differed substantially, ranging between 27 and 32 bp (7.9-9.3%). The single isolate from Jask differed by 1.66-1.82% (49-55 bp; 16.3-18.2%) from the novel isolates for O. canaliculata, and by 59 bp (20.1%) from Austrobilharizia sp. from Kuwait. Interspecific sequence divergence within Austrobilharzia was within the range of 21-63 bp (9.7-20.1%). However, the intergeneric divergence between the isolates for Ornithobilharzia and Austrobilharzia was somehow lower that the interspecific divergence for Austrobilharzia, i.e., 25-61 bp (7.7-18.4%). A single cox1 isolate for Austrobilharzia variglandis ex Larus sp. from Canada was not included in the sequence comparisons as it covers a distinct region of the cox1 gene and did not align with the remaining published isolates.
The aligned 28S dataset consisted of 76 terminals (2 newly-sequenced) and it was 1370 bp long, 78 of which were excluded prior to analyses. The cox1 dataset comprised 66 terminals and it was 1031 bp long. Analyses of the individual genes resulted in well-resolved trees (Fig. 2). The 28S rDNA hypothesis, presented in Fig. 2A, included representatives of all named and molecularly characterised species-level lineages except for the monotypic  under the GTR + I + Г model of sequence evolution. Analyses were run for 10,000,000 generation and 25% discarded as "burn-in". Posterior probability values are given above the branches; values. Nodes with < 0.95 posterior probability support have been collapsed. Branch length scale-bar indicates number of substitutions per site. Newly-generated sequences are indicated in colour indicated red and bold. Hosts of origin of individual sequences are indicated after the specimen's host name. Branches in blue indicate schistosomes with marine lifecycle. Shaded areas and taxa outlined with doted lines reflect on the current taxonomic framework of the family and also given on the right. www.nature.com/scientificreports/ Jilinobilharzia as molecular data currently do not exist (the single species, J. crecci Liu & Bai, 1976, has not been reported since its original description). Therefore, the ingroup taxa consisted of representative sequences of the families Schistosomatidae and the closely related Spirorchiidae (see Supplementary Table S1). The outgroup comprised representative of the Aporocotylidae and it was informed from previous phylogenies 25 . Our phylogenetic hypothesis recovered the spirorchiids in freshwater crocodilian and testudine hosts as the earliest diverging lineage. Spirorchiids with marine life-cycle clustered in a distinct clade basal to the Schistosomatidae. Members parasitic in marine testudines were identified as a distinct clade sister to all remaining schistosomes parasitic in birds and mammals. Schistosomes clustered into four distinct lineages: (i) an earlier diverging and strongly supported clade comprising the marine Ornithobilharzia and Austrobilharzia (ii) Macrobilharzia-a genus known from suliform birds which was resolved as a distinct lineage basal to the freshwater schistosomes, and two strongly-supported multi-taxa sister clades predominantly of (iii) mammalian and (iv) avian schistosomes. The mammalian schistosomes were further recovered as three distinct lineages: (i) Bivitelobilharzia-a genus including species parasitic in elephants and rhinoceros were recovered in a strongly-supported sub-clade sister to the main clade of mammalian schistosomes; (ii) a sub-clade of Schistosoma spp. with South East Asian distribution; and (iii) a clade comprising African representatives of Schistosoma. The North American mammalian representatives, Heterobilharzia and Schistosomatium were resolved as closer relatives to the large clade of avian schistosomes (Trichobilharzia + Marinabilharzia + Dendritobilharzia + Gigantobilharzia + Nasusbilharzia + Riverabilharzia). The remaining avian schistosomes clustered in two sister monophyletic clades with generally strong support for the major nodes. Bilharziella and Nasusbilharzia were recovered as earlier diverging to the two sister strongly-supported subclades of Gigantobilharzia + Dendritobilharzia + Marinabilharzia + Riverabilharzia, and Trichobilharzia + Allobilharzia + Anserobilharzia.
The cox1 tree was well-resolved and received strong support for most of the internal nodes (Fig. 2B). Taxa largely grouped in consistence with the 28S solution. Ornithobilharzia and Austrobbilharzia clustered into two distinct strongly-supported sister clades. The newly-sequenced isolate from Jask clustered in a clade with A. variglandis and A. terrigalensis; however, the isolate for O. canaliculata clustered with otherwise unidentified isolate labelled as Austrobilharzia sp. from Kuwait indicating a possible misidentification of the latter one. This was further confirmed by the high levels of genetic divergence in comparison with the other isolates of Austrobilharzia as indicated above.

Discussion
The present study is part of an effort to document the trematode diversity in P. cingulata (Gmelin, 1791), one of the most abundant snail species along the Iranian coast 23,24,26 . Sequence data for two species of marine avian schistosomes, Ornithobilharzia canaliculata (Rudolphi, 1819) and a putative new species of Austrobilharzia Johnston, 1917, are represented in a phylogenetic context together with other members of the family Schistosomatidae. This is the first report and molecular evidence for Ornithobilharzia canaliculata (Rudolphi, 1819) infecting P. cingulata as an intermediate host and it is the first partial molecular elucidation of its life-cycle. Our study adds to the diversity, host associations and phylogeny of the avian schistosomes with marine-based life-cycles, a group of schistosomes with great etiological importance.
Cercariae of the marine schistosomes are recognised as important etiological agents of human dermatitis [27][28][29] . Despite their importance to the public health, still very little is known about their diversity and evolution 7 as a consequence of largely under surveyed marine habitats for schistosomes worldwide 4 . This is in sharp contrast with the wealth of knowledge gathered about the mammalian and avian schistosomes with freshwater-based life-cycles, and information concerning the natural history of most marine schistosomes is scarce. The slow rates in recovering marine schistosomes, low species richness recorded in snail hosts and the convoluted taxonomy of the group, including separate taxonomic treatments of the distinct life-cycle stages, reflects the scarcity of data 30 . Matching sequence data for different life-cycle stages and across distant localities has accelerated life-cycles elucidations and host-parasite associations 4,5 . Although, the molecular systematics has had a major impact for the recent increase in discoveries and species delimitation, it has led to a plethora of putative new species and lineages of avian schistosomes for which only molecular data for their cercarial stages exist. Most of these putative species/species level lineages are of considerable importance due to their etiological significance. Their formal descriptions await as reliably identified adult stages are needed to help infer on their respective life-cycles and host-parasite associations.
Ornithobilharizia canalicata was originally described from Sterna galericulata in Brazil 11 . Later the species was reported from a wide range of gulls and terns across the Americas, Europe, Asia and New Zealand (see Table 1 for details). Larus dominicanus Lichtenstein and L. maculipennis Lichtenstein serve as the main hosts in the southern hemisphere; Larus delawarensis Ord, and L. occidentalis Audubon have been reported as hosts in North America and a total of 22 species of gulls and terns were reported as hosts across Europe and the Middle East. Despite the large number of definitive hosts, thus far the species was reported only from a single mollusc species, Lampanella minima (Gmelin), in North America. However, an experimental infection linking larval and adult stages has never been conducted. An important result from our study is the molecular confirmation of the conspecificity of our isolate from P. cingulata with the published isolate of an adult worm from North America. Matching sequence data for isolates from different life-cycle stages collected from disparate locations and times, provides unambiguous link between adult and larval stages from natural infections and accelerates species circumscription. The intercontinental distribution and the rather narrowly defined clade of gulls is instructive for studies on the transmission of avian zoonoses and the epidemiology of human cercarial dermatitis. The trans-continental distribution of Ornithobilharzia across the America, Europe and Asia is an explicit example that species dispersal is determined by the most vagile, bird host, involved in the trematode life-cycle. It is widely accepted that the distribution of the definitive host governs the larval trematode recruitment in the snail (first) Further, good documentation and re-evaluation of the morphological charters of the respective larval stages is urgently needed. Successful transmission of parasites with complex life-cycles requires an overlap of all hosts involved. The invertebrate first intermediate host has been recognised as one of the keys to the evolutionary expansions of the digenean trematodes. All schistosomes (marine and freshwater) are known to develop in gastropods. The basal position of the marine schistosomes (Austrobilharzia and Ornithobilharzia) has been considered as an indication for a successful ancestral marine-transmitted bird parasite transmission in colonising both freshwater snails and mammals 25 . The schistosomes emerging from marine heterobranch snails (Haminea and Siphonaria) and also recorded in penguins are a well-known example of secondary colonisation of marine habitats by the schistosomes 30,32,33 . Considering the snail intermediate hosts, in at least two instances, even congeneric schistosomes depend on markedly divergent gastropod lineages, i.e., pulmonates versus opisthobranchs or caenogastropods, indicative for an extensive host switching within the molluscan hosts 34 46 . Gyraulus has been reported as a natural host of the species in North America and New Zealand, while Anisus and Planorbis planorbis (L.) have been reported as hosts in Europe. Trichobilharzia jequitibaensis Leite, Costa, & Costa, 1978 is known to infect both lymnaeid and physid snails 47 . Austrobilharzia terrigalensis has been considered to utilise distinct snail hosts across its distributional range in Australia (Batillaria australis (Quoy & Gaimard)), North America (Cerithideopsis scalariformis (Say), and Ilyanassa obsoleta (Say) and the Pacific (Littorina pintado (W. Wood)). Intercontinental and trans-hemispheric distribution has been recently reported for Trichobilharzia querquedule 48,49 , however the species is known as a parasite specific to Physa spp. as an intermediate host.
In respect to their definitive hosts, a predominant part of the schistosomes is known as parasitises in birds. Currently a total of 13 genera are known as parasites in birds: Austrobilharzia Johnston, 1917 (6 species 5 and Trihobilharzia (c.35 species). Among them, species of Trihobilharzia have been subject of the most intensive research due to their recognition as leading etiological agents of human cercarial dermatitis ( 51 and references therein). The genus represents the most speciose among the bird schistosomes with about 35 species or species level lineages. However, about 65% of the remaining avian schistosomes, yet remain largely unstudied with fragmentary data on their diversity, biology and ecology. This is especially true in respect to the marine schistosomes for which life-cycle information lags behind their freshwater relatives. Despite the great efforts made so far in building up a comprehensive framework for the study of schistosome diversity, still considerable data are needed for assessing their true diversity due to the difficulty of directly relating larval and adult stages. Consistent efforts towards the use of integrative approach including collecting novel data from diverse host species and combining them thorough morphological examination, traditional systematics, classical taxonomy and phylogenetics have been proven as most valuable practice providing important information to the better understanding of the biodiversity and evolutionary relationships of the group 2-4 .
Our study strongly suggests that the biodiversity of the marine schistosomes is underestimated. The extensive discovery-based studies about schistosome diversity during the last two decades has revealed an immense diversity of avian schistosomes. However, unravelling their true diversity across hosts and geographic areas have been hindered by the difficulties of matching distinct life-cycle stages. There is still a need for more records of identifiable adult and larval schistosomes. Our study is a crucial step towards better understanding of important properties of marine schistosome biology and ecology, their patterns of diversification and distribution.

Materials and methods
Host and parasite collection. A total of 1745 adult P. cingulata (Gmelin, 1791) were sampled from 9 distinct locations along the Iranian coastline between December 2019 and February 2020 (Fig. 1). Samples comprised a minimum of 200 individual snails per locality, which were opportunistically collected by hand at the low tide from the intertidal zone. Snails were transferred alive to the laboratory, where they were measured (length and weigh) and labelled with a unique code given to each specimen. Each snail was then placed in an individual 50 ml beaker filled with filtered seawater and exposed to a warm light source for 3-4 h to simulate cercarial emergence. Beakers were screened under a stereomicroscope for the presence of cercariae indicating patent infections in the snail host. Prepatent infections were detected with snail dissections, which were conducted on the 3rd day of light stimulation. Both the released cercariae and schistosome sporocysts recovered from the host's tissue were washed with distilled water and preserved in molecular grade ethanol for DNA isolation and sequencing.
Sequence generation. Ethanol-preserved samples of pooled cercariae were subjected to DNA extraction and sequencing. Partial cox1 and 28S rDNA sequences were generated for the schistosome parasites recovered in order to achieve molecular identification and carry out reconstruction of their evolutionary relationships using published primers (28S (digl2 + 1500R 53,54 ; ECD2 + 900F 55,56 as internal sequencing primers; cox1: JB3 + JB4.5 57 or CO1-R 58 ). Contiguous sequences were aligned with MAFFT v.7 59,60 as an online execution. After alignment, sequences for cox1 were checked for stop codons using the echinoderm and flatworm mitochondrial code (translation table 9 61 ). All sequences were trimmed in order the first base to correspond to the first codon position in order to simplify position-coding in the downstream analyses.
Phylogenetic analyses. Phylogenetic analyses were performed on individual gene datasets using Bayesian inference (see Supplementary Table 1 for details on the taxa included in the analyses). Prior to analyses, the 'bestfitting' models of nucleotide substitution were estimated based on the Bayesian information criterion (BIC) in jModelTest v. 2.1.4 62 . BI analysis was carried out with MrBayes v. 3.2.7 63 on the CIPRES Science Gateway v.3.3 64 using Markov chain Monte Carlo (MCMC) searches on two simultaneous runs of four chains for 10 7 generations, sampling trees every 10 3 generations. The "burn-in" determined by stationarity of lnL assessed with Tracer v.1.5 65 was set for the first 25% of the trees sampled, and a consensus topology and nodal support estimated as posterior probability values 66