Allopatric and Sympatric Drivers of Speciation in Alviniconcha Hydrothermal Vent Snails

Abstract Despite significant advances in our understanding of speciation in the marine environment, the mechanisms underlying evolutionary diversification in deep-sea habitats remain poorly investigated. Here, we used multigene molecular clocks and population genetic inferences to examine processes that led to the emergence of the six extant lineages of Alviniconcha snails, a key taxon inhabiting deep-sea hydrothermal vents in the Indo-Pacific Ocean. We show that both allopatric divergence through historical vicariance and ecological isolation due to niche segregation contributed to speciation in this genus. The split between the two major Alviniconcha clades (separating A. boucheti and A. marisindica from A. kojimai, A. hessleri, and A. strummeri) probably resulted from tectonic processes leading to geographic separation, whereas the splits between co-occurring species might have been influenced by ecological factors, such as the availability of specific chemosynthetic symbionts. Phylogenetic origin of the sixth species, Alviniconcha adamantis, remains uncertain, although its sister position to other extant Alviniconcha lineages indicates a possible ancestral relationship. This study lays a foundation for future genomic studies aimed at deciphering the roles of local adaptation, reproductive biology, and host–symbiont compatibility in speciation of these vent-restricted snails.


Introduction
How novel species evolve in response to environmental pressures remains a fundamental question in evolutionary biology (Schluter 1998). New species can originate through a number of mechanisms, including geographic isolation (allopatric speciation) and formation of reproductive barriers due to localized natural selection (sympatric speciation) (Darwin 1859;Mayr 1942;Dobzhansky 1951;Coyne and Orr 2004;Schuler et al. 2016;Yang et al. 2017). The evolution of new species commonly occurs as a result of ecological or mutation-order speciation (Schluter 2009). Ecological speciation theory predicts that reproductive isolation should evolve between populations adapting to different environments (Mayr 1942) and is a consequence of divergent selection on functional traits enabling specialization for distinct ecological niches (Rundle and Nosil 2005;Schluter 2009). Habitat isolation can occur at different spatial scales that range from "macrospatial" (a form of allopatric speciation, where gene flow between populations is absent due geographic separation) to "microspatial" (where distinct populations co-occur geographically but rarely interbreed because of adaptations to different ecological niches) (Coyne and Orr 2004). By contrast, mutation-order speciation happens when reproductive isolation evolves due to random fixation of alleles in two populations that face the same selective pressures (Mani and Clarke 1990;Schluter 2009;Nosil and Flaxman 2011). Evolutionary divergence under this scenario is unlikely (though not impossible) if gene flow is present (Schluter 2009;Nosil and Flaxman 2011), whereas ecological speciation is common with or without gene flow (Schluter 2009;Nosil and Flaxman 2011;Feder et al. 2012;Campbell et al. 2018). Additional processes including inter-specific hybridization and host-microbe associations have received less attention but are now accepted as significant drivers of speciation (Schuler et al. 2016).
Our current understanding of speciation mechanisms in animals stems mostly from evolutionary research on terrestrial and freshwater organisms, whereas less is known about speciation in oceanic environments, in particular deep ocean basins (Miglietta et al. 2011). Processes that lead to evolutionary divergence in marine systems may differ from other environments given that physical barriers are less pronounced and deep-sea habitats might have been less affected by past climatic fluctuations (Miglietta et al. 2011;but see Vrijenhoek 2013), at least over much of the Cenozoic. Nonetheless, deep-sea hydrothermal vents provide an exception to other deep-sea environments due to their global but island-like distribution and their relative instability on ecological and geological time scales (Vrijenhoek 2013). The dense animal communities living in these habitats are exposed to highly variable and extreme environmental conditions, and many of the invertebrates are sustained through obligate symbiosis with chemosynthetic bacteria (Dubilier et al. 2008). The phylogenetically and ecologically closely related snail genera Alviniconcha and Ifremeria (Gastropoda: Abyssochrysoidea) are among the dominant inhabitants of Indo-Pacific hydrothermal vents. The genus Ifremeria is monotypic and apparently restricted to Western Pacific vents in the Manus, Lau, and North Fiji basins. In contrast, the broadly distributed genus Alviniconcha comprises six species : Alviniconcha kojimai and Alviniconcha boucheti co-occur in the Manus, Lau, and North Fiji basins; Alviniconcha strummeri is known from the southern Lau Basin; Alviniconcha hessleri and Alviniconcha adamantis appear to be endemic to the Mariana Back-Arc and Volcanic Arc, respectively; and Alviniconcha marisindica occurs along the Central Indian Ridge ( fig. 1 and table 1).
Alviniconcha and Ifremeria are nutritionally dependent on horizontally transmitted, chemosynthetic sulfide-oxidizing endosymbionts that are acquired from a diverse community of environmental bacteria and are harbored within vacuoles in the gill filaments (Windoffer and Giere 1997;Suzuki, Sasaki, Suzuki, Nogi, et al. 2006;Trembath-Reichert et al. 2019). Except for A. boucheti and A. marisindica, which are constrained to symbionts of the class Campylobacteria (formerly Epsilonproteobacteria; Waite et al. 2017), all other Alviniconcha species and Ifremeria nautilei typically host Gammaproteobacteria symbionts (Urakawa et al. 2005;Suzuki, Sasaki, Suzuki, Nogi, et al. 2005;Suzuki, Sasaki, Suzuki, Tsuchida, et al. 2005;Suzuki, Kojima, Sasaki, et al. 2006;Beinart et al. 2012). In zones of sympatry, the two snail genera occupy relatively narrow ranges of chemical and thermal conditions that are well within their physiological tolerance limits and that do not support maximal rates for chemoautotrophic growth (Sen et al. 2013). This observation suggests that interspecies competition is an important factor determining the realized niche of these snails (Sen et al. 2013). Beinart et al. (2012) proposed that host niche utilization is strongly influenced by the availability and metabolic capacity of locally dominant symbiont phylotypes. In their study of co-occurring Alviniconcha species from the Lau Basin, they showed that the distribution and abundance of host-symbiont combinations (holobionts) followed a latitudinal gradient in vent geochemistry, wherein A. boucheti holobionts dominated at northern localities characterized by high end-member concentrations of sulfide and hydrogen ($3.5-7 mM and $100-500 mM, respectively), A. kojimai holobionts dominated at midlatitude localities characterized by medium endmember concentrations ($2.5-4 mM and $50-100 mM, respectively), and A. strummeri holobionts dominated at southern localities characterized by low end-member concentrations of these reductants ($1-3 mM and 35-135 mM, respectively). These broad-scale patterns of FIG. 1. Map of sampling localities for the six Alviniconcha species analyzed in this study. Breusing et al. . doi:10.1093/molbev/msaa177 MBE microbe-associated habitat segregation suggest that the endosymbionts might have played a significant role in the diversification of these gastropods.
Despite their phylogenetic and ecological affinities, Ifremeria and Alviniconcha have undergone remarkably different degrees of speciation. In the present study, we used fossil-calibrated molecular phylogenies and population genetic analyses based on mitochondrial and nuclear markers to obtain insights into allopatric and sympatric mechanisms that led to the high degree of diversification in the genus Alviniconcha. To assess the role of allopatric isolation, we compared the phylogenetic divergence events among the Abyssochrysoidea with Mesozoic and Cenozoic geological events of the Indo-Pacific Ocean basins and examined genetic subdivision within several broadly distributed species.
To determine evidence for sympatric isolation, we compared the time-calibrated phylogenies with the diversity and phylogenetic relationships of the snails' chemosynthetic endosymbionts and assessed the potential for hybridization among co-occurring Alviniconcha species.

