Host hybridization as a potential mechanism of lateral symbiont transfer in deep‐sea vesicomyid clams

Abstract Deep‐sea vesicomyid clams live in mutualistic symbiosis with chemosynthetic bacteria that are inherited through the maternal germ line. On evolutionary timescales, strictly vertical transmission should lead to cospeciation of host mitochondrial and symbiont lineages; nonetheless, examples of incongruent phylogenies have been reported, suggesting that symbionts are occasionally horizontally transmitted between host species. The current paradigm for vesicomyid clams holds that direct transfers cause host shifts or mixtures of symbionts. An alternative hypothesis suggests that hybridization between host species might explain symbiont transfers. Two clam species, Archivesica gigas and Phreagena soyoae, frequently co‐occur at deep‐sea hydrocarbon seeps in the eastern Pacific Ocean. Although the two species typically host gammaproteobacterial symbiont lineages marked by divergent 16S rRNA phylotypes, we identified a number of clams with the A. gigas mitotype that hosted symbionts with the P. soyoae phylotype. Demographic inference models based on genome‐wide SNP data and three Sanger sequenced gene markers provided evidence that A. gigas and P. soyoae hybridized in the past, supporting the hypothesis that hybridization might be a viable mechanism of interspecific symbiont transfer. These findings provide new perspectives on the evolution of vertically transmitted symbionts and their hosts in deep‐sea chemosynthetic environments.


