Host Phylogeny Determines the Gut Microbial Landscape of Cephalopods

Background: Compared to vertebrate gut microbiomes, little is known about the factors shaping the gut microbiomes in invertebrates, especially in non-insect invertebrates. Class Cephalopoda is the only group in the phylum Mollusca characterized by a closed circulatory system and a well-differentiated digestive system to process their carnivorous diet. Despite their key phylogenetic position for comparative studies as well as their ecological and commercial importances, analyses of the cephalopod gut microbiome are limited. In this study, we characterized the gut microbiota of six species of wild cephalopods by Illumina MiSeq sequencing of 16S rRNA gene amplicons. Results: Each cephalopod gut consisted of a distinct consortium of microbes. Photobacterium and Mycoplasma were prevalent in all cephalopod hosts and were identied as core taxa. The gut microbial composition reected host phylogeny. The importance of host phylogeny was supported by a detailed oligotype-level analysis of operational taxonomic units assigned to Photobacterium and Mycoplasma, although Photobacterium typically inhabited multiple hosts, whereas Mycoplasma tended to show host-specic colonization. Further, we showed that class Cephalopoda has a distinct gut microbial community from those of other molluscan groups. The gut microbiota of the phylum Mollusca was determined by host phylogeny, diet, and environment (aquatic vs. terrestrial). Conclusion: We provide the rst comparative analysis of cephalopod and mollusk gut microbial communities. The gut microbial community of cephalopods is composed of the distinctive microbes and strongly associated with their phylogeny. The genera Photobacterium and Mycoplasma are core taxa in the cephalopod gut microbiota. Collectively, our ndings of this study provide evidence that cephalopod and mollusk gut microbiomes reect phylogeny, environment, and the diet of the host and these data can be suggested to establish future directions for invertebrate gut microbiome research.

Cephalopoda is one of the oldest and most successful groups in the phylum Mollusca [14]. Over 11,000 species of cephalopods have been documented, including about 800 extant species. They are found in all of the oceans and at most depths, from the surface to the deep sea [15]. Cephalopods are the only group of mollusks with a closed circulatory system, analogous to that of vertebrates [16], and they have the most advanced nervous system among invertebrates [17]. Furthermore, cephalopods have a welldifferentiated digestive system to process their carnivorous diet [18]. Thus, cephalopods are expected to have a unique gut microbiome that differs from those of other mollusks. Furthermore, they are major component in the marine food chain as predators of zooplankton, crustaceans, and small shes and as prey of vertebrate predators [19]. Cephalopods are an important shery resource for many countries [20].
Despite these unique biological properties and the ecological and commercial importance of cephalopods, few studies have evaluated the cephalopod gut microbiome. S Iehata, F Valenzuela and C Riquelme [21] described the Octopus mimus gut microbiome using a 16S rDNA clone library and Á Roura, SR Doyle, M Nande and JM Strugnell [22] characterized the gut microbiomes of wild and captive Octopus minor paralarvae. Notably, studies focused on the adult cephalopod gut microbiota using a highthroughput sequencing approach are lacking, thereby limiting our understanding of the co-evolution of animals and the gut microbiome.
Here, we characterized the gut microbiota of six species of wild cephalopods, cuttle sh (Sepia esculenta), beka squid (Loliolus beka), inshore squid (Uroteuthis edulis), Japanese ying squid (Todarodes paci cus), common octopus (Octopus vulgaris), and whiparm octopus (Octopus variabilis), by 16S rRNA gene sequencing using the Illumina MiSeq platform. We investigated associations between host phylogeny and the gut microbiome and compared our data with previously reported data for mollusks and other aquatic animals to investigate the factors that shape the gut microbial composition in the phylum Mollusca.