Gene Networks
Haplotype networks for the two mitochondrial genes, COI and 12S, indicated complete sorting among lineages for the six Alviniconcha species and their sister taxon I. nautilei ( fig. 2), as previously reported . The COI sequences from Alviniconcha and Ifremeria differed by at least 75 mutations, whereas pairs of Alviniconcha species differed by   (table 2). Note, however, that D values for the first ABBA-BABA test were close to -1, and it is possible that gene flow between A. kojimai and A. boucheti might have been detected with a larger genomic marker set.

Fossil-Calibrated Phylogeny and Species Tree
Fossil-calibrated Beast analyses based on ribosomal, mitochondrial, and histone marker genes revealed that the Abyssochrysidae is a relatively old family ( fig. 4) In general, the topology of the species tree for Alviniconcha agreed with that of the concatenated gene tree, although the position of A. adamantis indicated a potential sister relationship to A. boucheti and A. marisindica ( fig. 4). The second species tree including the nuclear markers EF1a, ATPSb, and ATPSa resulted in an unresolved polytomy and we were unable to obtain robust estimates for phylogenetic placements (supplementary fig. S7, Supplementary Material online).

Phylogenetics of Endosymbionts
Molecular Bayesian analyses of the full-length 16S rRNA revealed geographical partitioning for the endosymbionts of Alviniconcha and Ifremeria, suggesting that symbiont phylotypes associated with I. nautilei, A. boucheti, and A. kojimai were genetically distinct between the Lau, Manus, and North Fiji basins ( fig. 5). The endosymbionts of I. nautilei were most closely related to free-living Thiolapillus brandeum and the Gammaproteobacteria symbiont of A. hessleri. Together, these bacterial groups formed a neighboring clade to the "c-1" endosymbionts of A. kojimai (Beinart et al. 2012) and free-living Thiomicrospira and Hydrogenovibrio ( fig. 5A). In both the Lau and Manus basins, A. kojimai additionally hosted Campylobacteria lineages that were also found in A. boucheti. The endosymbionts of A. adamantis formed a distinct clade that was most closely related to the "c-Lau" group present in A. strummeri (Beinart et al. 2012).
Alviniconcha marisindica and A. boucheti were the only species that harbored predominantly Campylobacteria symbionts ( fig. 5B). Although A. marisindica contained a single 16S rRNA phylotype related to free-living Sulfurovum, A. boucheti associated with distinct Campylobacteria lineages that were related to free-living Sulfurimonas, with one exception: A. boucheti from American Samoa harbored a phylotype that was most closely related to the A. marisindica symbiont.
Character traces of symbiont type identified the Gammaproteobacteria endosymbiosis as the ancestral condition given that the evolutionary older host species harbored only Gammaproteobacteria ( fig. 6). By contrast, the Campylobacteria endosymbiosis appears to be a derived trait that evolved multiple times independently.