| INTRODUC TI ON
Host-microbe symbioses are universal phenomena that are now considered key drivers of evolutionary innovation (Archibald, 2015;Brucker & Bordenstein, 2012;McFall-Ngai, 2008;McFall-Ngai et al., 2013). Over the past several decades, it has been established that symbiotic associations led to the evolution of cellular organelles and eukaryotic cell life (Archibald, 2015), while recent studies have emphasized their role in the formation of species (Brucker & Bordenstein, 2012). Transmission of microbes contributes to maintenance of symbiotic relationships across host generations and differences in transmission modes have important implications for the evolution of both partners (Bright & Bulgheresi, 2013;Vrijenhoek, 2010). Under vertical transmission, symbionts are inherited through the maternal and/or-less frequently-the paternal germ line (e.g., Ebert, 2013;Moran & Dunbar, 2006;Watanabe, Yukuhiro, Matsuura, Fukatsu, & Noda, 2014). In the predominant case of uniparental maternal inheritance, symbiont and mitochondrial genomes are cotransmitted and are thus genetically and evolutionarily linked. Bottleneck effects during transovarial transmission strongly reduce the effective size and genetic diversity of the symbiont population within individual hosts, thereby increasing the fixation of slightly deleterious mutations (Rispe & Moran, 2000;Vrijenhoek, 2010). Since recombination with environmental bacteria is limited and certain symbiont gene functions become obsolete or are complemented by the host or secondary symbiotic microbes, vertically transmitted symbionts typically lose genes through drift and selection, resulting in significant reductions in genome size (Bennett & Moran, 2015;Moran, McCutcheon, & Nakabachi, 2008;Sloan & Moran, 2012). Apart from vertical transmission, symbionts can be transmitted horizontally, either through uptake of free-living strains in the environment or through direct transfer between hosts (Bright & Bulgheresi, 2013;Ebert, 2013;Vrijenhoek, 2010). Because symbionts are acquired from a potentially diverse mixture of bacterial strains each generation, horizontal transmission often results in genetic heterogeneity in the symbiont population and the absence of co-evolution between host and symbiont. In contrast to vertically transmitted symbionts, horizontally transmitted symbionts switch between intra-and extrahost life phases, which increases rates of recombination and selective pressures for retaining genes necessary for surviving outside the host environment (Vrijenhoek, 2010). In various cases, it has been shown that horizontal transmission can supplement the vertical transmission mode (Ebert, 2013), thereby providing opportunities for recombination that can counteract the ongoing genome degradation in vertically transmitted symbionts (Vrijenhoek, 2010). Despite growing research on diverse host-microbe relationships, the mechanisms of symbiont transmission and their evolutionary consequences are still poorly understood (Bright & Bulgheresi, 2013).
Deep-sea invertebrates that inhabit hydrothermal vents, hydrocarbon seeps and sites of organic enrichment have evolved intriguing symbioses that compensate for the absence of sunlight and the trophic benefits of photosynthesis. Associations with chemoautotrophic bacteria that derive energy from the oxidation of sulphides, hydrogen or methane can support lush invertebrate communities in these unusual habitats (Dubilier, Bergin, & Lott, 2008). Clams of the family Vesicomyidae belong to the key fauna in chemosynthesis-based ecosystems worldwide (Johnson, Krylova, Audzijonyte, Sahling, & Vrijenhoek, 2016). Lacking a functional digestive system, they rely nutritionally on their thiotrophic gammaproteobacterial symbionts that inhabit specialized cells in the gill tissue of their host.
Previous histological and molecular studies showed that vesicomyid symbionts are maternally inherited and can be grouped into two different phylogenetic clades that differ in their status of genome reduction (Kuwahara et al., 2011). Clade I symbionts have highly reduced genomes that lack crucial genes for DNA recombination and repair, whereas clade II symbionts retain functional copies of these genes and have slightly larger genome sizes (Kuwahara et al., 2007(Kuwahara et al., , 2011Shimamura et al., 2017).
Although maternal inheritance appears to be the main transmission route of these symbionts, rare occasions of horizontal transfer have been suggested given that host mitochondrial and symbiont 16S rRNA phylogenies are sometimes incongruent Stewart, Young, & Cavanaugh, 2008Vrijenhoek, 2010). Different mechanisms have been proposed to explain how lateral acquisition of symbionts could occur in vesicomyid clams (Stewart, Young, & Cavanaugh, 2008): (a) hybridization between host species including the presence of doubly uniparental inheritance, (b) acquisition from a stable free-living symbiont population, (c) direct transfer between hosts without the involvement of hybridization, for example through contact between eggs, contact between eggs and host tissue or uptake of symbionts that have been released from moribund clams.
Two recent studies hypothesized that direct transfer is the main mechanism leading to symbiont mixtures or displacements of native symbionts in vesicomyid clams. Decker, Olu, Arnaud-Haond, and Duperron (2013) reported that individual vesicomyid clams from the Gulf of Guinea can host mixtures of native and non-native symbiont phylotypes when distinct host species co-occur in the same seep habitat. These authors argued that physical proximity could promote symbiont exchanges among very distantly related clam taxa. Ikuta et al. (2016) recently showed that the symbiont of Phreagena okutanii (previously Calyptogena okutanii) spends part of its life attached to the surface of the host's eggs, thereby strengthening Stewart et al.'s (2008) argument that direct contact between eggs or eggs and host tissues can lead to lateral symbiont transfer between co-occurring clam species. While these studies considered the possibility of host hybridization or environmental symbiont acquisition unlikely, these hypotheses have not been directly addressed by previous analyses. Here, we present a new case of discrepant symbiont compositions in two eastern Pacific clams, Archivesica gigas and Phreagena soyoae, species that are easily distinguished based on mitochondrial and nuclear markers (Johnson et al., 2016). Using demographic inference models based on genome-wide SNPs as well as traditional DNA markers, we investigated the hypothesis that hybridization between the two species might be a mechanism of horizontal symbiont transfer in this system and examined the nature of this gene flow between the taxa.

| Sample collection and preparation
Clams were collected with remotely operated vehicles (ROVs) from eight eastern Pacific seep sites during cruises between 2000 and 2015 (Table 1; Figure 1). Upon recovery of the ROV, specimens were either dissected and frozen at -80°C or preserved in 95% ethanol. DNA was extracted from symbiont-bearing gill and symbiont-free foot or adductor muscle tissue with the QIAGEN DNeasy Blood & Tissue Kit according to manufacturer's instruc- tions. An RNA digestion step was included as advised in the protocol. We constructed the map of sampling localities with ggmap in rstudio (https ://cran.r-proje ct.org/web/packa ges/ggmap/ citat ion.html).

