Global 16S rRNA diversity of provannid snail endosymbionts from Indo-Pacific deep-sea hydrothermal vents.

Symbioses between invertebrate animals and chemosynthetic bacteria build the foundation of deep-sea hydrothermal ecosystems worldwide. Despite the importance of these symbioses for ecosystem functioning, the diversity of symbionts within and between host organisms and geographic regions is still poorly understood. In this study we used 16S rRNA amplicon sequencing to determine the diversity of gill endosymbionts in provannid snails of the genera Alviniconcha and Ifremeria, which are key species at deep-sea hydrothermal vents in the Indo-Pacific Ocean. Our analysis of 761 snail samples across the distributional range of these species confirms previous findings that symbiont lineages are strongly partitioned by host species and broad-scale geography. Less structuring was observed within geographic regions, probably due to insufficient strain resolution of the 16S rRNA gene. Symbiont richness in individual hosts appeared to be unrelated to host size, suggesting that provannid snails might acquire their symbionts only during a permissive time window in early developmental stages in contrast to other vent molluscs that obtain their symbionts throughout their lifetime. Despite the extent of our dataset, symbiont accumulation curves did not reach saturation, highlighting the need for increased sampling efforts to uncover the full diversity of symbionts within these and other hydrothermal vent species.

Introduction

Microbial symbioses are increasingly recognized as universal phenomena that impact virtually all levels of biological organization, from cellular to organismal to ecosystem scale (Bronstein, 2015). Growing evidence from various symbiotic partnerships suggests that microbial symbioses can expand the physiological and ecological capabilities of hosts and symbionts, which are predicted to be critical for ecosystem productivity, stability and biogeochemical cycling (Apprill, 2017; Beinart, 2019; Wilkins et al., 2019). Deep‐sea hydrothermal vents are probably some of the most enigmatic ecosystems that are sustained by microbial symbioses. In these systems, invertebrate animals live in association with chemoautotrophic bacteria that use chemical energy from venting fluids for the production of organic carbon, thereby providing food for their hosts (Dubilier et al., 2008; Sogin et al., 2021). Despite decades of research on this topic and the significance of chemosynthetic symbioses for ecosystem processes at hydrothermal vents, the diversity and distribution of symbionts within and across hosts and habitats remain underexplored, especially at large biogeographic scales.

Provannid snails of the sister genera Alviniconcha and Ifremeria are dominant animals in benthic communities at deep‐sea hydrothermal vents in the Indian and Western Pacific Ocean (Van Dover et al., 2001; Desbruyères et al., 2006). While the Western Pacific genus Ifremeria is represented by a single species, I. nautilei, that affiliates with methane‐ and/or sulfide‐oxidizing gammaproteobacterial symbionts (Windoffer and Giere, 1997; Borowski et al., 2002; Suzuki et al., 2006a), the genus Alviniconcha comprises five Western Pacific species (A. adamantis, A. boucheti, A. hessleri, A. kojimai, A. strummeri) and one Indian Ocean species (A. marisindica) that live in symbiosis with thiotrophic Gammaproteobacteria or Campylobacteria (Suzuki et al., 2006b; Johnson et al., 2015; Breusing et al., 2020). In both Alviniconcha and Ifremeria, the symbionts are assumed to be horizontally acquired and are harboured intracellularly within the host's gill tissue (Suzuki et al., 2006a, 2006b). Despite an environmental pathway for symbiont transmission, host and symbiont genera or species appear to exhibit a relatively strong selectivity in their partnerships towards each other (Beinart et al., 2012; Breusing et al., 2020), though host individuals are flexible in recruiting local strains of their specific symbiont phylotype(s) (Breusing et al., 2021a; Breusing et al., 2021b).

Most current analyses on the variation and structure of microbial symbionts within Alviniconcha and Ifremeria stem from studies in the Lau Back‐Arc Basin, while little is known about these patterns in populations from other spreading systems within the distributional range of these genera. Here, we compiled an extensive dataset of 761 snail samples from 10 geographic regions of the Indo‐Pacific Ocean (Fig. 1), some of which were previously unexplored, to assess the global diversity of chemosynthetic gill endosymbionts within Alviniconcha and Ifremeria through identification of 16S rRNA amplicon sequence variants (ASVs). Using ordination analyses and correlative statistics, we determined the influence of host species, host size, depth and geography on symbiont composition and distribution.