Symbiont Phylotypic Diversity and Partitioning
Amplicon sequencing of the 16S V4-V5 hypervariable regions generated an average of 21,302 paired-end reads/sample. These were reduced to $9,613 reads/sample after read merging. Sequence denoising resulted in 85 zero-radius operational taxonomic units (ZOTUs), 11 of which were verified as Alviniconcha or Ifremeria endosymbionts. Oligotyping of amplicon reads comprising these symbiont ZOTUs identified a total of 23 symbiont 16S rRNA oligotypes in the data set. For 337 snail samples, we were able to obtain both host species identities and 16S oligotypes. Only these individuals were, Each host individual contained one majority symbiont phylotype together with up to 12 minority variants, with A. hessleri exhibiting the highest and A. boucheti, A. kojimai, and A. strummeri exhibiting the lowest symbiont richness ( fig. 7). In general, symbiont phylotypes corresponded to host species, except in the case of A. kojimai and A. strummeri, where symbiont compositions largely overlapped (figs. 7 and 8). All host species, with the exception of A. boucheti, comprised mostly a combination of different Gammaproteobacteria lineages. Alviniconcha boucheti was dominated by Campylobacteria symbiont phylotypes, although in the Lau Basin this species also harbored Gammaproteobacteria symbionts that were typically associated with A. kojimai and A. strummeri. However, due to presence of the A. kojimai symbiont (Oligo1) in negative controls on the same well plate where these A. boucheti samples were processed, we cannot exclude the possibility that some of these observations resulted from well-to-well contamination. In accordance with our phylogenetic analyses, the symbionts associated with A. boucheti were partitioned across basins.
In most cases, we did not have sufficient sample sizes to reveal potential within-basin structuring of the Gammaproteobacteria or Campylobacteria symbiont phylotypes, with the exception of the A. hessleri symbionts, which seemed to form two latitudinal clusters ( fig. 8).
Partial Mantel tests supported the described patterns, indicating significant correlations between symbiont and host genetic distances (r ¼ 0.76; P ¼ 1e-04) as well as between symbiont genetic distance and geographic distance (r ¼ 0.55; P ¼ 1e-04). By contrast, no association between host genetic distance and geographic distance was found.