| Sanger sequencing of host and symbiont genes
The mitochondrial cytochrome c oxidase subunit I (COI), the nuclear histone 3 (H3) and the nuclear ADP/ATP translocase (ANT) genes were used for host species identification. PCR and sequencing protocols followed Johnson et al. (2016). Assembly of forward and reverse reads, multiple alignments and phasing of nuclear genes were done as in Breusing, Johnson, Tunnicliffe, and Vrijenhoek (2015). To identify the dominant symbiont lineage in the sampled clam species, we sequenced the full-length 16S rRNA using the universal eubacterial primers 27F/1492R (Lane, 1991). Gene amplifications and sequencing reactions were performed as in Vrijenhoek, Duhaime, and Jones (2007), while sequence analysis was done as described above.

| ezRAD sequencing and estimation of allele frequencies
The Sanger sequence analyses indicated that several clams with the A. gigas COI mitotype contained the P. soyoae symbiont 16S phylotype.
To determine whether hybridization had occurred between the two host species and could therefore possibly account for the observed symbiont switch, we developed a SNP panel from ezRAD sequencing of five putatively pure A. gigas (Gorda Ridge) and five putatively pure P. soyoae (Clam Bed) individuals. These sites were chosen as references as they each contained only one clam species without any evidence for genetic admixture or symbiont discrepancies. The composition of these SNPs was evaluated in four clams from Pedro's Whalefall in which host-symbiont discrepancies were found. Library preparation and sequencing was performed at the Huntsman Cancer Institute at the University of Utah and UC Davis. The library preparation protocol was adapted from the original methods described in Toonen et al. (2013) and is provided in full detail in Appendix S1. Sequencing of the 14 clams was done with a 2 × 125-150 bp paired-end protocol on Illumina HiSeq 2500 and 4000 instruments. Following quality checks with fastqc (https ://www.bioin forma tics.babra ham.ac.uk/proje cts/fastq c/), the raw reads were compared against draft genome assemblies of the clam symbionts and host mitochondrial and ribosomal genes (C. R. Young, unpublished data) as well as the PhiX sequencing control to remove potential contaminants and obtain a purely nuclear gene data set (Appendix S2). Unmapped paired-end reads were then trimmed and quality filtered with trimmomatic (Bolger, Lohse, & Usadel, 2014) and assembled in ddocent version 2.2.17 (Puritz, Hollenbeck, & Gold, 2014) following recommendations for assembly optimization at dDocent.com. Assembly parameters were adjusted as follows: clustering threshold = 0.9; minimum coverage of a read within an individual = 6; minimum number of individuals containing a unique sequence = 4.
Basic quality metrics and information about the sequencing data are given in Appendix S3. Exhaustive exploration of various settings in the ddocent SNP filtering pipeline all resulted in spurious patterns in population-specific allele frequency spectra (AFS) and poor convergence in downstream population genomic analyses. As recently shown by Warmuth and Ellegren (2019), traditional SNP calling from RADseq data can introduce bias in demographic inference. Based on these results, we used angsd version 0.920 (Korneliussen, Albrechtsen, & Nielsen, 2014) to estimate AFSs and other population genetic statistics in this study (Appendix S4). To remove low-quality sites from the analyses, we used a minimum mapping quality of 30 (minMapQ = 30), a minimum base quality of 20 (minQ = 20) and a minimum depth of 20 (setMinDepth = 20). We further adjusted mapping quality for excessive mismatches (C = 50), removed sites with missing data, excluded spurious and improperly paired reads and computed per-base alignment qualities (BAQ = 1) to resolve false variants that were caused by misalignments. Potentially paralogous regions were excluded by discarding reads that had multiple hits to the reference assembly and by considering only sites that had a maximum depth of 250 (which we chose as reasonable threshold based on the mean read depth distribution). Inferences were based on the folded AFS due to no outgroup information available before the analysis. The joint AFS between A. gigas and P. soyoae was calculated with the angsd subprogram realsfs and subsequently folded in ∂a∂i version 1.6.3 (Gutenkunst, Hernandez, Williamson, & Bustamante, 2009).