Fig. 1

Locations for Alviniconcha and Ifremeria species sampled in this study.

Results and discussion

Symbiont 16S rRNA diversity is partitioned by host species and geography

Our conservative analysis pipeline, which extends a previous study by Breusing et al. (2020) to now include seven species and 10 geographic areas, recovered 60 symbiont ASVs that were assigned to two campylobacterial (Sulfurovum, Sulfurimonas) and four gammaproteobacterial (Ca. Thiobios, Methylomonas, Thiolapillus, unclassified Thiomicrospiraceae) genera of provannid snail endosymbionts (Figs 2 and 3). Average pairwise identities within genera ranged from 95% to 99% (Sulfurovum: 95.4%; Sulfurimonas: 95.0%; Ca. Thiobios: 97.1%; Methylomonas: n.a.; Thiolapillus: 98.1%; unclassified Thiomicrospiraceae: 99.0%). In agreement with Breusing et al. (2020), ASVs were generally segregated by host species and broader geographic region (i.e. back‐arc basin, volcanic arc or mid‐ocean ridge), except for lineages within the unclassified Thiomicrospiraceae group which were shared between A. kojimai and A. strummeri (Fig. 2, 4; Appendix 1: Fig. S1). Based on PERMANOVAs and linear decomposition models the impact of host species and geography superseded the influence of DNA preservation, extraction and sequencing method (81.17% vs. 1.99% explained variation) and was significant even when corrected for confounding technical effects. In addition, there was no evident clustering of samples by methodology in multidimensional scaling, indicating that the observed patterns are true biological signals (Table 1; Appendix 1: Fig. S2).

Fig. 2

Fractional abundance plot of symbiont ASVs within individual snails according to Alviniconcha and Ifremeria species.

Fig. 3

Mid‐point rooted iqtree consensus phylogeny of ASVs within symbiont genera. Node labels indicate ultra‐fast bootstrap support values.

Fig. 4

(A) Principal coordinate analysis plot based on weighted UniFrac distances. Data were normalized to proportions before analysis. Numbers in brackets indicate sample sizes for each host taxon. (B) Alpha diversity within host species based on Shannon's and Simpson's diversity index.

Table 1

Results for linear decomposition models (LDM) and PERMANOVAs based on weighted UniFrac distances.

	LDM				PERMANOVA
Source of variation	df	F	VE (%)	p	F	p
Model 1
Geographic region	8	2.2861	16.49	0.0001	324.710	0.0010
Host	3	1.3426	25.82	0.0001	704.879	0.0010
Model 2
Vent	2	3.5363	30.39	0.0001	4959.805	0.0010
Host	2	6.1004	52.42	0.0001	6424.714	0.0010
Model 3
Methodology	1	0.1179	1.99	0.0001	264.664	0.0010
Vent	2	3.5363	29.79	0.0001	9919.611	0.0010
Host	2	6.1004	51.38	0.0001	12849.428	0.0010

Three different models were run to assess the effects of DNA preservation, extraction and sequencing method on patterns of symbiont diversity: (1) Model including the complete dataset and controlling for effects of methodology, (2) Model restricted to A. boucheti, A. kojimai and A. strummeri from the ELSC and controlling for effects of methodology, (3) Model restricted to A. boucheti, A. kojimai and A. strummeri from the ELSC and including methodology as main explanatory factor. Sources of variation are shown in sequential order tested in the model. Significant sources of variation are indicated in bold. df = degrees of freedom, F = F statistic, VE = explained variation, p = p‐value.

Like A. kojimai and A. strummeri, most other host species were associated with particular lineages of thiotrophic Gammaproteobacteria. Alviniconcha adamantis was affiliated with symbionts of the genus Ca. Thiobios, whereas A. hessleri and I. nautilei hosted distinct Thiolapillus symbiont ASVs. Many I. nautilei individuals further harboured a minority methanotrophic symbiont from the genus Methylomonas, especially at vent sites within the Eastern Lau Spreading Center (ELSC). Only Alviniconcha boucheti and A. marisindica were dominated by different region‐specific campylobacterial ASVs of the genera Sulfurimonas or Sulfurovum.