Discussion
In the present study, we performed fossil-calibrated phylogenetic reconstructions and population genetic analyses to illuminate mechanisms of speciation in provannid snails of the genus Alviniconcha, which are dominant members of deepsea hydrothermal vent ecosystems in the Western Pacific and Indian Ocean.
Major tectonic and climatic changes in the Indo-Pacific during the Mesozoic and Cenozoic probably played a central  Speciation in Hydrothermal Vent Snails . doi:10.1093/molbev/msaa177 MBE role for geographic isolation during the early evolution of Alviniconcha snails. Our fossil-calibrated phylogeny indicates that Ifremeria and Alviniconcha split 96-131 Ma close to the peak of the Aptian oceanic anoxic event (Li et al. 2008). Although anoxic events led to widespread declines of marine biodiversity, they also created novel opportunities for habitat MBE specialization and allopatric speciation within oxygendepleted deep-sea environments, such as hydrothermal vents (Rogers 2000). Similarly, expansion of oxygen minimum zones might have contributed to the evolution of the most recent common ancestor of extant Alviniconcha, which originated during the upper Eocene ($50 Ma) shortly after the Paleocene-Eocene Thermal Maximum, another period that was characterized by increasing anoxia in the deep sea. Around the same time, the Philippine Sea Plate developed above the Manus hotspot, and a wide deep-water passage connected the Indian and Pacific Oceans (Hall et al. 1995;Hall 2002; see figs. 14-25 in Hall 2002). This passage started to narrow around 30 Ma and eventually closed 25 Ma, when the Australian Plate subducted under the Philippine Sea Plate and created physical barriers to dispersal (Gaina and Müller 2007).
Gradual vanishing of the Indo-Pacific deep-water connection is consistent with the timing of the most recent common ancestors of the two major Alviniconcha clades and likely promoted their diversification: 1) A. marisindica and A. boucheti ($30-46 Ma) and 2) A. hessleri, A. kojimai, and A. strummeri ($18-31 Ma). A variety of along-arc splitting and rifting events in the Miocene to Pliocene subsequently led to the formation of the North Fiji Basin ($12 Ma), Mariana Back-Arc Basin ($7 Ma), and Lau Basin ($5 Ma) (Hall 2002;Crawford et al. 2003;Straub et al. 2015). Evolution of the North Fiji Basin and Mariana Back-Arc Basin corresponds well with the divergence of A. hessleri and A. kojimai ($10 Ma), suggesting a potentially allopatric origin of these lineages that might have allowed A. hessleri to invade deep-water vent sites of the Mariana Back-Arc, where it is nowadays an endemic species. Alviniconcha adamantis is the only other Alviniconcha species that occurs in the Mariana region, but it is restricted to shallow vents on volcanoes of the Izu-Bonin-Mariana Arc, which have been active in some form for 40 My (Stern et al. 2003). Interestingly, the phylogenetic position of A. adamantis was not well supported and formed a basal branch to the two other clades of Alviniconcha. Although we do not have enough data to determine its placement, the A. adamantis lineage might be ancestral to other Alviniconcha lineages. If true, this genus might have followed a gradual transition from shallow to deep-water vent sites, as has been proposed for bathymodiolin mussels (e.g., Jones et al. 2006;Miyazaki et al. 2010; but see Thubaut et al. 2013).
Although our results suggest that splits between some Alviniconcha species involved geographic subdivision, allopatric divergence alone might be insufficient to explain all instances of speciation. For example, A. boucheti, A. kojimai, and A. strummeri live sympatrically in the Lau Basin, and A. boucheti and A. kojimai also coinhabit vents in the North Fiji and Manus basins. For both of these species, genetic structure analyses indicated varying degrees of genetic

