Evolutionary history influences the microbiomes of a female symbiotic reproductive organ in cephalopods

ABSTRACT Many female squids and cuttlefishes have a symbiotic reproductive organ called the accessory nidamental gland (ANG) that hosts a bacterial consortium involved with egg defense against pathogens and fouling organisms. While the ANG is found in multiple cephalopod families, little is known about the global microbial diversity of these ANG bacterial symbionts. We used 16S rRNA gene community analysis to characterize the ANG microbiome from different cephalopod species and assess the relationship between host and symbiont phylogenies. The ANG microbiome of 11 species of cephalopods from four families (superorder: Decapodiformes) that span seven geographic locations was characterized. Bacteria of class Alphaproteobacteria, Gammaproteobacteria, and Flavobacteriia were found in all species, yet analysis of amplicon sequence variants by multiple distance metrics revealed a significant difference between ANG microbiomes of cephalopod families (weighted/unweighted UniFrac, Bray–Curtis, P = 0.001). Despite being collected from widely disparate geographic locations, members of the family Sepiolidae (bobtail squid) shared many bacterial taxa including (~50%) Opitutae (Verrucomicrobia) and Ruegeria (Alphaproteobacteria) species. Furthermore, we tested for phylosymbiosis and found a positive correlation between host phylogenetic distance and bacterial community dissimilarity (Mantel test r = 0.7). These data suggest that closely related sepiolids select for distinct symbionts from similar bacterial taxa. Overall, the ANGs of different cephalopod species harbor distinct microbiomes and thus offer a diverse symbiont community to explore antimicrobial activity and other functional roles in host fitness. IMPORTANCE Many aquatic organisms recruit microbial symbionts from the environment that provide a variety of functions, including defense from pathogens. Some female cephalopods (squids, bobtail squids, and cuttlefish) have a reproductive organ called the accessory nidamental gland (ANG) that contains a bacterial consortium that protects eggs from pathogens. Despite the wide distribution of these cephalopods, whether they share similar microbiomes is unknown. Here, we studied the microbial diversity of the ANG in 11 species of cephalopods distributed over a broad geographic range and representing 15–120 million years of host divergence. The ANG microbiomes shared some bacterial taxa, but each cephalopod species had unique symbiotic members. Additionally, analysis of host–symbiont phylogenies suggests that the evolutionary histories of the partners have been important in shaping the ANG microbiome. This study advances our knowledge of cephalopod–bacteria relationships and provides a foundation to explore defensive symbionts in other systems.

A ssociations with microorganisms have been crucial for the successful evolution of animals, often contributing essential roles in development and specialized functions (e.g., metabolism, digestion, protection from pathogens, and providing host camouflage) (1,2).Many cephalopods form symbiotic relationships with bacteria, for example, the specialized light organ of bobtail squid (Sepiolidae) and photophores in some loliginid squid host bioluminescent bacteria that are members of the Vibrionaceae and are used in a camouflage behavior known as counterillumination (3,4).Many cephalopods also have a second symbiotic organ called the accessory nidamental gland (ANG) that hosts a bacterial consortium and is part of the female reproductive system (5)(6)(7)(8)(9)(10).The ANG symbiosis has been studied in the Hawaiian bobtail squid Euprymna scolopes, where the ANG bacteria are recruited from the environment each generation (11)(12)(13).During oviposition, bacteria are deposited in the jelly coats of eggs, where they provide both antibacterial and antifungal protection for developing embryos (14)(15)(16).Antibacterial and antifungal activities have also been found in bacteria isolated from the ANG and eggs of the squids Doryteuthis pealeii (17) and Uroteuthis duvaucelii (18).
The basic structure of the ANG is similar across cephalopods.The organ consists of convoluted tubules composed of epithelial cells that surround extracellular bacteria (7,9).In E. scolopes, distinct tubules have different bacterial taxa (9).A few studies have described the development of the ANG in sepiolids (12,13), sepiids (19), and loliginids (20,21).In E. scolopes, the ANG begins to form around 4 weeks after hatching and reaches complete development at sexual maturity (approx. 2 months).It begins as a bi-lobed organ that is poised to recruit bacteria from the environment and eventually becomes a single uniform gland (12,13).
Although the ANG is distributed among multiple cephalopod species, whether the microbiome of this gland is conserved, as has been described for light organ/photophore symbioses (22,23), is unknown.Phylosymbiosis occurs when the host phylogeny reflects the host-associated microbiome community relationships at a moment in time and space and does not assume a stable evolutionary association (24,25).Though not universal, phylosymbiosis is observed in the microbiomes of many animals and insects, except, to date, for birds (26)(27)(28) and microscopic marine invertebrates (29).In marine organisms, phylosymbiosis has been observed in corals, sponges, ascidians (30), fishes (31) and marine mammals (32), and recently in the digestive tracts of some cephalo pods (33).In E. scolopes, complete ANG development depends upon the presence of environmental bacteria (12,13), and the symbiotic bacteria are vital in egg protection (14); however, whether the ANG exhibits phylosymbiosis is unknown.Examining the ANG microbiome of multiple cephalopod species for phylosymbiotic patterns may provide a first step to determining eco-evolutionary mechanisms that drive this pattern, as has been described in other host-microbe associations (25).
The ANG of E. scolopes has a stable microbiome of Alphaproteobacteria, Gammapro teobacteria, Verrucomicrobia, and Flavobacteriia (9,11).In other cephalopods, Alphapro teobacteria, Gammaproteobacteria, and Flavobacteriia were found in Doryteuthis pealeii of the Loliginidae family (7) and many Roseobacter spp.(Alphaproteobacteria) in Sepia officinalis of the Sepiidae family (6).Similar observations were made in 14 other cephalopod species belonging to families Loliginidae and Idiosepiidae (8).However, except for work published in E. scolopes, these previous studies relied on 16S rRNA gene clone libraries, fluorescence in situ hybridization, or characterizing a limited number of cultured strains.With high-throughput sequencing using 16S rRNA gene community analysis, we characterized the ANG bacterial diversity from four cephalopod families: Idiosepiidae (pygmy squid), Sepiidae (cuttlefish), Sepiolidae (bobtail squid), and Loliginidae (pencil squid), comprising 11 species and 7 genera to understand whether there is a core cephalopod ANG microbiome, if microbiota vary between different cephalopod groups and geographical regions, and whether phylosymbiotic signatures could be detected.This study lays the foundation for further investigations of mechanisms that drive the selection of ANG bacteria in cephalopods.

