The Microbial Community of Tetrodotoxin-Bearing and Non-Tetrodotoxin-Bearing Ribbon Worms (Nemertea) from the Sea of Japan

A potent marine toxin, tetrodotoxin (TTX), found in a great variety of marine and some terrestrial species, leaves intriguing questions about its origin and distribution in marine ecosystems. TTX-producing bacteria were found in the cultivable microflora of many TTX-bearing hosts, thereby providing strong support for the hypothesis that the toxin is of bacterial origin in these species. However, metagenomic studies of TTX-bearing animals addressing the whole microbial composition and estimating the contribution of TTX-producing bacteria to the overall toxicity of the host were not conducted. The present study is the first to characterize and compare the 16S rRNA gene data obtained from four TTX-bearing and four non-TTX-bearing species of marine ribbon worms. The statistical analysis showed that different nemertean species harbor distinct bacterial communities, while members of the same species mostly share more similar microbiomes. The bacterial species historically associated with TTX production were found in all studied samples but predominated in TTX-bearing nemertean species. This suggests that deeper knowledge of the microbiome of TTX-bearing animals is a key to understanding the origin of TTX in marine ecosystems.


Introduction
Tetrodotoxin (TTX) is a low molecular weight non-proteinaceous neurotoxin found in a variety of marine and some terrestrial species [1]. TTX is usually detected together with its analogues (TTXs). It causes seafood poisoning in the Indo-Pacific region. TTX poisonings associated with invasive TTX-bearing species were also reported in Europe, America, and Oceania [2]. Despite numerous data on TTX occurrence in various species including TTX-producing bacterial strains [1] and phytoplankton [3,4], the origin of the toxin and its distribution in marine ecosystems remain unclear. The hypothesis of the bacterial origin of TTX that prevails today has several weaknesses; the most serious one is that the amount of TTX produced by bacteria is small compared to its content in the respective host [2]. At the same time, all available studies of TTX-producing microflora dealt with its cultivable forms that can be maintained under laboratory conditions, while the microbial communities of TTX-bearing animals were largely ignored. However, it is not possible to determine or at least estimate the proportion of bacterial species contributing to host intoxication. Studying the microbiome of TTX-bearing animals by using metagenomic analysis provides a new approach to detecting potential TTX sources. A group of animals comprising both TTX-bearing and non-TTX-bearing species would be a good model for such studies.
Nemerteans, or ribbon worms, are mostly carnivorous marine worms belonging to a phylum of their own. TTX is detected in members of all major nemertean groups. These animals use a variety of

Results
We sampled the microbiome of 24 ribbon worms representing four TTX-bearing and four non-TTX-bearing species (Table 1). Six species used in the current research, including Cephalothrix simula, Tubulanus punctatus, Kulikovia alborostrata, Quasitetrastemma stimpsoni, Cerebratulus cf. marginatus and Micrura cf. bella, collected from the Peter the Great Bay (Sea of Japan, Russian Far East) in July and August 2018 were tested for TTX presence or absence using high-performance liquid chromatography in combination with tandem mass spectrometry (HPLC-MS/MS) analysis in the previous study [13]. TTX analysis of Parahubrechtia sp. and Hubrechtella juliae were conducted in current research (see Section 2.4 in the Results).