Alviniconcha adamantis
Alviniconcha hessleri Speciation in Hydrothermal Vent Snails . doi:10.1093/molbev/msaa177 MBE subdivision not only among but also within disparate backarc basins. However, based on recent biophysical models for the Western Pacific region (Mitarai et al. 2016), ocean circulation should support within-basin connectivity, especially for species that are assumed to have planktotrophic larvae with high dispersal potentials, such as Alviniconcha (War en and Bouchet 1993). These contrasting findings suggest that antagonistic ecological factors might exert strong selection against immigrants with maladapted phenotypes, thereby counterbalancing the homogenizing effect of dispersal and allowing speciation in the presence of gene flow. The present genetic data do not allow us to confidently distinguish whether the observed speciation patterns occurred during primary sympatry or resulted from secondary colonization following allopatry. At least the ancestors of A. kojimai and A. strummeri likely co-occurred in the predecessor of the modern Manus Basin and, thus, may have evolved sympatrically. Nevertheless, opening of the younger Lau and North Fiji basins could have allowed for brief geographic separation (<5-10 Ma), if species invaded these regions asynchronously. Hydrothermal vents are environmentally dynamic habitats that offer ample opportunities for niche segregation along steep thermal and geochemical gradients. Geochemical patchiness and the availability and/or fitness of symbiont phylotypes (rather than functional host anatomy) are believed to affect spatial partitioning among Alviniconcha species in the Eastern Lau Spreading Center (Beinart et al. 2012;Laming et al. 2020). Alviniconcha boucheti, which typically hosts Campylobacteria symbionts, was most abundant at northern vent localities characterized by high sulfide and hydrogen concentrations. On the other hand, A. kojimai and A. strummeri, which usually host Gammaproteobacteria symbionts, were more frequent at southern vent localities with lower sulfide and hydrogen concentrations. Functional differences among the symbiont types combined with genetically based host-symbiont specificities could explain the broadscale segregated distribution of these species and might have provided a mechanism for their diversification. Free-living chemosynthetic Campylobacteria are known to encode sulfide-oxidizing enzymes that are most effective at high H 2 S concentrations (Han and Perner 2015), providing a selective advantage over Gammaproteobacteria at sulfide-rich vents (Nakagawa and Takai 2008;Yamamoto and Takai 2011;Patwardhan et al. 2018). Recent physiological experiments and transcriptomic analyses imply that this might also be the case for the Lau Basin Alviniconcha symbionts (Breusing et al. 2020). Possible relationships between vent geochemistry, symbiont phylotypes, and host species have not been similarly explored in other back-arc basins. Nevertheless, the PACMANUS and Vienna Woods vent fields in the Manus Basin are known to exhibit physicochemical gradients similar to those found in the Lau Basin (Reeves et al. 2011), providing comparable opportunities for symbiontrelated niche differentiation. Our partial Mantel tests and NMDS analyses involving Alviniconcha species from the Lau and Manus basins, and the Mariana Back-Arc revealed a significant association of symbiont genotypes with host species and vent localities, suggesting that symbiont-mediated habitat segregation might be common in chemosynthetic vent snail symbioses.

Alviniconcha boucheti Vanuatu
Together, our host phylogenetic and bacterial phylotype analyses imply a possible role for endosymbionts in

MBE
Alviniconcha speciation. For instance, A. boucheti and A. marisindica comprise a monophyletic clade associated predominantly with Campylobacteria symbionts. A trace of character trait evolution indicated that Campylobacteria endosymbionts were acquired in the most recent common ancestor of these two host species, possibly prior to geological events that contributed to splitting of the A. boucheti-A. marisindica clade. It is possible that a switch to a new symbiont type offered novel ecological opportunities, thereby promoting adaptive divergence from the Gammaproteobacteria-dominated Alviniconcha lineages. Both A. boucheti and A. marisindica associate with symbionts related to free-living Sulfurovum, whereas A. boucheti also associates with symbionts of the genus Sulfurimonas, indicating that the Sulfurovum endosymbiosis might be the evolutionarily older condition in these two host species. Speciation due to symbiont-driven niche expansion has been proposed for deep-sea mytilid mussels, with similar time intervals of diversification as observed for provannid snails (Lorion et al. 2013). Interestingly, geographic partitioning of the A. boucheti symbiont phylotypes corresponded broadly with the genetic structure of the host populations. Although we cannot exclude the hypothesis that these covarying patterns were mediated by geographic separation, gene flow between basins is expected on evolutionary time scales (Mitarai et al. 2016); therefore, genetically based host-symbiont specificities might also shape the diversification of both partners. Presently, no substantive evidence exists for past hybridization between broadly co-occurring A. boucheti and A. kojimai (although future genomic studies may provide other insights into this question). Host-symbiont incompatibilities might contribute to hybrid dysgenesis, a scenario that has been described in insect-microbe symbioses (e.g., Feder et al. 1994;Bordenstein 2012, 2013). Although A. kojimai occurs in association with Campylobacteria symbionts, these observations were relatively rare, indicating a bias toward a Gammaproteobacteria endosymbiosis (Beinart et al. 2012). Microbial symbionts are currently considered the most disregarded element in animal and plant evolution (Shropshire and Bordenstein 2016). A recent study has shown that both transmission mode and symbiont function influence the degree of host dependence that can ultimately lead to major evolutionary transitions and innovations (Fisher et al. 2017). Although strong host dependence is usually correlated with vertical transmission mode, the same effect can be achieved in horizontally transmitted symbioses if symbionts provide an essential nutritional benefit (Fisher et al. 2017). Consequently, although further studies are needed to characterize the role of symbionts in host speciation, their impact is likely of particular significance in chemoand photosynthetic ecosystems where host fitness is tightly linked to symbiont metabolism, such as hydrothermal vents, cold seeps, and coral reefs.
An alternative, though not mutually exclusive, scenario for sympatric speciation in Alviniconcha could be rapid reproductive isolation and sexual selection. In many organisms, including marine snails of the genera Haliotis and Tegula, egg and sperm proteins are under strong positive selection and evolve at an accelerated rate that could ultimately split populations into reproductively isolated species (Galindo et al. 2003;Panhuis et al. 2006;Hellberg et al. 2012). Subsequent inter-specific gene flow could be reduced further due to assortative mating, which is often observed in gastropods with internal fertilization (Johannesson et al. 2008;Zahradnik et al. 2008). Although little is known about the reproductive biology of Provannidae, internal fertilization is the implied reproductive mode in this family (Gustafson and Lutz 1994).
Our study suggests that both geographical and ecological divergence (followed by putative reproductive isolation) played significant roles in the evolutionary diversification of the six extant Alviniconcha lineages at Indo-Pacific hydrothermal vents. It remains to be determined why the same phenomena did not similarly affect speciation within the genus Ifremeria, which is represented by only a single species. The bounded hypothesis on diversification proposes that the evolution of new species is constrained by competition for available niche space (e.g., Rabosky 2009;Price et al. 2014;Larcombe et al. 2018) and that an increase in the number of species within one clade may be compensated by a decline in another (Ricklefs 2010). Given that Ifremeria has similar ecological requirements as Alviniconcha (Podowski et al. 2010), it is possible that competitive exclusion and niche filling might have limited diversification within this genus. Future genomic and ecological studies may be helpful to illuminate some of these questions and to assess the relative importance of host-symbiont incompatibilities, local adaptation, and reproductive barriers in the evolutionary history of these gastropods.