Within geographic area, the gammaproteobacterial symbionts of A. kojimai and A. hessleri showed evidence for structuring by vent field (Appendix 1: Fig. S3), while no intra‐regional differentiation was observed or could be tested in symbionts of any other host species that we sampled from multiple localities (data not shown). However, this finding is likely an artefact of the limited resolution of the 16S rRNA marker gene. For example, recent metagenomic analyses indicate that symbiont populations of all host taxa from the Lau Basin are partitioned between vent sites (Breusing et al., 2021a; Breusing et al., 2021b). In contrast to the traditional view of microbial biogeography that poses that ‘everything is everywhere’ (Baas‐Becking, 1934), geographic subdivision of microbial symbionts appears to be common in a variety of marine symbioses, often exceeding that of the corresponding host populations (Ho et al., 2017; Gould and Dunlap, 2019; Davies et al., 2020; Ücker et al., 2021; Breusing et al., 2021a; Breusing et al., 2021b). Depending on the symbiotic system, these patterns might arise from local adaptation, contrasting dispersal limitations between hosts and symbionts, host ecological behaviour and/or differences in environmental transmission mode. Given the strong oceanographic barriers among back‐arc basins in the Western Pacific Ocean (Mitarai et al., 2016), the observed partitioning of host‐specific symbiont ASVs according to broader geographic area might be largely due to decreased symbiont dispersal opportunities (that are likely exacerbated by environmental differences between habitats). By contrast, symbiont structure within regions, where dispersal limitations appear to be weak (Mitarai et al., 2016), is probably predominantly driven by ecological factors, such as differences in depth or vent geochemistry (Breusing et al., 2021a; Breusing et al., 2021b). Indeed, in A. kojimai the observed partitioning of symbiont types by vent field was correlated with contrasting depth regimes (Appendix 2), which often aligns with gradients in fluid chemistry (Beinart et al., 2012). On the other hand, the strong latitudinal subdivision found for the Thiolapillus symbiont of A. hessleri might be explained by dispersal limitations as biophysical models indicate that the southern and northern parts of the Mariana Basin are largely isolated (Mitarai et al., 2016; Breusing et al., 2021a; Breusing et al., 2021b).

Our data suggest that other factors, such as host size, have a comparatively small influence on the diversity and composition of symbiont ASVs within host individuals. Despite significant associations of symbiont richness with host size, correlation coefficients were low, suggesting limited biological relevance of this factor on intra‐host symbiont diversity (Appendix 1: Fig. S4). These results were consistent independent of whether analyses were carried out across or within individual host species. For intra‐species analyses only correlations for A. kojimai and A. boucheti were significant, though weak (p ≦ 0.05; R 2 ≦ 0.09). In most cases individuals contained only one symbiont ASV in accordance with Sanger sequence analyses (Beinart et al., 2012, 2015), though in some individuals up to six ASVs were observed. Although our study lacks data from settling larvae and juveniles, these findings could indicate that symbiont acquisition in provannid snails follows a different process than in bathymodiolin mussels and is more similar to that in vestimentiferan tubeworms. Hydrothermal vent mussels remain competent for symbiont acquisition throughout their lifetime (Wentrup et al., 2014; Ansorge et al., 2019), which should favour increased symbiont diversity in older individuals as well as newly infected juveniles where symbiont sorting has not yet been completed. By contrast, vestimentiferan tubeworms obtain their symbionts exclusively in a narrow window after settlement during post‐larva metamorphosis (Nussbaumer et al., 2006). Symbiont diversity can thus be expected to be highest at that developmental stage, with little effect of host size on symbiont richness during later stages. Alternatively, our observations may indicate that 16S rRNA amplicon sequences do not provide enough strain‐level resolution to observe shifts in symbiont composition across development stages, and that metagenomic analyses of symbiont populations are necessary instead.