Core Microbiome of the Nemertean Species Studied
The final OTU table generated by removing mitochondrial and chloroplast sequences (but prior to rarefaction) (Supplementary Materials) was used to search for conserved OTUs across all 24 nemertean samples. Only OTUs present in more than 50% (n = 51) of the samples were considered (Table 3). OTU_NEM_179 assigned to the genus Cupriavidus (Betaproteobacteria) was found in all nemertean samples. Other most conserved OTUs, namely OTU_NEM_358 (Comamonadaceae), OTU_NEM_285 (not assigned to any phylum), and OTU_NEM_132 (OD1 phylum), were found in all nemertean species (more than 87% of the samples). Other broadly conserved sequences included members of the classes Alphaproteobacteria and Gammaproteobacteria (Proteobacteria), Clostridia and Bacilli (Firmicutes), and Actinobacteria (Actinobacteria). The comparison of the core nemertean OTUs to other GenBank sequences showed that they were taxonomically closest to environmental samples, including ones from seawater and marine sediments. To evaluate differences in the conserved microflora between TTX-bearing and non-TTX-bearing nemerteans, a Venn diagram representing their core OTUs was generated ( Figure 2). A total of 13 OTUs found mainly in the microflora of toxic nemerteans and 10 OTUs characteristic to non-toxic nemerteans were revealed (Tables 4 and 5). The microflora typical of the toxic nemerteans was dominated by members of the orders Alteromonadales, Pseudomonadales, Oceanospirillales, and Vibrionales (Gammaproteobacteria). Within Gammaproteobacteria, characteristic to the non-TTX-bearing nemerteans, only Pseudomonadales and Oceanospirillales were revealed, while members of other Proteobacteria classes were also found. The core microbiome of the non-TTX-bearing nemerteans displayed a greater diversity of Firmicutes compared to the TTX-bearing ones.