Sample Collection and DNA Methods
Snail samples were collected from 32 localities in the Western Pacific and Indian Oceans (table 1). DNA extraction, general polymerase chain reaction (PCR) conditions, amplicon purification, and DNA sequencing used methods that were reported for provannid snails (Johnson et al. 2010. Twelve primer sets were used to amplify three mitochondrial (COI, mt-12S, and mt-16S) and seven nuclear (H3, 28S-D1, 28S-D6, 18S, ATPSa, ATPSb, and EF1a) DNA targets from vent snail hosts as well as the 16S rRNA locus from the snails' chemosynthetic symbionts (table 3). The 16S rRNA gene of gill endosymbionts of two A. adamantis and five I. nautilei individuals was cloned prior to Sanger sequencing (Goffredi et al. 2007) with 24 colonies per individual. All PCR products were diluted in 40-50 ll of sterile water and purified with a Multiscreen HTS PCR 96 vacuum manifold system (Millipore Corp., Billerica, MA). PCR products were sequenced bidirectionally on an ABI 3130XL sequencer with BigDye Terminator v3.1 chemistry (Life Technologies Corp., Carlsbad, CA) and primers used in PCR. Samples from American Samoa were added late to the study and were only included in phylogenetic analyses of the symbionts. Forward and reverse sequences were trimmed, merged, and manually curated in Geneious Prime (http://www.geneious.com; last accessed July 24, 2020). Vent Snails . doi:10.1093/molbev/msaa177 MBE Heterozygous sites at nuclear loci were phased with Phase v2.1.1 (Stephens and Donnelly 2003) as in Breusing et al. (2015).
Structure v.3.5 (Pritchard et al. 2000) was applied to investigate the number of genetic clusters K in the total Alviniconcha sample set and to test for contemporary admixture in the co-occurring species A. boucheti, A. kojimai, and A. strummeri. Analyses were run for all six species together as well as for each species separately using 1) mitochondrial and nuclear data and 2) nuclear data only. Allelic information was provided as a number code. Inter-specific Structure analyses were based on five (COI, H3, EF1a, ATPSa, and ATPSb) and four genetic loci (H3, EF1a, ATPSa, and ATPSb), respectively. The same genetic markers were used in intra-specific analyses, except that we excluded the H3 locus given that it was essentially monomorphic within Alviniconcha populations. Model runs and settings for Structure generally followed Breusing et al. (2016) with ten repetitions for each K between 1 and 8 for inter-specific analyses and between 1 and 6 for intra-specific analyses. We subsequently used Clumpak v.1.1 (Kopelman et al. 2015) to integrate results across independent runs for each K value and determine the optimal clustering solution. We initially used the DK method (Evanno et al. 2005) to objectively infer the number of genetic clusters in the data set. However, as this method is biased toward identifying the uppermost hierarchical level of population structure (i.e., lowest parsimonious K value) (Cullingham et al. 2020), we visually inspected all K plots for biologically meaningful structuring and chose the highest K value where genetic subdivision could be observed (Meirmans 2015).
To distinguish hybridization from incomplete lineage sorting in the sympatric species A. boucheti, A. kojimai, and A. strummeri, we performed two ABBA-BABA tests on the concatenated nuclear gene alignments of 90 snail individuals using HybridCheck (Ward and van Oosterhout 2016). Gene flow was assessed between 1) A. boucheti (P3) and A. kojimai/ A. strummeri (P1 þ P2) and 2) A. strummeri (P3) and A. kojimai/A. hessleri (P1 þ P2), with A. adamantis and A. boucheti as outgroups, respectively. Sequences from I. nautilei were not used because of missing data. Jackknife sampling was performed on four blocks of 305 bp (i.e., the approximate size of each nonlinked locus).