Sample collection
ANGs from 11 species of squid and cuttlefish were collected from multiple locations (Table 1; Fig. 1).All samples were collected in accordance with appropriate animal protocols and local permit collections: E. scolopes, protocol A18-029 approved by the Institutional Animal Care and Use Committee, University of Connecticut; Euprymna tasmanica: permit number 200891 U.C. Merced; and Idiosepius pygmaeus: Queensland Department of Agriculture Fisheries and Forestry, permit number 170251.Dissection of E. tasmanica followed animal ethics guidelines at James Cook University and was approved by the James Cook University Animal Ethics Committee (JCU Animal Ethics Number A2189; Sepia officinalis were collected by certificated fishermen at the Ría de Vigo, Spain, and in the Pertuis Charentais, France.The individuals were transported in proper containers to the Experimental Culture Facilities of IIM-CSIC and UMR LIENSs, which are registered as "User and breeding center on animal experimentation" (approval number ES360570202001 and 173002, respectively).Transport, housing, and handling were carried out following the principles of animal welfare, within 2 h after fishing.Procedures for transportation, euthanasia, and dissection were carried out in accordance with the principles published in the European Animal Directive (2010/63/EU) for the protection of experimental animals and were approved by the National Competent Authority ethics committee from Spain (Research Project ES360570202001/17/EDUC FORM 07/CGM01) and France (APAFIS#20520-2019050614554709).All animals were euthanized using anesthetic such as MgCl 2 (3%; wt/vol) and ethanol (1%; vol/vol) or ethanol (2%; vol/vol) dissolved in seawater, and ANGs were carefully dissected using sterilized scissors.
The ANG samples of E. scolopes, Euprymna morsei, Euprymna berryi, and D. pealeii from 2015 were flash-frozen upon collection.The ANGs from the remaining samples were dissected almost immediately or within 2 h from the time of collection and stored in DNA/RNA Shield (Zymo Research, Irvine, CA, USA) during transportation to the University of Connecticut and subsequently stored at −80°C until processing.The samples in Fig. 1 were mapped with QGIS (34).Two samples of Doryteuthis opalescens ANGs were also collected from Oregon, USA (44°40′01.2″N,124°17′02.4″W)(CDFW Scientific Collection Permit SC-13563), but due to the insufficient replicates were not included in the analyses.We have separately described these two replicates (Fig. S1).

Validation of preservative
To confirm the efficacy of DNA/RNA Shield as a preservation buffer, we compared the ANG microbiome from E. scolopes where one half of the ANG was stored at −80°C immediately after dissection and the other half was placed in DNA/RNA Shield at room temperature for 1 week (n = 3).Samples were processed, and the microbial composition and beta diversity of the halves were compared as described below.No significant differences in the ANG microbiome based on the preservation method were observed [Bray-Curtis permutational multivariate analysis of variance (PERMANOVA): q-value =0.82].

DNA extraction, sequencing, and analysis
DNA extractions were performed as previously described for E. scolopes (11).Briefly, ANG samples and five extraction controls (containing only filter sterilized artificial seawater) were transferred to 1.5-mL tubes and homogenized with a sterile pestle in ultra-pure water to dissociate the host tissue.Samples were centrifuged at 5,000 rcf to remove host cell debris, and the supernatant with bacterial cells was lysed with bead-beating and processed for DNA extraction with the Qiagen DNEasy Blood and Tissue kit.The DNA was quantified with the Qubit dsDNA high-sensitivity assay (ThermoFisher Scientific).Negligible DNA (<0.01ng/µm) was found in the controls, which were also processed for PCR and sequencing.
The V4 region of the bacterial 16S rRNA gene (primer 515f GTGYCAGCMGCCGCGGTAA and 806r GGACTACNVGGGTWTCTAAT) was amplified by PCR (35) and confirmed with agarose gel electrophoresis.No bands were detected in negative controls.PCR products were sequenced on a Miseq (Illumina Inc, USA) at Microbial Analysis, Resources and Services, University of Connecticut.The full amplification and sequencing protocol are available at http://www.earthmicrobiome.org/protocols-and-standards.
The raw sequences were denoised and demultiplexed using QIIME2.2020(36).The amplicon sequence variants (ASVs) were generated using the DADA2 method (37) and aligned against the trained Greengenes reference database at a 99% or higher identity (38) (May 2013 release).Some ASVs of Opitutae and Leisingera were not assigned by Greengenes and were manually annotated with NCBI BLAST (optimized for megablast).Chloroplast, mitochondria, and sequences unassigned at the kingdom level were filtered out.The number of reads in the negative controls ranged from 0 to 15 reads/sample.Among the negative controls, one ASV assigned as Rhodospirillaceae was common Colors correspond to different cephalopod families.

TABLE 1
The taxonomy and species names of the cephalopod samples collected along with their site of collection, number of replicates for each site, and the year the samples were collected between the negative controls (>50% of reads) and all replicates of Sepia officinalis (>40% of reads).However, we do not interpret the reads in S. officinalis as spurious since the number of reads far outweigh the negative controls by more than 3,000 reads.The other ASVs in negative controls were assigned as Bacilli and Rhodobacterales, but these specific ASVs were not found in any ANG samples.The QIIME2 artifacts were converted to phyloseq objects using qiime2R (v0.99.6).The final amplicon sequence variant (ASV)-abundance matrix was analyzed by two different methods: (1) normalized by transform_sample_counts and (2) rarefied by rarefy_even_depth, both to the minimum number of reads in the sample set (1,538) in R 3.3.2with the phyloseq (v1.26.1) package (39).Rarefaction curves were generated with alpha_rarefaction in QIIME2 (Fig. S2).No difference in alpha and beta diversity analyses was observed with the two methods; thus, the statistical analyses shown here are with transformed data.The data were analyzed for differences between families of cephalo pods and between species within the same family.We calculated alpha diversity with observed, Shannon, and inverse Simpson metrics using estimate_richness and plot_rich ness with phyloseq.Pairwise significances between comparisons were performed with the Wilcoxon test, and group comparisons were performed with the Kruskal-Wallis test with Bonferroni-corrected P-values.Beta diversities were analyzed with Bray-Curtis and weighted and unweighted UniFrac (40).The effect of host taxonomy on beta diversity differences was calculated with PERMANOVA with beta-group-significance (QIIME2) with 999 permutations.In addition to PERMANOVA in QIIME2, pairwise.adonis(v0.4) was used to analyze pairwise differences between groups.For diversity comparisons within families, samples were subset and normalized to their minimum number of reads (sepiolids = 7,671, loliginids = 6,289, and sepiids = 1,538).Other R packages used for making figures include ggplot (41), RColorBrewer (42), dplyr (43), ggpubr, scales, and reshape2 (44).

Comparison of Euprymna scolopes ANGs from present and past studies
ANG samples from E. scolopes in this study were combined with a previous study by reference (45).Pairwise comparison of samples from this study and reference (45) found no significant difference (Bray-Curtis PERMANOVA: P-value =0.3).Thus, regardless of the version of QIIME used [QIIME2.2020 in this study and QIIME1 in reference (45)], we observed no differences in microbiomes of E. scolopes.The ANG microbiome described by reference (45) using QIIME1 observed a higher abundance of Flavobacteriia (0.05%-10%).This extended analysis of additional ANG samples was not included in phylosym biosis analyses to ensure relatively similar sample sizes.The included E. scolopes samples were randomly selected from the larger available data set (11,45).Of those samples selected, re-analyses of ANG microbiomes from different geographic locations on Oahu, Hawaii, were not significantly different (Bray-Curtis PERMANOVA: P-value =0.879).

Phylogenetic similarity of Opitutae and Ruegeria sequences
Multiple sequences identified as Opitutae and Ruegeria were found in the sepiol ids.These sequences were submitted to NCBI with accession numbers (OQ305981-OQ305990).We used ARB SILVA's ACT (alignment, classification, and tree service) to determine the sources (host-associated or environmental) of closely related sequences (95% similarity cut-off) to sequences from sepiolid ANGs.Pelagiococcus litoralis was used as the outgroup for the Opitutae tree and Leisingera cerula as the outgroup in the Ruegeria tree.FastTree with default parameters (GTR+CAT model) was used to generate the tree in the Newick format.

Phylosymbiosis
Phylosymbiosis was analyzed using methods as previously described (46).ASV abundan ces were first averaged for the replicates, and microbial cladograms were constructed using Bray-Curtis and weighted UniFrac distances to test for patterns of phylosymbiosis (26).The host tree phylogeny was obtained from Tanner et al. (47), modified to retain only the species used in this project using drop.tipwith ape package (48).The host trees were modified to include sepiolid phylogenies resolved by Sanchez et al. (49).We considered ASVs present in at least 50% of the replicates and with at least 30% relative abundance in each host species as dominant and were identified with QIIME2 plugin q2-coremicrobiome.The topological congruences of the host and microbiome trees were measured with normalized matching cluster and Robinson-Foulds metrics with 10,000 random trees using a previously published script (26).Both metrics produce a score ranging from 0 (completely congruent) to 1 (completely incongruent).Mantel tests (vegan package) with 1,000 permutations were run to evaluate the Pearson and Spearman correlation of host divergence times and microbiome distance and host geographic distance and microbiome distance.The host divergence times were obtained from references (47,49).The geographic distance matrix was generated using the Geographic distance matrix generator (v1.2.3).The distance decay plots were plotted with betapart (v1.5.4) package.

ANG bacterial communities differ between cephalopod families
The total number of amplicon reads from all 70 samples was 2,153,636, ranging from 1,538 to 123,544 reads per sample, with an average read depth of 22,879.Rarefac tion curves, plotting observed ASVs against the sequencing depth from 1,538 reads, displayed a plateau for most species, indicating that the sequencing depth was adequate to capture the bacterial diversity (Fig. S2).The normalized data set had 1,221 unique ASVs.Significant differences were observed with Wilcoxon tests in alpha diversity measures between families of cephalopods (Fig. 2A).There was no significant difference in richness and evenness between species of loliginids and sepiids, but within sepiolids, E. scolopes had a significantly higher diversity (Kruskal-Wallis rank sum test, Bonferroni-corrected P > 0.05, Fig. S3).ANG communities from the loliginids had the highest richness while those from the sepiolids had the lowest (Fig. 2A).
PERMANOVA of beta diversity metrics (weighted UniFrac, unweighted UniFrac, and Bray-Curtis) revealed that the overall microbiomes of ANGs between families of cephalopods were significantly different (Fig. 2B and C; Table 2).All comparisons of species/families with unweighted UniFrac were significantly different (q < 0.01).Ordination plots of principal component analyses with weighted UniFrac and Bray-Curtis distances also showed close clustering of samples within the same family (Fig. 2B  through E).There was an overlap between the ANG microbiomes of I. pygmaeus (sole representative of Idiosepiidae) and some sepiolids.Pairwise PERMANOVA with weighted UniFrac of I. pygmaeus with each sepiolid revealed similarity to the ANG microbiomes of E. tasmanica and E. morsei, but not other sepiolids (Table 2).Significant differences were observed with weighted UniFrac and Bray-Curtis distances within loliginids and sepiids (Table 2).A smaller significant difference (q-value =0.03) within sepiolids was observed with weighted UniFrac and a larger difference with Bray-Curtis (q-value =0.001,Table 2).Pairwise comparison within sepiolid species revealed a significant difference between all species using Bray-Curtis (q-value =0.001), in weighted UniFrac (q-value =0.032), and between combinations of Eumandya parva, Euprymna berryi, and E. scolopes (Table S1).Other pairwise comparisons of sepiolid species were not significant using weighted UniFrac (q-value >0.11,Table S1).These analyses indicate that while most cephalopod species harbored distinct ANG communities, some of the more closely related sepiolids harbored similar ANG microbiomes.
The year of collection and location of the same species did not affect the microbiome composition.For example, no difference in the ANG microbiome was observed between D. pealeii collected in 2015 and 2019 (Table 2).Also, the most abundant ASVs (Rhodobac teraceae, Alphaproteobacteria) remained consistently present in all D. pealeii samples regardless of the time of collection (Fig. 4).There was also no significant difference with weighted UniFrac in the ANG microbiomes of S. officinalis collected from coasts of France and Spain, separated by 748 km (Table 2).However, S. officinalis samples did vary by collection site when analyzed via unweighted UniFrac (q-value =0.002, R 2 = 37%) and Bray-Curtis (q-value =0.05,Table 2), and alpha diversity was also significantly higher in S. officinalis collected from France than those collected from Spain (P = 0.01, Fig. S3).
Two samples of Doryteuthis opalescens ANG (analyzed separately due to the low sample size) had a high abundance of Alphaproteobacteria (22% ± 3%) and Gammapro teobacteria (60% ± 15%) but lacked other taxa like Bacteroidia and Flavobacteriia that were found in D. pealeii ANGs (Fig. S1).One ASV identified as an Alphaproteobacteria (accession number: OQ435377) was common to all the Doryteuthis samples (Fig. S1).This shared ASV was further identified with the NCBI database as a Hyphomicrobiales of genus Anderseniella.
For each cephalopod species, 5-10 ASVs made up ~50% of the relative abundance of their microbiomes (Fig. 4).Although some ASVs were shared between species, no single ASV was found common across all samples.For example, an Alphaproteobacteria ASV was shared between the loliginids D. pealeii and U. duvaucelii.Also, a Ruegeria ASV was shared between E. morsei and E. parva, Illumatobacter sp. was found in majority of E. berryi, E. tasmanica, and E. parva, and an Opitutae ASV was found in most E. morsei, E. berryi, and E. tasmanica samples.Overall, closely related bacteria were found between sepiolids.However, similar observations were not made with different Sepia species (S. esculenta and S. officinalis), which had significantly different microbiomes (Fig. 4; Table 2).
The alignment of Opitutae and Ruegeria ASVs found in sepiolids, I. pygmaeus, and S. esculenta with other closely related bacteria from the ARB databases (SILVA) (50) revealed that these ASVs were associated with both hosts and environmental samples (Fig. 5).The Opitutae ASVs from sepiolids were more closely related to each other than other sequences in the ARB databases (Fig. 5).Sepiolid Opitutae ASVs were also related to uncultured Verrucomicrobia from the ANG of the European squid Loligo vulgaris (Fig. 5).The ASVs from this study had a >99% sequence identity with ANG bacteria from E. scolopes from previous studies (Fig. 5) (9).Phylogenetic trees of the sepiolid Ruegeria ASVs classified them as closely related to different species such as R. lacuscaerulensis and R. atlantica (Fig. 5).These varying results may suggest that the Opitutae in different sepiolids are more closely related to each other or, more likely, that there are insufficient cultured representatives available for Opitutae sequences compared to those for Ruegeria spp.

DISCUSSION
In this study, we surveyed the microbiomes of a female-specific symbiotic organ, the ANG, in multiple host species from four cephalopod families.Overall, the cephalopod ANG microbiomes have a lower bacterial diversity compared to the mammalian gut and oral microbiomes (52)(53)(54), marine sponges (55), and corals (30).High abundances of Alphaproteobacteria were found in all species of cephalopods similar to the observations made in previous studies of single host species (7,8,11) while lower abundances of Gammaproteobacteria were also found in all host species.However, at lower bacterial taxonomic levels, the ANG bacteria varied between host families and between host species.We also provide further evidence that ANGs from mature adults were consistent over time and space.The ANG microbiomes of D. pealeii were similar from samples collected in 2015 and 2019, and S. officinalis from waters of France and Spain were similar (Table 2).These findings agree with observations made for E. scolopes (11,45) where ANG communities were consistent between collection years and from different sampling locations on Oahu, HI, USA.Overall, these combined studies suggest a consistent species-specific ANG microbiome.A consistent microbiome across time and space may imply a strong selection pressure to acquire specific ANG bacteria during colonization.
The alpha diversity of ANG symbionts between cephalopod families varied significantly.The loliginids have a richer and more even bacterial community compared to the sepiids and sepiolids (Fig. 2).This difference may be due to many reasons, such as the range of host migration, the necessity for increased protection from diverse pathogens, or other factors.The loliginids have a large migratory range, so hosts may have increased exposure to different microbial niches in different environments (56)(57)(58).Beta diversity analyses showed significant differences in microbiomes between families (Table 2).Bray-Curtis and weighted and unweighted UniFrac analyses suggest that microbiomes between sepiids, loliginids, sepiolids, and I. pygmaeus were different.Further, micro biomes within sepiids and loliginids were also different, which may suggest unique microbiomes between species in these families.Though the microbiomes of the ANGs from closely related sepiolids were more similar and clustered together with weigh ted UniFrac, they were different with the Bray-Curtis analysis (Table S1; Fig. 2).Since weighted UniFrac considers taxonomic relationships of bacterial communities and Bray-Curtis does not (59), similar microbiomes within sepiolids as analyzed with weighted UniFrac suggest that related but not identical bacterial communities are present in the ANGs of this group.
Most of the sepiolids possessed a relatively high abundance of Ruegeria (Alphaproteo bacteria in the Roseobacter clade) and Opitutae (Verrucomicrobia).Species of Ruegeria and Opitutae made up ~50% of the relative abundance of most sepiolid ANG microbiomes (Fig. 3 and 4).The presence of these taxa in the ANG of E. scolopes was also confirmed using fluorescence in situ hybridization (9).A previous study found that the relative abundance of Opitutae had a negative correlation with squid size in E. scolopes, while the relative abundance of Alphaproteobacteria had a positive correlation to squid size, suggesting that members of the Opitutae may initially colonize the ANG in this species (12).Follow-up work also suggests a role for these bacteria in initiating the development of the ANG in E. scolopes (13).The high Opitutae prevalence in all the sepiolids and I. pygmaeus suggests a conserved function for this bacterial family in these cephalopod groups (Fig. 3 and 5).To date, no member of this group has been successfully cultured, but the data from this study suggest that the Opitutae may play an important role in the ANG symbiosis of some cephalopods.
Cephalopods evolved ~540 Mya in the Cambrian, and species used in this study diverged ~120-15 Mya (47).With this diverse data set, we found that recently diverged cephalopod species had similar microbiomes compared to those that diverged more than 30 Mya (Fig. 6; Table S1).The sepiids diverged approximately 120 Mya, followed by idiosepiids (~95 Mya), sepiolids, and loliginids (~80 Mya).Within the cephalopod families, the sepiolids analyzed for this study diverged from each other ~15-30 Mya (Table S1), the sepiids S. officinalis and S. esculenta ~50 Mya, and the loliginids ~80-100 Mya (47,49).The Mantel test and the dendrogram-based approaches are thought to have high sensitivity in detecting the host phylogenetic correlation, with the Mantel test being more sensitive (59).The Mantel test here indicated an increase in host divergence times correlated with an increase in microbiome divergence.The differences observed in the microbiomes of species within the loliginid and sepiid families are likely due to the increased divergence times (Fig. 6A and C).A similar pattern was found in Nasonia sp.wasps, where host species that diverged 0.4 Mya had a more similar microbiome than those that diverged 1.0 Mya (60).Differences between the sepiid ANG microbiomes may also be explained by variable egg-laying behaviors.Sepia officinalis deposits ink into their egg cases during laying, whereas S. esculenta does not (61).As cephalopod ink is known to have antimicrobial properties (62,63), functional differences in egg protection may also explain the selection of different microbiomes found in these two cuttlefish species.Phylosymbiosis was also observed in the gut microbiomes of other cephalopod species, but these putative symbionts and their function are largely unknown (33).
Phylosymbiosis of multiple organs within the same host has only been examined for a few animals, e.g., the gut and/or skin of mammals (64,65), bats (52), and fishes (31).However, the strength of phylosymbiosis varies depending on the organ type and host exposure to the environment (59).Some cephalopods also provide two different internal symbiotic organs (the light organ and the ANG) with specialized roles.The light organ or photophore in many cephalopods houses bioluminescent bacteria belonging to the Vibrionaceae (Gammaproteobacteria) in sepiolids and loliginids (4) that provide camouflage via counterillumination against nocturnal predators (66,67).Like the ANG, bacteria in light organs also have phylogenetic congruence with their squid hosts (67).Colonization experiments of different strains of Vibrio fischeri revealed that native or closely related strains are more effective at colonizing the light organ (67,68), and the requirements for light organ colonization are well studied in E. scolopes (23).Thus, bacterial adaptation may also have a role in phylosymbiosis.Further experiments such as the transplantation of ANG bacteria between species of cephalopods may reveal the functional benefits of phylosymbiosis.Recent studies that describe the development of the ANG in E. scolopes suggest that such gnotobiotic experiments under controlled conditions may be possible to test these hypotheses (13).The potential to examine phylosymbiosis between these light organs and ANGs, which in contrast to the gut microbiome are posited to have only a restricted open period for colonization, could provide a unique insight into varying selective pressures moving forward.
Host-microbe relationships that are conserved over evolutionary time may indicate a significant functional role (25).There are many reasons that phylosymbiosis can arise such as coevolution, cospeciation, ecological drift, or shifts in the host geographic range (25).Though phylosymbiosis is observed in many species, its functional significance for associations is still unclear (25).Transfer experiments of microbiota between Nasonia wasp species demonstrated an impact on larval growth and adult survival, which suggests that selective pressure underpins phylosymbiosis in the case of Nasonia spp.(24,46).While both the ANG and light organ symbioses are defensive, the selective pressures associated with the loss of function of these symbioses are difficult to test in the laboratory.Because the light organ and ANG are not nutritionally obligate symbioses for the host, the lack of an organ does not affect host survivability.However, in animals raised under aposymbiotic and symbiotic conditions, bacterial colonization is necessary for morphogenesis and/or full development of the organs (3,13,23).The ANG and its bacteria may have a similar relationship like that of the light organ and vibrio species in multiple cephalopods as both these relationships are important for host defense.
Defensive symbioses are common in both terrestrial and aquatic environments (69).Bacterial symbionts, metabolites, and/or biosynthetic gene clusters responsible for antimicrobial activity have been identified from the ANG and/or eggs of E. scolopes (14-16), D. pealeii (17), and U. duvaucelii (18).Compounds and strains of Alphaand Gammaproteobacteria and Flavobacteriia isolated from E. scolopes demonstrated various degrees of inhibition against three strains of the fungus Fusarium keratoplasticum, the yeast Candida albicans, and other marine bacteria (14)(15)(16).The fungal-inhibiting species included Ruegeria sp., which here was found in all sepiolids, I. pygmaeus, and S. officinalis.Leisingera sp.isolated from E. scolopes eggs produces the antibacterial compound indigoidine ( 16), but species from this genus were not found in any other sepiolids.Pseudoalteromonas spp.from E. scolopes also strongly inhibited fungi and other bacteria (14,15) and here were also detected in I. pygmaeus, E. morsei, E. parva, and D. pealeii.Members of Actinobacteria are known for their potent secondary metabolites (70) and were found here in multiple cephalopod ANGs including E. berryi, E. tasmanica, E. parva, and the sepiids.Future studies will explore the antimicrobial potential of ANG symbionts from other cephalopods and potentially culture new and uncharacterized strains.
We also evaluated the hypothesis that geography may contribute to differences observed in ANG microbiomes between host families.While we do not have many multiple representatives of the same species collected from different geographical locations, S. officinalis collected from France and Spain were separated by 730 km, yet their ANG microbiomes were very similar, even at the ASV level.Many cuttlefish, including S. officinalis, exhibit a genetic pattern of isolation-by-distance, where the populations have increased genetic differences with increased geographies (<300 km), i.e., they have low mobility (71,72).This could mean that populations of S. officinalis from France were different from those in Spain.Previous research also found that two genetically isolated populations of E. scolopes had very similar ANG bacterial composi tions (45).These findings also support the hypothesis that geography may not significantly influence the selection of bacteria at the broader taxonomic level.In sepiolids, species or strains of bioluminescent vibrios in the light organ correlated to the geo graphical distribution of the hosts (73).Similarly, the ANG microbiome data suggest that related but distinct microbiomes are selected by closely related cephalopods.The species of Ruegeria and Opitutae were found in the ANGs of sepiolids, and the ASV differences may be attributed to differences in geography or strain availability during ANG colonization.The similar microbiome in sepiolids from different geographic regions separated by 1,000-8,000 km also suggests that the influence of geography may not be significant at the broader taxonomic levels but may impact microbiome composi tion at the species/strain level.E. morsei and E. berryi were collected from the same location, but these two species had distinct (Bray-Curtis analyses) but similar (weighted UniFrac analyses) microbiomes (Table S1).If ANG bacteria in E. morsei and E. berryi are recruited during a discreet colonization window as has been described for E. scolopes (12,13), then geographic differences during host development may account for varying bacterial diversity.Future studies raising different sepiolids in the same substrate and water conditions may help answer this question.Given the role of ANG bacteria in egg protection, bacteria with similar antimicrobial function may be selected by the host.Functional redundancy in the infant human gut for Bifidobacterium longum involved in milk metabolism has been observed (74).Furthermore, like the similarities observed for the microbiome species in the mammalian gut (53,74), host factors shared by different cephalopod species may help select similar symbionts.Future studies that compare the full-length 16S gene along with metagenomes may reveal strain diversity and functional redundancy in ANG microbiomes between cephalopod species.
Cephalopods obtain their ANG bacteria via horizontal transmission each genera tion; thus, environmental bacterial communities are important to consider (11)(12)(13)21).Although we did not collect environmental seawater samples in this study, previous work showed that the majority of operational taxonomic units (OTUs) found in the ANG of E. scolopes were also detected in the seawater and sediment of the host's natural environ ment (11).Similar microbiomes within S. officinalis collected from different locations and within sepiolids may suggest a strong selective pressure for certain bacteria.Life histories may also influence how different cephalopods obtain their ANG bacteria.
Cephalopod families also have different life history patterns.Some sepiolids and sepiids bury themselves in sand, whereas loliginids live a pelagic lifestyle (61,75,76).However, sepiolids and sepiids do not have similar microbiomes, so lifestyle is unlikely to be the sole contributor to bacterial selection.
The 16S rRNA microbiome analyses of the ANG microbiome from 11 species of cephalopods revealed shared bacterial communities at higher taxonomic levels.Finer analyses demonstrated distinct microbiomes associated with cephalopod families and species.The high bacterial diversity associated with the ANG defensive symbiosis may provide a unique opportunity to screen for potentially novel antimicrobial compounds and other secondary metabolites.We observed differences in the microbiomes as the host evolutionary distance increased.Sepiolids collected from geographically distant regions separated by thousands of kilometers had similar microbiomes, suggesting a conserved functional role of the ANG bacterial community and an evolutionary relationship between hosts and symbionts.Moreover, a growing number of cephalopod genomics and symbiont metagenomic/transcriptomic resources (77-79) will shed light on the evolution of the ANG and its functional role in hosts.

FIG 1
FIG 1 Representative photos of the ANGs from different families of cephalopods (left) and sampling locations (right).The ANG is situated inside the white inset box.Scale bars: Sepiidae 1 cm; Idiosepiidae 0.1 cm; Sepiolidae 0.5 cm; Loliginidae 0.5 cm.Right: sampling locations and species from which ANGs were collected.

FIG 3
FIG 3 Percent relative abundance of taxa with ASVs greater than 0.1% of normalized data set.Representative figures of cephalopods are above the names of the cephalopod families (where Idioindicates Idiosepiidae)."n" represents the number of biological replicates from each species.

FIG 4
FIG4 Bubble plot of the top 44 most abundant ASVs that made up approximately 40%-50% of each cephalopod species.The ASVs were identified to the closest taxonomic level assigned.Similar colors of the bubbles represent ASVs from the same class/phylum.The size of the bubble reflects the relative abundance of the ASV in each sample.The bar plot at the bottom corresponds to the percent relative abundance of the ASVs in the bubble plot for each sample.

FIG 5
FIG 5 Phylogenetic tree of abundant Opitutae spp.and Ruegeria spp. in ANG microbiomes.Maximum likelihood phylogenetic tree from FastTree (Newick format) of 240 bp of 16S rRNA gene sequences identified as (A) bacterial class Opitutae spp. of phylum Verrucomicrobia in sepiolids, Sepia esculenta, and Idiosepius pygmaeus and (B) genus Ruegeria spp. in sepiolids and I. pygmaeus.The symbols represent the sources of the sequences.

FIG 6
FIG 6 Comparison of the dendrograms of hosts and associated ANG microbiomes.The microbiome trees were generated with (A) weighted UniFrac distances and (C) Bray-Curtis distances.The topological congruence of the microbiome and host trees were calculated with nMC.(B, D) The rate of microbiome divergence across phylogenetic distance (B: weighted UniFrac, D: Bray-Curtis) with divergence times in Mya.The pink line represents distance decay.The P-values and Mantel r statistics are shown in the bottom right corners of the plot.

TABLE 2
Comparison of ANG microbiomes with permutational analysis of variance of weighted UniFrac and Bray-Curtis distances a