Community Diversity of the Nemertean Microbiomes Studied
To statistically analyze the diversity and bacterial richness of the microbial communities studied, alpha and beta diversity metrics implemented in QIIME 2 were used. The average number of observed OTUs, Shannon index, and Faith's phylogenetic diversity index were used to determine the bacterial diversity of the nemertean species studied (Figure 3). Statistical testing showed no significant differences between the microbial communities of different nemertean species in the average number of observed OTUs (F = 1.57, p > 0.05, ANOVA; Figure 3a) and Faith's phylogenetic diversity index (F = 1.46, p > 0.05, ANOVA; Figure 3b), while the Shannon diversity index was significantly higher in Parahubrechtia sp. and Q. stimpsoni compared to the other species (F = 4.1, p = 0.009, ANOVA; Figure 3c). The Parahubrechtia sp. and Q. stimpsoni microbiomes did not differ significantly in the Shannon diversity index (p>0.05, Student's t-test). Thus, Parahubrechtia sp. and Q. stimpsoni had greater bacterial OTUs richness and evenness compared to the rest nemertean species. The differences in the microbiome diversity between the nemertean species studied was effectively demonstrated by principal coordinates analysis (PCoA) based on the Bray Curtis dissimilarity distance matrix (pseudo-F = 4.5, p = 0.001, PERMANOVA, 999 permutations in each test; Figure 4a). The samples from different hosts appeared to form at least three significantly different clusters, with the T. punctatus samples standing separately from the rest species, those from M. cf. bella and K. alborostrata forming another group and the other species grouping together. Interestingly, samples from the same species tended to cluster together except those from C. simula, which varied greatly in the relative abundance of OTUs. The same analysis applied to the cluster comprising C. cf. marginatus, H. juliae, Parahubrechtia sp., and Q. stimpsoni revealed a significant difference between these species (pseudo-F = 3, p = 0.001, PERMANOVA, 999 permutations in each test). PCoA based on Jaccard distance showed statistically significant differences between the samples from different nemerteans in the number of unique OTUs relative to the total number of OTUs (pseudo-F = 1.7, p = 0.001, PERMANOVA, 999 permutations in each test; Figure 4b). Only the samples from of C. cf. marginatus, Parahubrechtia sp., and K. alborostrata that formed a compact cluster were relatively similar in the membership of OTUs. The C. cf. marginatus and Parahubrechtia sp. samples were relatively close to two C. simula samples, assigned to the same cluster. The samples from the rest nemerteans were relatively different. To analyze OTU phylogeny, unweighted and weighted UniFrac distance matrices were used. PCoA based on either distance matrix revealed a significant difference between the nemertean species studied (Figure 5a; weighted UniFrac, pseudo-F = 3.8, p = 0.001, PERMANOVA, 999 permutations in each test; b. unweighted UniFrac, pseudo-F = 1.8, p = 0.001, PERMANOVA, 999 permutations in each test). By the relative abundance of observed microorganisms, based on the weighted UniFrac PCoA, the nemertean samples were clustered into two groups, with the K. alborostrata and M. cf. bella separated from the rest species (Figure 5a). By the relative occurrence of microorganisms, based on the unweighted UniFrac PCoA, significant differences were revealed both between different nemertean species and between samples from the same species (Figure 5b).

TTXs Analysis in Parahubrechtia sp. and Hubrechtella juliae
HPLC-MS/MS analysis of the Parahubrechtia sp. and H. juliae samples showed no chromatographic peaks characteristic of TTXs.

Discussion
To evaluate microbiomes, TTX-positive and TTX-negative nemerteans from different phylogenetic groups were chosen for the present study. HPLC-MS/MS analysis held in the previous research [13] revealed TTX presence in Cephalothrix simula, Tubulanus punctatus, Kulikovia alborostrata, Quasitetrastemma stimpsoni and absence of the toxin in Cerebratulus cf. marginatus and Micrura cf. bella. In the current study, no TTX was revealed in Parahubrechtia sp. and Hubrechtella juliae species.
The nemertean microflora has been poorly studied. In particular, the metagenomic approach was only used in one study characterizing two samples from toxic C. simula [12]. In the nemerteans presently studied as well as in C. simula [12], the same bacterial phyla, namely, Proteobacteria, Firmicutes, Bacteroidetes, and Actinobacteria, predominated. Some similarities were also observed in the composition of the most abundant orders; in particular, Rickettsiales, Alteromonadales, Pseudomonadales, and Vibrionales were common for all the hosts examined. On the genus level, Alteromonas, Pseudomonas, and Vibrio were relatively abundant both in C. simula [12] and the species presently studied. In the cultivable microflora associated with C. simula [22] and Q. stimpsoni, K. alborostrata, H. juliae, and Malacobdella grossa [20], members of the genera Pseudoalteromonas, Pseudomonas, Vibrio, Alteromonas, Shewanella, and Bacillus were found. The same genera were detected in the nemerteans presently examined. The predominance of Proteobacteria as well as Bacillus in these species is in accord with the data on the bacterial composition in marine ecosystems. This indicates that bacteria of these groups are part of the associative microflora of many marine invertebrates [23].
Though 457 OTUs were detected in the current research, only 13 core OTUs were shared among all nemerteans (Table 3). Such a small number of core OTUs may be accounted for by different feeding preferences of the nemerteans studied. The diet of ribbon worms includes semi-fluid, soft, or partly digested parts of arthropods, annelids (mostly oligochaetes and polychaetes), nematodes with a small proportion of living and dead mollusks, and fish [24,25]. The Cephalothrix species are also known to feed on other nemerteans and practice cannibalism [24]. The comparison of the core microbiome of the TTX-bearing and non-TTX-bearing nemerteans presently examined revealed more unique OTUs shared between the toxic species. Interestingly, members of the family Vibrionaceae, often associated with TTX production, were typical of the microflora of TTX-bearing nemerteans. Among them, Listonella anguillarum (=Vibrio anguillarum), a pathogenic bacteria causing hemorrhagic septicemia disease in marine and freshwater fish, bivalves, and crustaceans [26], was found in all toxin-bearing nemerteans, with its greatest abundance in Q. stimpsoni. Another common fish pathogen, occasionally affecting humans, Acinetobacter johnsonii [27] (Moraxellaceae), was found in almost all nemerteans except for C. simula and was most abundant in K. alborostrata.
The statistical analysis confirmed our previous conclusion that the nemerteans studied harbor distinct bacterial communities. Samples from the same species generally were more similar in bacterial composition compared to those from different species. The only exception was C. simula, whose three samples varied greatly in the microbial composition. Parahubrechtia sp. and Q. stimpsoni displayed significantly greater OTUs richness and evenness compared to the rest species. K. alborostrata and M. cf. bella were more similar in the relative abundance of observed OTUs as revealed by both Bray Curtis and weighted UniFrac PCoA analysis. Parahubrechtia sp., C. cf. marginatus, Q. stimpsoni, and H. juliae were assigned to the same cluster by Bray Curtis PCoA indicating their similarity in the microbiome composition. These results do not clearly demonstrate the relationship between hosts phylogeny and bacterial community. Microbial communities usually show more similarity between closely related hosts than between distant ones [28]. Among marine invertebrates, it was shown for sponges in a number of studies [29][30][31]. In the present study, K. alborostrata and M. cf. bella, both from the family Lineidae, had similar microbial communities. Parahubrechtia sp. and T. punctatus (Tubulanidae) on the one hand and C. cf. marginatus and K. alborostrata (Lineidae) on the other displayed sizeable dissimilarity in their microbiomes.
To date, members of 31 bacterial genera claimed in TTX production were isolated from various TTX-bearing organisms, as well as marine and freshwater environments [1]. Carroll et al. [18] were the first to report the relationship between the presence of Vibrio bacteria and TTX synthesis in seven species of the British nemerteans. It is worth noting that although the genus Vibrio, widespread in marine ecosystems, was commonly reported as the main bacterial TTX-producer, its ability to produce the toxin was questioned. Strand et al. [32] showed that using liquid chromatography-mass spectrometry (LC-MS) and bioassay for TTX detection in bacterial and nemertean samples might lead to false-positive results and should be supplemented with tandem mass spectrometry (MS/MS). The authors revealed a compound produced by Vibrio alginolyticus isolated from the ribbon worm Lineus longissimus. It was similar to TTX in molecular weight but differed from it in retention time and MS/MS fragmentation pattern and exhibited no toxic effect. However, recently, Turner et al. [12] confirmed the TTX production for V. alginolyticus isolated from the TTX-negative nemertean Tubulanus annulatus using LC-MS/MS. The present study showed that Vibrio is abundant in the nemertean microflora and is predominant in TTX-bearing species. V. alginolyticus was not found in the nemerteans presently examined. However, as was mentioned above, L. (=Vibrio) anguillarum, also claimed to produce TTX [33], was revealed in all TTX-bearing nemerteans. Interestingly, other members of Gammaproteobacteria, known to produce TTX, including Pseudomonas, Shewanella, Alteromonas, Marinomonas, and Pseudoalteromonas, were found in the majority of the nemertean species presently studied and were most abundant in TTX-bearing K. alborostrata and Q. stimpsoni. Recently, Turner et al. [12] isolated Pseudomonas luteola, capable of TTX synthesis, from C. simula. By using anti-TTX immunocytochemistry, Melnikova et al. [20] showed that the Pseudoalteromonas sp. strain 1942 isolated from H. juliae displayed TTX immunoreactivity. The bacterial strain Bacillus sp. 1839 isolated from C. simula was also reported to produce TTX first by using confocal laser scanning microscopy with polyclonal anti-TTX antibodies [19], and later on, by means of LC-MS/MS analysis [21]. In the present study, Bacillus species were more abundant in the microflora of non-TTX-bearing nemerteans, especially, in H. juliae and Parahubrechtia sp. and in TTX-bearing Q. stimpsoni. Moreover, the similarity search between a TTX-producing Bacillus sp. 1839 16s rRNA sequence and the present data revealed 16 sequences with 97% similarity with this strain in a Q. stimpsoni sample. We expected C. simula, the most toxic nemertean known to date, to display the greatest abundance of bacteria historically associated with TTX production. However, these bacteria occurred in all nemerteans presently studied with the predominance in the TTX-bearing species. In concordance with this, TTX-producing bacteria were found in TTX-negative nemerteans [12,21]. Nevertheless, one may conclude that TTX-positive nemerteans tend to accumulate more TTX-producing strains in their microflora compared to TTX-negative ones.

Sample Collection
Samples of Kulikovia alborostrata (Takakura, 1898), Micrura cf. bella (Stimpson, 1857), Cerebratulus cf. marginatus Renier, 1804, Cephalothrix simula (Iwata, 1952), Tubulanus punctatus Takakura, 1898, Quasitetrastemma stimpsoni (Chernyshev, 1992), Parahubrechtia sp., and Hubrechtella juliae (Chernyshev, 2003) were collected in July and August 2019 in the Spokoynaya (42.7090N, 133.1809E) and Vostok (42.9021N, 132.7385E) Bays of the Peter the Great Bay (Sea of Japan, Russian Far East). All nemertean specimens were kindly identified by Dr. Alexey V. Chernyshev, an expert in the nemertean zoology. Representatives of the majority of the species studied were collected in rhizoids of the brown algae Saccharina sp. at a depth of up to 1 m. Parahubrechtia and Hubrechtella were collected by trawling in mud at a depth of 5-7 m, and Cerebratulus, in sandy mud at a depth of up to 2 m. The water temperature in the collection sites was 17-21.7 • C. For 16s rRNA sequencing analyses, three individuals of each species were sampled. The samples were immediately transported to the laboratory, rinsed with sterilized seawater, placed in 100% ethanol, and stored at −20 • C until processing. For toxin testing, rinsed specimens were frozen at −20 • C without ethanol preservation.

DNA Extraction and 16S rRNA Sequence Processing
Ethanol-preserved nemertean samples of the Kulikovia, Micrura, Cerebratulus, Cephalothrix, and Tubulanus species were dissected into 20 mg cross sections used for DNA extraction. Specimens of the Parahubrechtia, Hubrechtella, and Quasitetrastemma species, taken for the same purpose, were not dissected due to their small size. Sections or whole worms were homogenized in a 1.5 mL Eppendorf LoBind tube containing 350 µL of lysis buffer with Proteinase K using a hand homogenizer for 3 min.
Using an E.Z.N.A. Mollusc DNA Kit (Omega Bio-tek, Norcross, GA, USA), total genomic DNA was extracted from each sample. The extraction was performed following the manufacturer's protocol with one exception: vortexing steps were replaced by gentle hand rotation of the samples. DNA purity and concentration were assessed using an UV5Nano spectrophotometer (Mettler Toledo, USA). The amount of DNA per sample varied from 39.7 to 509.8 ng/µL. All samples were stored at −20 • C until used.
For metagenomic analysis based on V3-V4 region 16S rRNA gene amplicons, total genomic DNA was amplified using Herculase II Fusion DNA Polymerase Nextera XT Index Kit V2 and sequenced on an Illumina Miseq platform with 2 × 301 bp paired-end reads (Macrogen Inc., Seoul, South Korea). The demultiplexing and quality control of the Illumina MiSeq dataset was processed by the Macrogen Inc.

Sequence Processing and Taxonomic Assignment
The quality of the reads was filtered using BBDuk software to remove Illumina adapters and artifacts (parameters for the pass: "ktrim = r k = 27 qtrim = rl trimq = 20 minoverlap = 24 minlength = 150"; http://jgi.doe.gov/data-and-tools/bb-tools/). Bioinformatics analysis was conducted using the Quantitative Insights into Microbial Ecology toolkit (QIIME) 2 version 2019.7 (https://qiime2.org) with proper built-in plugins [34]. The VSEARCH algorithm [35] was applied to merge the forward and reverse reads. Merged reads with an overlap length below 400 bp or over 570 bp were discarded. The quality-filtered sequences were dereplicated and combined into a single FASTA file for open-reference operational taxonomic unit (OTU) picking at 97% sequence identity. A representative sequence from each OTU was aligned against the Greengenes core set (database version 13.8) [36] using the feature-classifier plugin from QIIME2 with the classify-sklearn method [37,38]. The sequences that were not successfully classified by the Greengenes classifier were manually BLAST searched to specify taxonomy. Chimeric sequences were removed using the VSEARCH. Chloroplast and mitochondria sequences were excluded from downstream analyses. Core microbiomes were identified by the presence of an OTU in >50% (n = 51) of the nemertean samples.
Sequence files and metadata for all samples used in this study have been deposited in NCBI Sequence Read Archive repository under BioProject PRJNA604891 (accession numbers SRR11027358-SRR11027381).

Statistical Analysis
Alpha and beta diversity metrics were calculated using QIIME 2 [34] on a rarefied dataset to prevent statistical artifacts due to different sequencing depths. Alpha diversity was measured with the number of observed OTUs, Shannon diversity index, and Faith's phylogenetic diversity index. The sampling depth was set to 14590. The statistical significance of between-species differences in alpha diversity was evaluated using one-way analysis of variance (ANOVA) at p < 0.05. Beta diversity was measured with the Jaccard distance, Bray-Curtis distance, weighted and unweighted UniFrac distances, and clustered by PCoA. The statistical significance of sample clustering on these distance metrics was estimated using permutational multivariate analysis of variance (PERMANOVA) test at p < 0.05. Only OTUs found in at least 4 samples were included in the beta diversity analyses. The sampling depth was set to 7400.

Toxin Testing in Parahubrechtia sp. and Hubrechtella juliae
Extracts of Parahubrechtia sp. and H. juliae were prepared as described in Vlasenko et al., 2018 [13]. Briefly, samples were homogenized in 0.1% acetic acid in 70% methanol solution 1:10 (vol/vol) and ultrasonically treated. The homogenates were centrifuged, the supernatants were taken out and evaporated, and the resulting precipitates were dissolved in 0.1% aqueous acetic acid solution to the liquid-solid ratio of 1 mL/g, and filtered through centrifugal concentrators. The extracts were analyzed for tetrodotoxin (TTX) and its associated analogues (TTXs) by using HPLC-MS/MS following Bane et al. (2016) with modifications described in Vlasenko et al., 2018 [13]. TTX concentration was calculated using the calibration curve of standard TTX solution series (Alomone Labs Ltd, Jerusalem, Israel). The concentration of TTX analogues was calculated following the procedure proposed by Chen et al. [39]. The toxin detection criteria were the S/N ratio of the precursor MRM transition peak > 3, the relative intensity of the fragment ion peak > 4%, the order of toxins elution following Bane et al. [40]. The method was validated using standard TTX solutions in MRM mode. The linearity range varied from 0.6 to 100 ng/mL, the recovery range from 1 to 100 ng/mL of TTX was 98.4%, the LoQ was 0.6 ng/mL, the limit of detection (LoD) was 0.2 ng/mL, and the relative standard deviation varied within 4.5-14.6%.

Conclusions
The present study is the first to characterize the composition of the bacterial communities of TTX-bearing and non-TTX-bearing marine ribbon worms belonging to different taxa. The nemerteans studied proved to differ considerably in the microbiome composition. We also addressed the relationship between the presence of TTX-producing bacteria and TTX content in nemerteans, and we concluded that the TTX-bearing nemerteans presently studied tend to have more TTX-producing strains in their microflora. The present study provides the basis for further research on the microbial ecology and TTX distribution in Nemertea. Identifying TTX source(s) in one group of marine invertebrates may shed light on the TTX origin in marine ecosystems as a whole. Author Contributions: D.I.M. performed the toxin assays and DNA isolation from nemertean samples, analyzed the data, and wrote the manuscript. T.Y.M. conceived the study; performed sequence processing, their taxonomic assignment, and statistical analysis; and prepared figures. All authors contributed to the collecting of samples and editing of the manuscript and approved the final draft. All authors have read and agreed to the published version of the manuscript.
Funding: Not applicable.