Fossil-Calibrated Phylogeny of the Abyssochrysidae
All Alviniconcha species were included in a fossil-calibrated phylogenetic reconstruction of the Abyssochrysidae with the program Beast v.2.5.7 (Bouckaert et al. 2014). Fossil calibrations included the divergence of Neptunea amianta and N. antiqua, which have probably been separated since the earliest appearance of the genus in the Late Eocene, 33-37 Ma (Amano 1997), as well as Provanna and Desbruyeresia, which are known from Cenomanian (Cretaceous) seep deposits in Japan, 93-100 Ma (Kaim et al. 2008). These estimates can be considered a minimum age of the split between the genera. As sequences for the nuclear DNA markers EF1a, ATPSa, and ATPSb could not be obtained for fossil samples, only ribosomal and mitochondrial markers plus the nuclear H3 locus were used for analysis (table 3). New sequences for all species of Alviniconcha from several localities were added to existing alignments from Johnson et al. (2010Johnson et al. ( , 2015 with default parameters using Muscle executed in Geneious Prime (http://www.geneious.com; last accessed July 24, 2020). Fragments of the 16S mtRNA and 12S mtRNA loci that had been previously excluded due to alignment ambiguities with more distantly related taxa were informative and therefore included within the genera Alviniconcha and Ifremeria. For the rest of the Abyssochrysidae, these regions are represented as missing data. For each locus and species, the most common allele was chosen as representative sequence for phylogenetic reconstructions. JModelTest v2 (Darriba et al. 2012) was used to identify the best fitting evolutionary model for each gene and to partition the data for phylogenetic analyses. Concatenated gene trees were estimated based on a relaxed, lognormal clock with a Yule tree prior, a GTR þ I þ C substitution model as well as unlinked base frequencies and rates. The Markov chain Monte Carlo chains were run for 200 million generations, with parameters sampled every 5,000 generations. We used Tracer v1.7.1 (Rambaut et al. 2018) to determine if effective sample sizes were adequate for all parameters and if analyses had reached convergence. Node dates were calibrated as a normal distribution around the timing of genera splits. A StarBeast v2.5.7 (Bouckaert et al. 2014) species tree was estimated including only the genera Alviniconcha and Ifremeria as replicate individuals for more distantly related taxa were missing. These estimated trees included at least one individual from each population per species or replicates from a single population. Species trees were calculated based on 1) the same loci used for the concatenated gene tree as well as 2) all available genetic markers (i.e., including the nuclear genes EF1a, ATPSa, and ATPSb). Phylogenetic reconstructions were performed with loci unlinked by partition, a GTR þ I þ C site model with empirical frequencies, an uncorrelated lognormal clock, an analytical population size integration model, and a Yule tree prior. Markov chain Monte Carlo analyses were run for 200 million generations, with a tree sampled every 10,000 generations, a preburnin of 1,000 and 10 initialization steps. All analyses were run multiple times and results assessed in Tracer, FigTree v1.4.4 (http://tree.bio. ed.ac.uk/software/figtree/; last accessed July 24, 2020), and Densitree v2.2.1 (Bouckaert 2010).

Phylogenetic Analyses of Snail Endosymbionts
New 16S rRNA endosymbiont sequences from the gills of A. adamantis and I. nautilei were compared against the NCBI nonredundant database. The new sequences were then aligned with Sina v1.2.11 (Pruesse et al. 2012) against the global Silva 16S rRNA alignment with available 16S rRNA endosymbiont sequences from GenBank for Alviniconcha species, I. nautilei, and other representative chemosynthetic taxa (Beinart et al. 2012(Beinart et al. , 2015. Phylogenetic trees for the Gammaproteobacteria and Campylobacteria symbionts were constructed with MrBayes v.3.2.7a (Altekar et al. 2004) under a GTR þ I þ C substitution model in the Cipres Science Gateway v3.3 (Miller et al. 2010). Two sets of three heated chains and one cold chain were run for 1,100,000 generations and sampled every 100 generations. The first 100,000 generations were discarded as burnin. Reference sequences of Nautilia lithotrophica (NR_025497) and Beggiatoa alba (NR_041726) were used as outgroups for the Campylobacteria and Gammaproteobacteria trees, respectively. Ancestral endosymbiosis traits were identified with Mesquite v.2.75 (Maddison WP and Maddison DR 2011). High-Throughput 16S V4-V5 Amplicon Sequencing and Analysis The 16S V4-V5 hypervariable region was investigated to assess symbiont phylotypic variation in relation to host species and geographic location. After overnight lysis in proteinase K, DNA from RNALater (Thermo Fisher Scientific, Inc.) preserved gill samples was extracted with the Zymo Quick DNA 96 Plus kit (Zymo Research, Inc.) and further cleaned with the ZR-96 Clean-up kit (Zymo Research, Inc.) following the manufacturer's protocols. Barcoded amplicon libraries (515 samples and 62 negative controls) targeting a $375bp fragment were prepared using the primer pairs 515F/ 926R (Walters et al. 2015) and sequenced on an Illumina MiSeq platform with a 2Â 250-bp protocol at the Argonne National Laboratory (Lemont, IL). Host species identity was determined based on morphology (Ifremeria) or molecular barcoding of the mitochondrial COI gene (Alviniconcha).
After quality-checking the demultiplexed reads with FastQC v0.11.5 (Andrews 2010; http://www.bioinformatics. babraham.ac.uk/projects/fastqc/; last accessed July 24, 2020), we clipped adapters with Trimmomatic v0.36 (Bolger et al. 2014) and followed the Usearch v11 (Edgar 2010) pipeline for analysis. Reads were merged allowing lengths >300 bp, and then filtered applying a maximum error rate of 0.001, a quality threshold of 20, and a minimum length of 300 bp. The fast-x_uniques and unoise3 commands were used to dereplicate, denoise, and cluster the reads into ZOTUs, whereas the otutab command was used to generate the OTU table. Taxonomy was assigned to each ZOTU in Qiime2 (https:// qiime2.org; last accessed July 24, 2020) with a Naïve Bayes classifier that we trained against the target 16S V4-V5 region using the Silva 132 99% reference database. Ambiguous annotations were resolved using BlastN searches (Altschul et al. 1990). Recent studies have shown that well-to-well contamination accounts on average for 2.37% of reads in high biomass samples that were processed at the Argonne National Laboratory (Minich et al. 2019). Consequently, for further analysis, we removed ZOTUs that had <2.37% abundance in a sample, ignored all ZOTUs that were not annotated as vent snail endosymbionts and we excluded all samples with <1,000 symbiont reads. To identify any unrecovered phylotype variation, we analyzed the filtered symbiont read set with the Oligotyping v2.0 method (Eren et al. 2013) based on 24 informative nucleotide positions. Only oligotypes with at least 2.37% abundance in a sample (-a), 100 supporting reads (-M), and occurrence in more than one individual (-s) were considered. We further excluded reads that had base quality scores <20 in the nucleotide positions of interest (-q).
Nonmetric multidimensional scaling and diversity analyses for the final symbiont oligotype data set were performed with the PhyloSeq package in R v3.5.2 (McMurdie and Holmes 2013; R Core Team 2018) based on weighted Unifrac distances on log-transformed read counts. We also performed partial Mantel tests using the Vegan package (Oksanen et al. 2019) to identify relationships between Alviniconcha host and symbiont genetic distances and geography based on Pearson's product-moment correlations. Pairwise host Speciation in Hydrothermal Vent Snails . doi:10.1093/molbev/msaa177 MBE genetic distances were calculated for mitochondrial COI sequences based on Kimura's two-parameter model (Kimura 1980). For geographic distances, we determined the geodesics between vent sites using the Geosphere package (Hijmans 2019).

Supplementary Material
Supplementary data are available at Molecular Biology and Evolution online.