| Phylogenetic and population genetic analyses
We used the program popart version 1.7 (http://popart.otago. ac.nz/) to create phylogenetic networks for the symbiont 16S rRNA TA B L E 1 Geographic coordinates, depths, dive numbers, sample sizes (N) and host species for the investigated clam sites

| Demographic inference
We used the program ima2 (Hey, 2010)  Complementary to the ima2 approach, the program ∂a∂i version 1.6.3 (Gutenkunst et al., 2009) was used to infer the demographic history of the two clam species from the folded joint allele frequency spectrum estimated from the RADseq data set. This approach is considerably more flexible than ima2 with respect to demographic models and the variability of rates among different genomic regions.
As in the ima2, we defined two populations based on the respective genotypic signature. We tested seven different models of evolution- was chosen as best fit. AIC weights (Burnham & Anderson, 2002;Stewart et al., 2008) were used to express relative support among the set of models that we examined.

| Sanger sequenced genes: Haplotype diversity and differentiation in hosts and symbionts
Phylogenetic networks for a 518-bp fragment of the host mito-

chondrial COI gene revealed a clear segregation of haplotypes into
A. gigas-and P. soyoae-specific clades (Figure 2

| Gene flow
Both ∂a∂i and ima2 analyses provided evidence for divergence with gene flow between the two clam species, supporting models of asymmetric migration from P. soyoae into A. gigas (Figures 4 and 5; Table 3). Despite the large phylogenetic distance between the two clam genera, the ima2 analyses could not approximate the time of population splitting or ancestral population size accurately, which indicates a lack of information to constrain these parameters due to the limited number of loci examined. ∂a∂i favoured the secondary contact model with heterogeneous gene flow (SC2M) as most likely scenario of the speciation process (Table 3;  While the SC2M model fits the data significantly better than any other model, all models with two classes of gene flow parameters (IM2M and AM2M) were better fits to the data than those without (Table 3), and better predicted the AFS observed in A. gigas and P. soyoae (Figure 5), suggesting that accounting for differential introgression rates across the genome is useful to predict certain characteristics of our data.  hybridization, (b) environmental acquisition from a free-living symbiont population or (c) host-to-host transfer, for example through direct contact between symbiont-bearing eggs or uptake of symbiont cells that have been released from a dying clam individual. Although the host-to-host transfer hypothesis has been favoured by multiple authors (Decker et al., 2013;Ikuta et al., 2016;Ozawa et al., 2017;Stewart et al., 2008), the possibility of interspecific hybridization has never been investigated.

| D ISCUSS I ON
Hybridization between different species is an important evolutionary process that can provide fundamental insights into the molecular mechanisms of reproductive isolation and adaptation.
The outcomes of interspecific gene flow can be seen as a continuum of two extremes: erosion of species barriers through merging of gene pools (Allendorf, Leary, Spruell, & Wenburg, 2001)  In this study, we examined the symbiont composition and host hybridization hypothesis in the two clam species, A. gigas and P. soyoae, from cold seep sites in the Pacific Ocean. Our barcoding analyses indicated that several individuals with the A. gigas mitotype contained the P. soyoae-specific symbiont phylotype at localities where both host species either co-occurred (Ben's Seep and Extrovert Cliff) or where A. gigas was the only taxon found (Pedro's Whalefall). Demographic inference provided evidence that asymmetric gene flow between P. soyoae and A. gigas did occur in the evolutionary history of the two species. Both ima2 and ∂a∂i analyses F I G U R E 5 Observed and fitted joint folded allele frequency spectra as calculated in ∂a∂i. The figure shows the AFS of A. gigas (x-axis, nine individuals) plotted against the AFS of P. soyoae (y-axis, five individuals). The colour scheme indicates the frequencies of minor alleles in each population across all polymorphic sites. The SC2M model was the scenario with the highest likelihood, and its AFS is shown in comparison with the other tested models favoured a model of mostly unidirectional migration from P. soyoae to A. gigas over models of strict isolation. Based on genomic information, ∂a∂i identified a secondary contact scenario with differential gene introgression as best fit to the observed joint allele frequency spectrum. The SC2M model is the most likely demographic scenario, among those that we examined, underlying the evolutionary divergence of these clam species. Interestingly, however, we did not find any evidence for admixed individuals in the investigated clam populations, as might be expected under a recent secondary contact scenario. This observation could be due to insufficient sampling. It is also possible that there have been enough generations that the genetic disequilibria from a recent hybridization event have nearly Although open questions about the demographic history of A. gigas and P. soyoae remain, our data suggest that interspecific hybridization could be a mechanism for horizontal symbiont trans-  (Bordenstein, O'Hara, & Werren, 2001;Bordenstein & Werren, 2007;Brucker, & Bordenstein, 2012Jaenike, Dyer, Cornish, & Minhas, 2006;Vala, Breeuwer, & Sabelis, 2000). Under this scenario, the P. soyoae symbiont would have an unknown fitness advantage that favours switching to a new host.
Incongruent compositions of host mitochondrial and symbiont genomes due to hybridization must involve some form of paternal cotransmission. One possibility is occasional inheritance of mitochondria and symbionts through the paternal germline (paternal leakage). Paternal leakage is a common phenomenon that occurs in a variety of different taxa at low frequency (reviewed in Breton & Stewart, 2015), which would agree with our observation that introgressed symbionts are rare. The second possibility is doubly uniparental inheritance (DUI) of mitochondria, which would involve a regulated system of paternal cotransmission as suggested previously for Vesicomya sp. mt-II (Stewart et al., 2008). Although DUI has not been described in vesicomyid clams, it is known to occur in some Veneroida (Gusman, Lecomte, Stewart, Passamonti, & Breton, 2016;Zouros, 2013). To determine whether DUI is present in vesicomyids, it would be necessary to sequence mitochondrial genomes from sexed individuals and identify if genetic differences exist between female and male mitochondria.
Despite supporting hybridization as a potential mechanism of lateral symbiont transfer, our data do not reject the alternative hy-  Abbreviations: AIC, Akaike information criterion; AICw, AIC weights; Ln, model likelihood.
and P. soyoae belong to vesicomyid symbiont Clade I which is characterized by highly reduced genomes without essential genes for an extracellular lifestyle (Kuwahara et al., 2011), so that it is improbable that a free-living population of these symbionts exists in the environment. Second, if physical proximity promoted lateral symbiont acquisition as argued by Decker et al. (2013), we would expect to find more individuals with non-native symbionts in sympatric populations of A. gigas and P. soyoae. Furthermore, non-native symbiont phylotypes were only observed in A. gigas, but not in P. soyoae.
If host-to-host transfer was the underlying mechanism leading to symbiont switching, foreign phylotypes should be observed in both species, unless natural selection acted against the P. soyoae host × A. gigas symbiont combination (as mentioned above).
Our results raise several interesting hypotheses that can be addressed in future studies. In this study, we only investigated the symbiont 16S rRNA gene. To disentangle host-to-host transfer from retention of an introgressed strain, it will be necessary to sequence and compare whole genomes of symbionts from clam individuals that contain the native and foreign P. soyoae phylotypes. A symbiont that was transferred historically via host hybridization can be expected to be highly divergent from the native phylotype, while a symbiont that was transmitted via contemporary host-to-host transfer should be relatively similar. Since we sequenced the 16S rRNA gene directly, we were not able to uncover potential symbiont mixtures in individual host animals, given that this sequencing approach is biased towards the most abundant phylotype (Zimmermann et al., 2014).
Furthermore, symbiont types could be variable across the host gill, as seen in other taxa from chemosynthetic environments (Duperron et al., 2005;Zimmermann et al., 2014). Studies that involved cloning and pyrosequencing techniques showed that in some cases divergent symbiont lineages can co-occur in a single clam host (Decker et al., 2013;). To examine whether and how different symbiont lineages coexist in A. gigas and P. soyoae, whole-genome analyses based on high-throughput sequencing techniques will be useful. Comparative genomic and population genomic approaches will help to illuminate the genomic consequences of occasional lateral symbiont acquisition in deep-sea vesicomyid clams and lead to a better understanding of how host-symbiont interactions shape the evolution of both partners.

ACK N OWLED G EM ENTS
We thank the crews of the RV Point Lobos, Rachel Carson and Western Flyer and pilots of the ROV Ventana, Tiburon and Doc Ricketts for carefully collecting innumerable clams. We would also like to thank Jim Barry for letting us raid his freezers to collect even more clams. Thank you as well to Shana Goffredi who was responsible for many of the early clam collections and also provided useful discussions about symbiosis. We also want to thank Brian Dalley who developed and tested the ezRAD library preparation protocol for us. Funding was provided by grants of the German Research  Bioinformatic code for the ezRAD analyses is provided in the Appendix S1-S6.