Sampling
Cuttle sh, beka squid, inshore squid, Japanese ying squid, common octopus, and whiparm octopus were captured from offshore waters in the Republic of Korea. All samples were transferred to the laboratory directly in the living state and sacri ced in the laboratory by an adequate anesthetic method.
The dorsal mantle length and weight of each individual were determined. Samples were then dissected to obtain the stomach, cecum, and other digestive organs. In Supplementary Table S1, detailed metadata for the cephalopod samples are presented.
Identi cation of cephalopod hosts by Cytochrome oxidase I sequencing Cephalopod subjects were initially subjected to basic taxonomic identi cation based on morphological characteristics. For the detailed identi cation of cephalopods, genomic DNAs were extracted from the esh of specimens aseptically. A fragment of each tissue sample was suspended in 750 ml of lysis buffer and homogenized by FastPrep-24 (MP Biomedicals, Santa Ana, CA, USA) with glass beads (0.5 mm diameter) for 45 s at 5.0 m/s. Standard phenol-chloroform DNA extraction was preformed after lysis. The extracted DNAs were PCR-ampli ed using cytochrome c oxidase subunit I (CO1) primers designed for diverse metazoan invertebrates [67]. PCR products were puri ed using the QIAquick PCR Puri cation Kit (Qiagen, Hilden, Germany) following the standard protocol and bidirectionally sequenced using an automated DNA analyzer system (PRISM 3730XL DNA Analyzer; Applied Biosystems, Foster City, CA, USA) and the BigDye Terminator Cycle Sequencing Ready Reaction Kit (Applied Biosystems).
The sequence fragments were assembled using SeqMan (DNASTAR). The assembled CO1 gene sequences were compared with other CO1 gene sequences in the nucleotide collection (nr/nt) in GenBank by a BLAST search (Additional le 2: Supplementary DNA extraction and sequencing of bacterial 16S rRNA genes From dissected intestinal tracts of collected cephalopods, the cecum was primarily used to investigate gut microbial communities of cephalopods. The luminal contents of dissected cecal samples were also collected and used to extract microbial DNAs. To maximize microbial cell lysis for DNA extraction, the cecum and luminal contents were homogenized by shaking in a sterile screw tube containing zirconia beads (2.3 mm, 0.1 mm diameter) and glass beads (0.5 mm diameter) using FastPrep-24 (MP Biomedical) for 50 s. After lysis, genomic DNAs from the homogenized gut samples were extracted using the Qiagen DNA Stool Mini Kit (Qiagen). The V3-4 hypervariable region of the 16S rRNA gene was ampli ed with the primers 341F and 805R, and four independently ampli ed products for each sample were pooled and puri ed using the QIAquick PCR Puri cation Kit (Qiagen) to minimize bias. Libraries were prepared using the Nextera XT DNA Library Preparation Kit for Illumina MiSeq platform (Illumina, San Diego, CA, USA). The prepared DNA libraries were sequenced by certi ed service provider (Macrogen, Seoul, Korea) using the Illumina MiSeq platform with 2 × 300 bp reads, according to the manufacturer's instructions.

Sequence analysis
Raw 16S rRNA sequence data were processed using QIIME 1.9.1 [72]. Paired-end sequence reads were assembled with default parameters and minimally quality ltered with a Phred quality score threshold of 20. Data were then error-ltered using a de novo chimera removal algorithm, USEARCH (v. 7.0.1090) [73] High-quality sequence reads were assigned to OTUs by an open-reference OTU picking protocol [74] using the QIIME toolkit, where the UCLUST [73], OTU picking algorithm was applied to search sequences against the Greengenes reference database at a 97% sequence similarity threshold [75]. A representative sequence for each OTU was aligned with the Greengenes reference using PyNAST [76]. For bacterial taxonomic assignment, Ribosomal Database Project (RDP) classi er (version 2.3; https://rdp.cme.msu.edu/classi er/classi er.jsp) was used with 80% as a con dence value threshold [77]. An even-depth rare ed OTU table matrix (6000 sequences) was constructed to calculate various αand β-diversity indices and for further microbial analyses using QIIME pipelines. Sequences belonged to Mycoplasma and Photobacterium, were re-clustered with minimum entropy decomposition (MED) for sensitive discrimination of closely related organisms [26].

Network-based analysis of Mycoplasma and Photobacterium
Network maps of Mycoplasma and Photobacterium were generated using QIIME and visualized using Cytoscape (version 3.4.0) [78,79]. Brie y, the even-depth rare ed MED tables constructed with Mycoplasma and Photobacterium and converted to Cytoscape format using a QIIME script (make_otu_network.py). In the converted MED network maps, samples and MEDs represented nodes of the network and these nodes were connected by edges, indicating the abundance of the MED in the samples. Edge-weighted spring embedded models were derived for network arrangement.

Comparison of gut microbiomes of cephalopods and various animal
Sequence data for the sea slug (Elysia chlorotica) and eastern oyster (Crassostrea virginica) gut microbiome were obtained from the MG-RAST server (mgp561 and mgp1994, respectively; http://metagenomics.anl.gov) [35,36]. Sequence data for Hawaiian land snail (Auriculella ambusta) gut microbiome was downloaded from NCBI Sequence Read Archive (SRP047488, https://www.ncbi.nlm.nih.gov/sra) [34]. Sequence data for the blood clam and marine sh gut microbiome were used from our unpublish data. Since the targeted region and applied sequencing technologies differed among experiments, we assigned taxonomic characteristics against the identical reference database using RDP classi er. After unaligned sequences were discarded, the even-depth rare ed OTU table was generated and used for further analyses. Non-phylogenetic distance metrics (binary Jaccard and Bray-Curtis dissimilarities) were calculated and visualized by a principal coordinates analysis (2D PCoA).

Statistical analysis
All statistical analyses were performed using GraphPad Prism (v. 8.43; GraphPad, San Diego, CA, USA). Two-tailed Mann-Whitney U-tests were used to compare the β-diversity indices between multiple groups. Analysis of similarities (ANOSIM) and multivariate ANOVA based on similarities (adonis) tests with the βdiversity matrix were performed using the QIIME pipeline (compare_categories.py)[80, 81]. Statistical signi cance for both tests was determined based on 10,000 permutations.

Results
Characteristics of the cephalopod gut microbiota After sequence quality ltering (and excluding sequences found fewer than 15 times in the entire sample), 3,661,327 high-quality reads from 30 samples were generated with a mean sample depth of 122,044 and a standard deviation of 20,693 reads.
After rarefaction, 76,381 high-quality sequences were clustered into 1,835 OTUs at a 97% sequence identity threshold (357 ± 103 OTUs per sample). The phylogenetic diversity index, an alpha diversity measure, was used to estimate bacterial species richness (Additional le 1: Supplementary Fig. S1). The Chao1 metric reached a plateau after 75000 reads, suggesting that the depth of coverage was su cient to capture nearly all of the biological diversity within samples (Additional le 1: Supplementary Fig. S2).
Cephalopod gut microbial community re ects host phylogeny Taxonomic pro les clearly show that cephalopod gut microbiotas share core taxa but also include unique components (Fig. 1). A beta diversity analysis supported this observation. The cephalopod gut microbial communities were clustered according to host species in a PCoA of unweighted UniFrac distances (Fig. 1b). Furthermore, average interspeci c variation was signi cantly higher than intraspeci c variation in the beta diversity matrix (Fig. 1c). Similarity of the gut microbial composition was higher within the order Octopoda compared to within the orders Sepiida and Teuthoidea. A PCoA based on UniFrac unweighted distances indicated that species belonging to the order Octopoda were similar, with signi cantly lower intra-order beta diversity than those for other orders (Fig. 1b and d).
In vertebrates, phylogeny and host diet in uence the structure of gut microbial communities [23,24]. We evaluated whether there is a similar correlation between host phylogeny and gut microbial composition in cephalopods. We constructed a host phylogenetic tree based on mitochondrial housekeeping genes reported by JE Uribe and R Zardoya [25]. We also generated an unweighted-pair-group method with arithmetic-mean (UPGMA) tree based on 16S rDNA gene sequences from the gut microbial communities (Fig. 2a). The host phylogenetic tree and the UPGMA tree of the gut microbiomes showed identical topologies, indicating that the gut microbial communities in cephalopods are closely related to the host phylogeny. As shown in Fig. 2b, a heat map indicated that host CO1 similarity and microbial dissimilarity are strongly correlated with cephalopod host species. Furthermore, the order Octopoda showed high intraorder CO1 and microbial composition similarity.
A majority of OTUs were matched to the genera Mycoplasma and Photobacterium, which were regarded as the core taxa of the cephalopod gut microbiota (48.3% and 23.8%, respectively; Additional le 1: Supplementary Fig. S3). Although OTUs belonging to these genera were differentially distributed according to host phylogeny, the limited taxonomic resolution rendered an OTU-level analysis ineffective.
Furthermore, sequences included in major OTUs are overestimated during taxonomic strati cation, distorting the sequence distribution. To overcome these obstacles, we decomposed the OTUs assigned to identical genera (Mycoplasma and Photobacterium) and re-clustered the sequences into ne-scale units using nucleotide entropy by the minimum entropy decomposition (MED) method, an unsupervised oligotyping approach [26]. The OTUs belonging to Mycoplasma and Photobacterium were resolved into 228 oligotypes. Distortion in the sequence distribution was reduced for oligotypes (Additional le 1: Supplementary Fig. S4). We performed a network analysis using oligotypes to evaluate the distribution of the core taxa with better taxonomic resolution (Fig. 3). The distributions of oligotypes among hosts were consistent with the aforementioned results for the core OTUs and showed host-speci c connections. In the case of Mycoplasma, oligotypes were divided into three sub-clusters according to host: cuttle sh and Japanese ying squid; beka squid and inshore squid; whiparm octopus and common octopus. The majority of Photobacterium oligotype nodes were connected to multiple hosts. There was also a striking difference in co-speciation patterns between the genera Mycoplasma and Photobacterium in the oligotype-level phylogenetic analysis (Additional le 1: Supplementary Fig. S5). In Mycoplasma, we found that most oligotypes colonized a single host species. Oligotypes assigned to Photobacterium with earlier diverged, were found in multiple host species, whereas those that diverged more recently were hostspeci c.
Sexual dimorphism is not associated with the gut microbial community in cephalopods Male and female cephalopods have different growth rates, body sizes [27], and organ structures [28,29]. Moreover, there is evidence that sexual differences affect the diet and ecological niche [30,31] of cephalopods. Because the host diet and habitat are critical determinants of the microbiota [32,33], we investigated whether the cephalopod gut microbiota is shaped by sex (Additional le 1: Supplementary   Fig. S6a). Sex did not affect the gut microbial community at the levels of class (Additional le 1: Supplementary Fig. S6b) and order, except in the case of the order Teuthida (Additional le 1: Supplementary Fig. S6c).
Host phylogeny, diet, and habitat shape the gut microbiota of mollusks We further compared the gut microbiota of cephalopods and other mollusks and identi ed the relative contributions of various environmental and genetic factors to the microbial community composition. We obtained data for the gut microbiomes of four mollusk species from public databases, including the Hawaiian land snail (Achatinella mustelina) [34] and emerald sea slug (Elysia chlorotica) [35] belonging to the class Gastropoda and oyster (Crassostrea virginica) [36] and blood clam (Anadara broughtonii) (our unpublished data) belonging to the class Bivalvia. A marine sh gut microbiome (62 species, our unpublished data) was also included in the analysis for comparison between mollusks and vertebrates.
Each mollusk class showed a largely distinct gut microbial composition (Fig. 4a). The core taxa in cephalopods, Mycoplasma and Photobacterium, were not predominant in other mollusks or shes. Moreover, a PCoA based on binary Jaccard distances (Fig. 4b) showed that cephalopods form a single cluster distinct from vertebrates and other molluscan groups, although each species in Cephalopoda had a unique microbiome. The classes Gastropoda and Bivalvia each formed a single cluster as well. Emerald sea slug was most similar to cephalopods (Fig. 4b). The class Gastropoda, including the emerald sea slug, was the most recently diverged branch within the class Cephalopoda (Fig. 4c). These results suggested that the molluscan gut microbiota re ects host phylogeny. Interestingly, the Hawaiian land snail belonging to class Gastropoda had a completely different gut microbiota from those of other marine mollusks, despite their phylogenetic similarity ( Fig. 4d and Additional le 1: Supplementary Fig.  S7a). We then compared the gut microbiotas of marine sh to those of mollusks to determine the role of the marine environment. The gut microbiotas of marine mollusks were more similar to those of marine shes than to those of the Hawaiian land snail, despite their evolutionary divergence ( Fig. 4d and Additional le 1: Supplementary Fig. S7a). This These results indicate that gut microbiota of marine animal is potentially affected by the aquatic environment. Among marine mollusks, diet was also a strong determinant of the gut microbiota. The cephalopod gut microbiota could be categorized according to the host diet in a binary Jaccard analysis (Additional le 1: Supplementary Fig. S7b).

Discussion
Studies of invertebrate microbiomes have largely focused on insects, although Mollusca is the second largest invertebrate phylum [11]. In particular, few studies have focused on cephalopod gut microbiomes [21,22], with a lack of studies of multiple species and analyses of the factors determining the community structure. Thus, we characterized the gut microbiota of six wild cephalopod species (cuttle sh, beka squid, inshore squid, Japanese ying squid, common octopus, and whiparm octopus).
Based on our comparative analysis of 16S rRNA gene sequences obtained by Illumina MiSeq sequencing, we identi ed the genera Mycoplasma and Photobacterium as the core taxa in the cephalopod gut microbiota. These genera are the predominant taxa in the digestive tracts of wild Chilean octopus [21] and aquacultured common octopus [22]. Mycoplasma is an obligate parasitic bacterial group and is a component of the gut microbiome of many marine animals, such as Norway lobster [37], jelly sh [38], and Atlantic salmon [39]. Their roles in the intestinal ecosystem are typically recognized as pathogenic or opportunistic bacteria in vertebrates [40][41][42]; however, little is known about their roles in invertebrates, other than a report of a potential symbiotic Mycoplasma in scorpion [43]. The genus Photobacterium is well known for its bioluminescence [44] and pathogenicity [45,46]; however, their phylogeny and taxonomy are not clearly elucidated [47]. Members of the genus Photobacterium show ecological diversity, including taxa that are symbiotic [48][49][50] or parasitic [51,52] with sea animals, free-living in seawater [53] and in saline lake water [54], and even piezophilic [55]. Light production is a common feature of many genera in Vibrionaceae, and Photobacterium is one of the most extensively studied groups [56,57]. In this study, the genus Photobacterium was particularly abundant in beka squid (58.0%) and inshore squid (75.9%), members of the suborder Myopsida, containing Hawaiian bobtail squid (Euprymna scolopes). The Hawaiian bobtail squid is famous for its light-associated symbiosis and symbiont-speci c immune tolerance with the bioluminescent bacterium Aliivibrio scheri [58,59], once assigned to the genus Photobacterium [47]. Although beka squid and inshore squid do not have bioluminescence, the dominance of Photobacterium in Myopsida hosts suggests that there is a general symbiotic relationship between Myopsida hosts and Vibrionaceae bacteria.
Our results revealed that host phylogeny is re ected in the gut microbiota of cephalopods, indicating that the host phylogeny can be predicted by the gut microbial community. Phylogenetic analyses based on many housekeeping genes in Cephalopoda have yielded contradictory results, making evolutionary relationships within the class di cult to de ne [60, 61]. The cephalopods included in our study have clearly resolved phylogenetic positions based on well-supported consensus trees [25,60,62]. Our cephalopod gut microbial composition-based UPGMA tree matched the consensus tree for the host species. Accordingly, the gut microbial composition is a potential target for studies of cephalopod phylogeny.
In microbial community analyses by 16S amplicon sequencing, sequences are typically clustered into OTUs based on similarity, with a typical threshold of 97%. This clustering process is bene cial for downstream analyses; however, with respect to the operational de nition of a species, 3% dissimilarity is only a rough approximation. There is a risk that closely related species could be identi ed as a taxonomic unit in the clustering process. Furthermore, OTU-based analyses show a limited resolution for analyses below the genus level. The MED method overcomes a number of the limitations of the OTU-based approach. MED provides a computationally e cient means to partition marker gene datasets into MED nodes, which represent homogeneous OTUs. We used the MED approach to perform a network analysis at the within-genus level. The oligotyping analysis revealed different co-evolutionary histories between two major cephalopod species. The distribution of oligotypes of Mycoplasma was concentrated with host-speci c colonization; however, a large number of Photobacterium oligotypes were located in multiple cephalopod species. Based on these results, Mycoplasma colonization in cephalopods was frequently related to host-speci c evolution or biological activities, while Photobacterium colonized cephalopods more broadly, and interactions with Photobacterium might be essential for survival or adaptation of cephalopod species. This nding agrees with the Atlantic cod gut microbiome study [49].
We expected the cephalopod gut microbiota to differ between sexes based on differences in the growth rate, body size, diet, and space niche between male and female octopuses. However, we did not detect a signi cant difference in gut microbial composition according to sex. This lack of a difference has several potential explanations. First, we used cephalopod samples that are similar in size and were collected simultaneously at the same location, thereby minimizing the effects of body size, space niche, and similar variables. Second, there are con icting results regarding the difference in octopus diet between sexes [63, 64]. Therefore, a meta-analysis is needed to clarify the differences between male and female cephalopods. This is the rst comparative analysis of the cephalopod and mollusk gut microbiota. We identi ed three factors that in uence the gut microbiota of cephalopods and mollusks: host phylogeny, habitat type, and diet. The host phylogeny was the most prominent determinant of the gut microbiota. Although all cephalopod hosts in our study had similar diets and living environments, gut microbial compositions were distinguished by host phylogeny. This was supported by our beta-diversity, phylogenetic, and network analyses. With respect to the living environment, the gut microbiota of aquatic mollusks was more similar to that of sh than to that of terrestrial mollusks, suggesting that environmental conditions overwhelm other factors. However, further research is still needed, including analyses of mollusk samples in a wider range of environments. Lastly, the mollusk gut microbiota was distinguished by diet in our beta-diversity analysis. However, diet was a host phylogeny-dependent factor in our study. Therefore, to assess the independent effect of diet on the gut microbiota of mollusks, follow-up studies are needed.
We found that the features of cephalopod and mollusk gut microbial communities were quite similar to the common features of the vertebrate gut microbiota, which is also affected by host phylogeny [7], evolutionary divergence time [65], living environment [5], and diet [66]. The shared characteristics of the microbiomes suggest that insights from studies of the vertebrate gut microbiota can be applied to invertebrate studies; this can help establish directions for invertebrate gut microbiome research. In addition, new ndings based on invertebrate gut microbiome studies have the potential to be applied to vertebrate and human research.

Conclusions
Taken together, we performed the rst comparative analysis of the cephalopod gut microbiota by a highthroughput sequencing approach. We revealed that each species in class Cephalopoda has a unique gut microbiota. We identi ed the genera Mycoplasma and Photobacterium as core taxa in the gut microbiota of cephalopods. Furthermore, we found that the cephalopod gut microbial community composition was determined by host phylogeny, and host phylogeny is also an important determinant of the gut microbiota of marine mollusks. Diet and habitat also contributed to the mollusk gut microbiota.

Consent for publication
Not applicable.

Availability of data and material
The newly generated 16S rRNA sequence datasets are available in the European Nucleotide Archive (ENA) of EMBL-EBI under the accession number PRJEB27490. The cytochrome oxidase subunit 1 (CO1) gene sequences used for identifying host species have been submitted to NCBI GenBank (https://www.ncbi.nlm.nih.gov/genbank) under accession numbers MH542436-MH542464 (under the title "Factors shaping invertebrates gut microbiota: host phylogeny, habitat, and diet are involved in shaping of gut microbiota of Cephalopoda, Mollusca").