Symbiont richness differs between host species and individuals

Despite low impact of host size, Alviniconcha and Ifremeria exhibited notable variability in symbiont diversity, both among individuals and species (Fig. 4B). These patterns could result from differences in the availability and composition of free‐living symbiont lineages at the time of infection, subsequent mutations inside the host and/or host selection on particular strains. Among host taxa, A. adamantis and A. marisindica showed the lowest symbiont diversity, which is probably due to the fact that these species were each sampled from only a single vent site and were represented by relatively few individuals (Fig. 4B). Interestingly, A. hessleri displayed some of the highest alpha diversities, with up to six ASVs within single host individuals, despite its restricted geographic distribution and small sample size compared to some of the other Alviniconcha species included in our analyses. Maybe the wide variation of geochemical conditions in the Mariana Back‐Arc Basin (Trembath‐Reichert et al., 2019) allows for a greater range of micro‐niches, which could promote diversity in the free‐living symbiont pool. In this case, symbionts within this host species might have a higher functional diversity that could favour coexistence of multiple strains, as has recently been reported for bathymodiolin mussels, where hosts can carry up to 16 symbiont strains due to variation in metabolic gene content (Ansorge et al., 2019). Alternatively, some of the observed variations might reflect intra‐host mutations of a single or a few symbiont phylotypes post‐infection. In the absence of genomic data, this explanation seems likely as all A. hessleri symbiont ASVs were very similar to each other, with an average of 99.4% pairwise sequence identity.

Symbiont richness is not saturated

Although we analysed symbiont 16S rRNA composition in over 700 snail individuals, symbiont discovery did not reach saturation in our dataset (Fig. 5). The number of ASVs within A. hessleri and I. nautilei, which both host symbionts of the genus Thiolapillus, was closest to reaching a plateau, while ASV accumulation curves for all other species showed a steady increase (Fig. 5). This is an interesting finding given that A. hessleri and I. nautilei were sampled across a relatively restricted area compared to some of the other species (Appendix 1: Table S1). For other taxa that were represented by few individuals and geographic locations (e.g. A. adamantis, A. marisindica), but also those with widespread distributions (e.g. A. kojimai, A. boucheti), increased sampling efforts will probably reveal a currently hidden diversity of symbiont ASVs in the future. Consequently, while our dataset does not allow comparisons of diversification between symbiont genera or species at this time, more ASVs especially for some of the gammaproteobacterial taxa (e.g. unclassified Thiomicrospiraceae, Ca. Thiobios) will likely be recovered given the prevalence of gammaproteobacterial symbioses in provannid snails and other vent invertebrates (Dubilier et al., 2008).

Fig. 5

Symbiont ASV accumulation curves.

Conclusions

Here, we characterized the global diversity of chemosynthetic gill endosymbionts associated with species within the genera Alviniconcha and Ifremeria. As predicted by previous work, we found that each host species harboured 1–2 species‐ or genus‐level symbiont phylotypes. However, we were able to further assess strain‐level symbiont composition and diversity within and between individual snails by employing amplicon analysis of the 16S rRNA gene. In all host species, ASV accumulation curves indicated that the full diversity of symbionts associated with Alviniconcha and Ifremeria remains to be characterized. In most cases, symbiont ASV composition and richness were related to geographic range, with most ASVs detected in species where we sampled a large number of individuals across >10 geographically distant vent fields (e.g. A. kojimai and A. boucheti). An exception to this was A. hessleri, which had high symbiont richness and inter‐region symbiont structure despite a smaller sample size and much more modest geographic range, suggesting that these are not the only factors dictating symbiont composition and diversity. A more complete appraisal of the taxonomic and functional diversity of symbionts associated with Alviniconcha and Ifremeria will be critical to our understanding of the ecology and evolution of these genera, which have been assessed as ‘Endangered’ or ‘Vulnerable’ on the IUCN Red List (https://www.iucnredlist.org) due to imminent risks from deep‐seabed mining activities at hydrothermal vents in the Indian and Pacific oceans.

Experimental procedures

Sample collection and amplicon library preparation

Animal samples were obtained with remotely or human‐operated vehicles from 23 Indo‐Pacific vent localities that encompassed the global distributional range of species within the genera Alviniconcha and Ifremeria (Appendix 2; Fig. 1). Upon recovery of the samples, endosymbiont‐bearing gill tissue was dissected and frozen or stored in RNALater™ (Thermo Fisher Scientific, Waltham, MA, USA) at −80°C. DNA was purified with the Zymo Quick DNA 96 Plus and ZR‐96 Clean‐up kits (Zymo Research, Irvine, CA, USA) or the Qiagen DNeasy Blood & Tissue kit (Qiagen, Hilden, Germany). 2 × 250 bp paired‐end amplicon libraries for the 16S rRNA V4–V5 region were constructed with the 515F/926R primer pair (Walters et al., 2015) and sequenced to an average of 34 844 total reads on Illumina MiSeq and NovaSeq platforms at the Argonne National Laboratory (Lemont, IL, USA) and Novogene (Beijing, China) respectively (Appendix 2). Host species were identified through shell morphology (Laming et al., 2020) and subsequent sequencing of the mitochondrial COI gene with universal primers (Folmer et al., 1994; Geller et al., 2013).

Identification of amplicon sequence variants

We used the usearch v11 denoising pipeline (Edgar, 2010) to decompose merged, adapter‐clipped paired‐end reads into ASVs, imposing a merge length of 300–400 bp, a maximum error rate of 0.001 and a minimum base quality of 20. The taxonomic identity of each variant was determined in Qiime2 (https://qiime2.org) with a Naïve Bayes classifier trained against the V4–V5 region extracted from the silva 132 99% reference database as well as through blast+ searches against the NR database (Camacho et al., 2009). Only ASVs that had a match to a previously verified Alviniconcha or Ifremeria gill endosymbiont sequence were considered for further analysis. To assess potentially unrecovered variation in the symbiont dataset we applied the oligotyping v2.0 method (Eren et al., 2013). ASVs with less than 2.37% abundance in a sample were excluded to account for sample cross‐contamination (Minich et al., 2019). Phylogenetic relationships among ASVs were determined with the iqtree (Minh et al., 2020) plugin for Qiime2 based on 10 independent runs with each 5000 ultrafast bootstrap samples. Ultrafast bootstrap trees were optimized through the nearest neighbour interchange procedure with a perturbation strength of 0.2 and a stopping criterion of 200 trees.

16S rRNA diversity analyses

We used the phyloseq package in R v4.0.3 (McMurdie and Holmes, 2013; R Core Team, 2020) to assess symbiont 16S rRNA variation within and between hosts and geographic regions, excluding samples with less than 1000 reads to ensure statistical robustness. For alpha and beta diversity analyses symbiont abundances were normalized to proportions (McKnight et al., 2019). Metric and non‐metric multidimensional scaling plots were constructed based on weighted UniFrac distances. To verify that the distribution of ASV diversity is representative of real biological patterns and not technical artefacts from differences in methodology, we performed linear decomposition models (LDMs) and a modified version of PERMANOVA with the LDM package in R, as these methods have been shown to be relatively robust to variance in group dispersion (Hu and Satten, 2020). Analyses were run on both the full dataset and a data subset including only samples of Alviniconcha from the ELSC which were processed with a mixture of methods. PERMANOVAs and LDMs were conducted with 1000 and 10 000 maximum permutations respectively, with methodology included as either confounding variable or main explanatory factor. Relationships between number of ASVs and host size were determined based on Spearman rank correlations with the ggpubr package (Kassambara, 2020).

Data Availability

All bioinformatic scripts and final files for analysis are available on GitHub under https://github.com/cbreusing/Provannid_16S_SSU_meta-analysis. Raw 16S rRNA amplicon reads have been deposited in the Sequence Read Archive under BioProjects PRJNA473256, PRJNA473257, PRJNA610289, PRJNA610290, PRJNA763784 and PRJNA767887, while host COI sequences are available in GenBank under accession numbers listed in Appendix 2.

Appendix S1: Supporting information.

Click here for additional data file.

Appendix S2: Supporting information.

Click here for additional data file.