# Abundant toxin-related genes in the genomes of beneficial symbionts from deep-sea hydrothermal vent mussels

Lizbeth Sayavedra, Manuel Kleiner, Ruby Ponnudurai, Silke Wetzel, Eric Pelletier, Valérie Barbe, Nori Satoh, Eiichi Shoguchi, Dennis Fink, Corinna Breusing, et al.

<figure><img src="image_1.png" /></figure>

## To cite this version:

Lizbeth Sayavedra, Manuel Kleiner, Ruby Ponnudurai, Silke Wetzel, Eric Pelletier, et al. Abundant toxin-related genes in the genomes of beneficial symbionts from deep-sea hydrothermal vent mussels. eLife, 2015, 4, pp.e07966. 10.7554/eLife.07966. cea-04573299

HAL Id: cea-04573299

https://cea.hal.science/cea-04573299

Submitted on 13 May 2024

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>
# Abundant toxin-related genes in the genomes of beneficial symbionts from deep-sea hydrothermal vent mussels

Lizbeth Sayavedra, Manuel Kleiner, Ruby Ponnudurai, Silke Wetzel, Eric Pelletier, Valérie Barbe, Nori Satoh, Eiichi Shoguchi, Dennis Fink, Corinna Breusing, et al.

<figure><img src="image_1.png" /></figure>

## To cite this version:

Lizbeth Sayavedra, Manuel Kleiner, Ruby Ponnudurai, Silke Wetzel, Eric Pelletier, et al. Abundant toxin-related genes in the genomes of beneficial symbionts from deep-sea hydrothermal vent mussels. eLife, 2015, 4, pp.e07966. 10.7554/eLife.07966. cea-04573299

HAL Id: cea-04573299

https://cea.hal.science/cea-04573299v1

Submitted on 13 May 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
<figure><img src="image_3.png" /></figure>

# Abundant toxin-related genes in the genomes of beneficial symbionts from deep-sea hydrothermal vent mussels

Lizbeth Sayavedra¹, Manuel Kleiner¹†, Ruby Ponnudurai²†, Silke Wetzel¹, Eric Pelletier³,⁴,¹¹, Valerie Barbe³, Nori Satoh⁵, Eiichi Shoguchi⁵, Dennis Fink¹, Corinna Breusing⁶, Thorsten BH Reusch⁶, Philip Rosenstiel⁷, Markus B Schilhabel⁷, Dörte Becher⁸,⁹, Thomas Schweder²,⁸, Stephanie Markert²,⁸, Nicole Dubilier¹,¹⁰, Jillian M Petersen¹*‡

¹Max Planck Institute for Marine Microbiology, Bremen, Germany; ²Institute of Pharmacy, Ernst-Moritz-Arndt-University, Greifswald, Germany; ³Genoscope - Centre National de Séquençage, Commissariat à l'énergie atomique et aux énergies alternatives, Evry, France; ⁴Metabolic Genomics Group, Commissariat à l'énergie atomique et aux énergies alternatives, Evry, France; ⁵Marine Genomics Unit, Okinawa Institute of Science and Technology, Onna, Japan; ⁶Evolutionary Ecology, GEOMAR Helmholtz Centre for Ocean Research Kiel, Kiel, Germany; ⁷Institute of Clinical Molecular Biology, Kiel, Germany; ⁸Institute of Marine Biotechnology, Greifswald, Germany; ⁹Institute of Microbiology, Ernst-Moritz-Arndt-University, Greifswald, Germany; ¹⁰University of Bremen, Bremen, Germany; ¹¹University of Évry-Val d'Essonne, Evry, France

**Present address:** ‡Department of Microbiology and Ecosystem Science, Division of Microbial Ecology, Research Network Chemistry Meets Microbiology, University of Vienna, Vienna, Austria

**Competing interests:** The authors declare that no competing interests exist.

**Funding:** See page 22

**Received:** 09 April 2015

**Accepted:** 14 September 2015

**Published:** 15 September 2015

**Reviewing editor:** Axel A Brakhage, Friedrich Schiller University Jena and Hans-Knöll-Institut, Germany

© Copyright Sayavedra et al. This article is distributed under the terms of the [Creative Commons Attribution License](http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use and redistribution provided that the original author and source are credited.

**Abstract** *Bathymodiolus* mussels live in symbiosis with intracellular sulfur-oxidizing (SOX) bacteria that provide them with nutrition. We sequenced the SOX symbiont genomes from two *Bathymodiolus* species. Comparison of these symbiont genomes with those of their closest relatives revealed that the symbionts have undergone genome rearrangements, and up to 35% of their genes may have been acquired by horizontal gene transfer. Many of the genes specific to the symbionts were homologs of virulence genes. We discovered an abundant and diverse array of genes similar to insecticidal toxins of nematode and aphid symbionts, and toxins of pathogens such as *Yersinia* and *Vibrio*. Transcriptomics and proteomics revealed that the SOX symbionts express the toxin-related genes (TRGs) in their hosts. We hypothesize that the symbionts use these TRGs in beneficial interactions with their host, including protection against parasites. This would explain why a mutualistic symbiont would contain such a remarkable 'arsenal' of TRGs.

DOI: 10.7554/eLife.07966.001

## Introduction

Mussels of the genus *Bathymodiolus* dominate deep-sea hydrothermal vents and cold seeps worldwide. The key to their ecological and evolutionary success is their symbiosis with chemosynthetic bacteria that provide them with nutrition (von Cosel, 2002; Van Dover et al., 2002). *Bathymodiolus* mussels host their symbionts inside specialized gill epithelial cells called bacteriocytes (Cavanaugh et al., 2006; Petersen and Dubilier, 2009).

Their filtering activity exposes *Bathymodiolus* mussels to a plethora of diverse microbes in their environment. Despite this, they are colonized by only one or a few specific types of chemosynthetic symbionts. Some mussel species associate exclusively with sulfur-oxidizing (SOX) symbionts that use

<footnote>*For correspondence: jmpeters@mpi-bremen.de</footnote>

<footnote>†These authors contributed equally to this work</footnote>
<figure><img src="image_3.png" /></figure>

# Abundant toxin-related genes in the genomes of beneficial symbionts from deep-sea hydrothermal vent mussels

Lizbeth Sayavedra¹, Manuel Kleiner¹†, Ruby Ponnudurai²†, Silke Wetzel¹, Eric Pelletier³,⁴,¹¹, Valerie Barbe³, Nori Satoh⁵, Eiichi Shoguchi⁵, Dennis Fink¹, Corinna Breusing⁶, Thorsten BH Reusch⁶, Philip Rosenstiel⁷, Markus B Schilhabel⁷, Dörte Becher⁸,⁹, Thomas Schweder²,⁸, Stephanie Markert²,⁸, Nicole Dubilier¹,¹⁰, Jillian M Petersen¹*‡

¹Max Planck Institute for Marine Microbiology, Bremen, Germany; ²Institute of Pharmacy, Ernst-Moritz-Arndt-University, Greifswald, Germany; ³Genoscope - Centre National de Séquençage, Commissariat à l'énergie atomique et aux énergies alternatives, Evry, France; ⁴Metabolic Genomics Group, Commissariat à l'énergie atomique et aux énergies alternatives, Evry, France; ⁵Marine Genomics Unit, Okinawa Institute of Science and Technology, Onna, Japan; ⁶Evolutionary Ecology, GEOMAR Helmholtz Centre for Ocean Research Kiel, Kiel, Germany; ⁷Institute of Clinical Molecular Biology, Kiel, Germany; ⁸Institute of Marine Biotechnology, Greifswald, Germany; ⁹Institute of Microbiology, Ernst-Moritz-Arndt-University, Greifswald, Germany; ¹⁰University of Bremen, Bremen, Germany; ¹¹University of Évry-Val d'Essonne, Evry, France

**Present address:** ‡Department of Microbiology and Ecosystem Science, Division of Microbial Ecology, Research Network Chemistry Meets Microbiology, University of Vienna, Vienna, Austria

**Competing interests:** The authors declare that no competing interests exist.

**Funding:** See page 22

**Received:** 09 April 2015

**Accepted:** 14 September 2015

**Published:** 15 September 2015

**Reviewing editor:** Axel A Brakhage, Friedrich Schiller University Jena and Hans-Knöll-Institut, Germany

© Copyright Sayavedra et al. This article is distributed under the terms of the [Creative Commons Attribution License](http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use and redistribution provided that the original author and source are credited.

**Abstract** *Bathymodiolus* mussels live in symbiosis with intracellular sulfur-oxidizing (SOX) bacteria that provide them with nutrition. We sequenced the SOX symbiont genomes from two *Bathymodiolus* species. Comparison of these symbiont genomes with those of their closest relatives revealed that the symbionts have undergone genome rearrangements, and up to 35% of their genes may have been acquired by horizontal gene transfer. Many of the genes specific to the symbionts were homologs of virulence genes. We discovered an abundant and diverse array of genes similar to insecticidal toxins of nematode and aphid symbionts, and toxins of pathogens such as *Yersinia* and *Vibrio*. Transcriptomics and proteomics revealed that the SOX symbionts express the toxin-related genes (TRGs) in their hosts. We hypothesize that the symbionts use these TRGs in beneficial interactions with their host, including protection against parasites. This would explain why a mutualistic symbiont would contain such a remarkable 'arsenal' of TRGs.

DOI: 10.7554/eLife.07966.001

## Introduction

Mussels of the genus *Bathymodiolus* dominate deep-sea hydrothermal vents and cold seeps worldwide. The key to their ecological and evolutionary success is their symbiosis with chemosynthetic bacteria that provide them with nutrition (von Cosel, 2002; Van Dover et al., 2002). *Bathymodiolus* mussels host their symbionts inside specialized gill epithelial cells called bacteriocytes (Cavanaugh et al., 2006; Petersen and Dubilier, 2009).

Their filtering activity exposes *Bathymodiolus* mussels to a plethora of diverse microbes in their environment. Despite this, they are colonized by only one or a few specific types of chemosynthetic symbionts. Some mussel species associate exclusively with sulfur-oxidizing (SOX) symbionts that use

<footnote>*For correspondence: jmpeters@mpi-bremen.de</footnote>

<footnote>†These authors contributed equally to this work</footnote>
**eLife digest** Although bacteria are commonly associated with causing illness, many are actually beneficial to the organism they live in or on. The phenomenon of one species helping another to survive is known as symbiosis.

Animals thrive at hydrothermal vents in the deep sea because of their partnerships with symbiotic bacteria. The bacteria use the geochemical energy found at hydrothermal vents to convert carbon into sugars, thus providing their animal hosts with essential nutrients. Unlike the symbiotic communities that associate with humans and other mammals, in which thousands of bacterial species co-exist, deep-sea mussels associate with just one or two species of symbiotic bacteria. This relative simplicity is ideal for investigating how the intimate associations between animals and bacteria work.

Genes contain the instructions cells and organisms need to survive, and so one way that researchers investigate symbiosis is by studying the genes of the organisms involved. Such studies of beneficial bacteria are beginning to reveal that the molecular mechanisms involved in symbiosis are remarkably similar to those responsible for the harmful effects produced by some bacteria.

By performing genetic sequencing on the symbiotic bacteria from deep-sea mussels, Sayavedra et al. have discovered that the bacteria have an unusually large number of toxin-like genes, and that all of these genes are active in the bacteria when they are inside host mussels. This was unexpected, as the bacteria are known to benefit their mussel hosts. The toxin-like genes from the symbiotic bacteria are similar to toxins found in the bacteria that cause diseases such as cholera and the plague in humans and other animals.

Sayavedra et al. suggest that the symbiotic bacteria have 'tamed' these toxins to use them in beneficial interactions with their host. For example, some of the toxins could help the bacteria and mussels to recognize and interact with each other, and others could help to protect the mussel host from its natural enemies. The next step will be to test these ideas, which will be challenging as the mussels cannot be bred in the laboratory.

DOI: 10.7554/eLife.07966.002

reduced sulfur compounds and sometimes hydrogen as an energy source, and carbon dioxide as a carbon source. Some have only methane-oxidizing (MOX) symbionts that use methane as an energy source and carbon source. Some mussel species host both types in a dual symbiosis (Fisher et al., 1993; Distel et al., 1995; Duperron et al., 2006; Dubilier et al., 2008; Petersen et al., 2011). In all species except one, a single 16S rRNA phylotype for each type of symbiont (SOX or MOX) is found in the gills (Dubilier et al., 2008). There are more than 30 described *Bathymodiolus* species, and most associate with a characteristic symbiont phylotype, which is not found in other species (Duperron et al., 2013).

Although these associations are clearly very specific, the molecular mechanisms that underpin this specificity are still unknown. No chemosynthetic symbiont has ever been obtained in pure culture. Therefore, molecular methods for investigating uncultured microbes have been essential for understanding their biodiversity, function, and evolution (reviewed by Dubilier et al., 2008).

The *Bathymodiolus* symbionts are assumed to be horizontally transmitted, which means that each new host generation must take up their symbionts from the surrounding environment or co-occurring adults (Won et al., 2003b; Kadar et al., 2005; DeChaine et al., 2006; Fontanez and Cavanaugh, 2014; Wentrup et al., 2014). To initiate the symbiosis, hosts and symbionts must have evolved highly specific recognition and attachment mechanisms. Once they have been recognized, the symbionts need to enter host cells and avoid immediate digestion, just like other intracellular symbionts such as *Burkholderia rhizoxinica* and *Rhizobium leguminosarum*, or pathogens such as *Legionella*, *Listeria*, or *Yersinia* (Hentschel et al., 2000; Moebius et al., 2014). Indeed, like many intracellular pathogens, the *Bathymodiolus* symbionts seem to induce a loss of microvilli on the cells they colonize (*Cossart and Sansonetti, 2004; Bhavsar et al., 2007; Haglund and Welch, 2011; Wentrup et al., 2014). Finally, the symbionts achieve dense populations inside the host cells (e.g., *Duperron et al., 2006; Halary et al., 2008*). Therefore, they must be able to avoid immediate digestion by their hosts. Although the mechanisms of host cell entry and immune evasion have been extensively studied in pathogens and plant-microbe associations such as the rhizobia-legume symbiosis, far less is known about the mechanisms beneficial symbionts use to enter and survive within animal host cells.
**eLife digest** Although bacteria are commonly associated with causing illness, many are actually beneficial to the organism they live in or on. The phenomenon of one species helping another to survive is known as symbiosis.

Animals thrive at hydrothermal vents in the deep sea because of their partnerships with symbiotic bacteria. The bacteria use the geochemical energy found at hydrothermal vents to convert carbon into sugars, thus providing their animal hosts with essential nutrients. Unlike the symbiotic communities that associate with humans and other mammals, in which thousands of bacterial species co-exist, deep-sea mussels associate with just one or two species of symbiotic bacteria. This relative simplicity is ideal for investigating how the intimate associations between animals and bacteria work.

Genes contain the instructions cells and organisms need to survive, and so one way that researchers investigate symbiosis is by studying the genes of the organisms involved. Such studies of beneficial bacteria are beginning to reveal that the molecular mechanisms involved in symbiosis are remarkably similar to those responsible for the harmful effects produced by some bacteria.

By performing genetic sequencing on the symbiotic bacteria from deep-sea mussels, Sayavedra et al. have discovered that the bacteria have an unusually large number of toxin-like genes, and that all of these genes are active in the bacteria when they are inside host mussels. This was unexpected, as the bacteria are known to benefit their mussel hosts. The toxin-like genes from the symbiotic bacteria are similar to toxins found in the bacteria that cause diseases such as cholera and the plague in humans and other animals.

Sayavedra et al. suggest that the symbiotic bacteria have 'tamed' these toxins to use them in beneficial interactions with their host. For example, some of the toxins could help the bacteria and mussels to recognize and interact with each other, and others could help to protect the mussel host from its natural enemies. The next step will be to test these ideas, which will be challenging as the mussels cannot be bred in the laboratory.

DOI: 10.7554/eLife.07966.002

reduced sulfur compounds and sometimes hydrogen as an energy source, and carbon dioxide as a carbon source. Some have only methane-oxidizing (MOX) symbionts that use methane as an energy source and carbon source. Some mussel species host both types in a dual symbiosis (Fisher et al., 1993; Distel et al., 1995; Duperron et al., 2006; Dubilier et al., 2008; Petersen et al., 2011). In all species except one, a single 16S rRNA phylotype for each type of symbiont (SOX or MOX) is found in the gills (Dubilier et al., 2008). There are more than 30 described *Bathymodiolus* species, and most associate with a characteristic symbiont phylotype, which is not found in other species (Duperron et al., 2013).

Although these associations are clearly very specific, the molecular mechanisms that underpin this specificity are still unknown. No chemosynthetic symbiont has ever been obtained in pure culture. Therefore, molecular methods for investigating uncultured microbes have been essential for understanding their biodiversity, function, and evolution (reviewed by Dubilier et al., 2008).

The *Bathymodiolus* symbionts are assumed to be horizontally transmitted, which means that each new host generation must take up their symbionts from the surrounding environment or co-occurring adults (Won et al., 2003b; Kadar et al., 2005; DeChaine et al., 2006; Fontanez and Cavanaugh, 2014; Wentrup et al., 2014). To initiate the symbiosis, hosts and symbionts must have evolved highly specific recognition and attachment mechanisms. Once they have been recognized, the symbionts need to enter host cells and avoid immediate digestion, just like other intracellular symbionts such as *Burkholderia rhizoxinica* and *Rhizobium leguminosarum*, or pathogens such as *Legionella*, *Listeria*, or *Yersinia* (Hentschel et al., 2000; Moebius et al., 2014). Indeed, like many intracellular pathogens, the *Bathymodiolus* symbionts seem to induce a loss of microvilli on the cells they colonize (*Cossart and Sansonetti, 2004; Bhavsar et al., 2007; Haglund and Welch, 2011; Wentrup et al., 2014). Finally, the symbionts achieve dense populations inside the host cells (e.g., *Duperron et al., 2006; Halary et al., 2008*). Therefore, they must be able to avoid immediate digestion by their hosts. Although the mechanisms of host cell entry and immune evasion have been extensively studied in pathogens and plant-microbe associations such as the rhizobia-legume symbiosis, far less is known about the mechanisms beneficial symbionts use to enter and survive within animal host cells.
The symbiosis between *Vibrio fisheri* bacteria and *Euprymna scolopes* squid is one of the few beneficial host-microbe associations where the molecular mechanisms of host-symbiont interaction have been investigated. A number of factors are involved in initiating this symbiosis such as the symbiont-encoded 'TCT toxin', which is related to the tracheal cytotoxin of *Bordetella pertussis* (*McFall-Ngai* et al., 2013). A few studies of intracellular insect symbionts have shown that they use type III and type IV secretion systems to establish and maintain their association with their host (reviewed by *Dale and Moran, 2006; Snyder and Rio, 2013*). These secretion systems are commonly used by intracellular pathogens to hijack host cell processes, allowing their entry and survival within host cells (e.g., *Hueck, 1998; Steele-Mortimer et al., 2002; Backert and Meyer, 2006*). An example is the *Sodalis* symbionts of aphids and weevils, which use a type III secretion system for entry to the host cell and are thought to have evolved from pathogens (*Dale et al., 2001; Clayton et al., 2012*). The virulence determinants of their pathogenic ancestors might therefore have been co-opted for use in beneficial interactions with their insect hosts.

In contrast to the *Sodalis* symbionts of insects and the *Vibrio* symbionts of squid, the *Bathymodiolus* SOX symbionts are not closely related to any known pathogens. Moreover, because they fall interspersed between free-living SOX bacteria in 16S rRNA phylogenies, they are hypothesized to have evolved multiple times from free-living ancestors (Figure 2—figure supplement 1) (*Petersen et al., 2012*). Comparative genomics is a powerful tool for identifying the genomic basis of beneficial and pathogenic interactions, particularly if the symbionts or pathogens have close free-living relatives that do not associate with a host (e.g., *Galagan, 2014; Ogier et al., 2014; Zuleta et al., 2014*). Genomes of closely related free-living and symbiotic relatives of *Bathymodiolus* SOX symbionts were recently published. Their closest free-living relatives are marine SOX bacteria called SUP05, which are abundant in the world's oceans, particularly in oxygen minimum zones (OMZs) and hydrothermal plumes (*Sunamura et al., 2004; Lavik et al., 2009; Walsh et al., 2009; Anantharaman et al., 2012; Wright et al., 2012; Petersen and Dubilier, 2014*). The *Bathymodiolus* SOX symbionts and SUP05 bacteria form a monophyletic clade together with the SOX symbionts of vesicomyid clams based on 16S rRNA gene phylogenies (Figure 2—figure supplement 1) (*Distel et al., 1995; Petersen et al., 2012*). Closed genomes are available for the symbionts of two clam species (*Kuwahara et al., 2007; Newton et al., 2007*).

All members of this monophyletic group (the mussel and clam symbionts, and SUP05) share similar core metabolic features. They are all capable of autotrophic growth, and all use reduced sulfur compounds as an energy source (*Newton et al., 2008; Walsh et al., 2009*). They can differ in auxiliary metabolic capabilities such as hydrogen oxidation, nitrate reduction, or mixotrophy (*Newton et al., 2008; Petersen et al., 2011; Anantharaman et al., 2012; Murillo et al., 2014*). However, the major difference between these organisms is their lifestyle: SUP05 bacteria are exclusively free-living. The clam symbionts are exclusively host-associated, are vertically transmitted, and have reduced genomes. The *Bathymodiolus* symbionts appear to have adapted to both niches, as they have a host-associated stage and are assumed to also have a free-living stage.

The goal of this study was to identify the genomic basis of host-symbiont interactions in *Bathymodiolus* symbioses. We used high-throughput sequencing and binning techniques to assemble the first essentially complete draft genomes of the SOX symbionts from *Bathymodiolus* mussels. We used comparative genomics of the symbionts' genomes to those of their close free-living and obligate symbiotic relatives to reveal genes potentially involved in *Bathymodiolus* host-symbiont interactions. We used phylogenetics and bioinformatic prediction of horizontally acquired genes to investigate the origins of these genes. Finally, we used transcriptomics and proteomics to determine whether potential host-symbiont interaction genes are being expressed by the symbionts in their host.

## Results

### Draft genome sequences of *Bathymodiolus* symbionts

We sequenced the genomes of the SOX symbionts from three *Bathymodiolus* individuals: two were *Bathymodiolus azoricus* from the Menez Gwen vent field on the northern Mid-Atlantic Ridge (MAR) (Figure 1). We refer to these as BazSymA and BazSymB. The third mussel individual was an undescribed *Bathymodiolus* species (BspSym), from the Lilliput hydrothermal vent on the southern MAR (SMAR) (Figure 1). Symbiont draft genomes from each individual were almost complete (see 'Materials and methods'). Despite different sequencing and assembly strategies, the draft genomes
The symbiosis between *Vibrio fisheri* bacteria and *Euprymna scolopes* squid is one of the few beneficial host-microbe associations where the molecular mechanisms of host-symbiont interaction have been investigated. A number of factors are involved in initiating this symbiosis such as the symbiont-encoded ‘TCT toxin’, which is related to the tracheal cytotoxin of *Bordetella pertussis* (*McFall-Ngai* et al., 2013). A few studies of intracellular insect symbionts have shown that they use type III and type IV secretion systems to establish and maintain their association with their host (reviewed by *Dale and Moran, 2006*; *Snyder and Rio, 2013*). These secretion systems are commonly used by intracellular pathogens to hijack host cell processes, allowing their entry and survival within host cells (e.g., *Hueck, 1998*; *Steele-Mortimer et al., 2002*; *Backert and Meyer, 2006*). An example is the *Sodalis* symbionts of aphids and weevils, which use a type III secretion system for entry to the host cell and are thought to have evolved from pathogens (*Dale et al., 2001*; *Clayton et al., 2012*). The virulence determinants of their pathogenic ancestors might therefore have been co-opted for use in beneficial interactions with their insect hosts.

In contrast to the *Sodalis* symbionts of insects and the *Vibrio* symbionts of squid, the *Bathymodiolus* SOX symbionts are not closely related to any known pathogens. Moreover, because they fall interspersed between free-living SOX bacteria in 16S rRNA phylogenies, they are hypothesized to have evolved multiple times from free-living ancestors (Figure 2—figure supplement 1) (*Petersen et al., 2012*). Comparative genomics is a powerful tool for identifying the genomic basis of beneficial and pathogenic interactions, particularly if the symbionts or pathogens have close free-living relatives that do not associate with a host (e.g., *Galagan, 2014*; *Ogier et al., 2014*; *Zuleta et al., 2014*). Genomes of closely related free-living and symbiotic relatives of *Bathymodiolus* SOX symbionts were recently published. Their closest free-living relatives are marine SOX bacteria called SUP05, which are abundant in the world’s oceans, particularly in oxygen minimum zones (OMZs) and hydrothermal plumes (*Sunamura et al., 2004*; *Lavik et al., 2009*; *Walsh et al., 2009*; *Anantharaman et al., 2012*; *Wright et al., 2012*; *Petersen and Dubilier, 2014*). The *Bathymodiolus* SOX symbionts and SUP05 bacteria form a monophyletic clade together with the SOX symbionts of vesicomyid clams based on 16S rRNA gene phylogenies (Figure 2—figure supplement 1) (*Distel et al., 1995*; *Petersen et al., 2012*). Closed genomes are available for the symbionts of two clam species (*Kuwahara et al., 2007*; *Newton et al., 2007*).

All members of this monophyletic group (the mussel and clam symbionts, and SUP05) share similar core metabolic features. They are all capable of autotrophic growth, and all use reduced sulfur compounds as an energy source (*Newton et al., 2008*; *Walsh et al., 2009*). They can differ in auxiliary metabolic capabilities such as hydrogen oxidation, nitrate reduction, or mixotrophy (*Newton et al., 2008*; *Petersen et al., 2011*; *Anantharaman et al., 2012*; *Murillo et al., 2014*). However, the major difference between these organisms is their lifestyle: SUP05 bacteria are exclusively free-living. The clam symbionts are exclusively host-associated, are vertically transmitted, and have reduced genomes. The *Bathymodiolus* symbionts appear to have adapted to both niches, as they have a host-associated stage and are assumed to also have a free-living stage.

The goal of this study was to identify the genomic basis of host-symbiont interactions in *Bathymodiolus* symbioses. We used high-throughput sequencing and binning techniques to assemble the first essentially complete draft genomes of the SOX symbionts from *Bathymodiolus* mussels. We used comparative genomics of the symbionts’ genomes to those of their close free-living and obligate symbiotic relatives to reveal genes potentially involved in *Bathymodiolus* host-symbiont interactions. We used phylogenetics and bioinformatic prediction of horizontally acquired genes to investigate the origins of these genes. Finally, we used transcriptomics and proteomics to determine whether potential host-symbiont interaction genes are being expressed by the symbionts in their host.

## Results

### Draft genome sequences of *Bathymodiolus* symbionts

We sequenced the genomes of the SOX symbionts from three *Bathymodiolus* individuals: two were *Bathymodiolus azoricus* from the Menez Gwen vent field on the northern Mid-Atlantic Ridge (MAR) (Figure 1). We refer to these as BazSymA and BazSymB. The third mussel individual was an undescribed *Bathymodiolus* species (BspSym), from the Lilliput hydrothermal vent on the southern MAR (SMAR) (Figure 1). Symbiont draft genomes from each individual were almost complete (see ‘Materials and methods’). Despite different sequencing and assembly strategies, the draft genomes
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Sampling sites. Map showing the sampling sites of *Bathymodiolus* mussels at hydrothermal vents along the Mid-Atlantic Ridge (red stars). *B. sp.* is found at Lilliput (BspSym), *Bathymodiolus azoricus* at Menez Gwen (BazSymA and BazSymB) and Lucky Strike. The details of the sampling sites are described in *Supplementary file 1E*. The map was produced with GeoMapApp 3.3.</figcaption></figure>

DOI: 10.7554/eLife.07966.003

were 90.7–97.7% complete (Table 1). The total assembly sizes were between 1.7 and 2.3 Mbp, on 52 to 506 contigs (Table 1). Each draft genome contained one copy of the 16S rRNA gene. The BazSymB assembly only contained an 829 bp fragment of the 16S rRNA gene; however, we PCR amplified and sequenced this gene from the DNA used to generate the metagenome. The 16S rRNA genes from the two *B. azoricus* symbionts were 100% identical and were 99.3% identical to BspSym. The core metabolic potential of the *Bathymodiolus* SOX symbionts is described in Appendix 1-Symbiont metabolism. A detailed description of the genomes is beyond the scope of this article and will be published elsewhere.

**General genome comparison**

*Bathymodiolus* symbiont genomes were more similar to the free-living SUP05 than to the clam symbionts in terms of size and GC content (Table 1). Analysis of codon usage showed that all three *Bathymodiolus* SOX symbiont genomes had a greater proportion of genes that may have been acquired through recent horizontal gene transfer (HGT) compared to the clam symbionts *Candidatus Vesicomyosocius okutanii*, and *Candidatus Ruthia magnifica* (Table 1). BspSym, the genome that assembled into the fewest contigs, lacked synteny compared to SUP05 and the clam symbionts, as shown by whole genome alignment (Figure 2—figure supplement 2). Genome alignment of BazSymA and BazSymB was not attempted because the assemblies were highly fragmented. The possibility of incorrect genome assembly for BspSym was ruled out for four regions by PCR amplification of sequences spanning the regions without synteny. For confirmation, one PCR product was Sanger-sequenced and found to be identical to the draft genome assembly of BspSym. These regions without genome synteny therefore most likely represent true genome reshuffling in *Bathymodiolus* symbionts.

The *Bathymodiolus* symbiont genomes had more mobile elements compared to their closest relatives (*Supplementary file 1A*). *Bathymodiolus* symbionts had between 13 and 23 transposases and three to five integrases. SUP05 had 14 transposases and one integrase. We did not find any
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Sampling sites. Map showing the sampling sites of *Bathymodiolus* mussels at hydrothermal vents along the Mid-Atlantic Ridge (red stars). *B. sp.* is found at Lilliput (BspSym), *Bathymodiolus azoricus* at Menez Gwen (BazSymA and BazSymB) and Lucky Strike. The details of the sampling sites are described in *Supplementary file 1E*. The map was produced with GeoMapApp 3.3.

DOI: 10.7554/eLife.07966.003</figcaption></figure>

were 90.7–97.7% complete (Table 1). The total assembly sizes were between 1.7 and 2.3 Mbp, on 52 to 506 contigs (Table 1). Each draft genome contained one copy of the 16S rRNA gene. The BazSymB assembly only contained an 829 bp fragment of the 16S rRNA gene; however, we PCR amplified and sequenced this gene from the DNA used to generate the metagenome. The 16S rRNA genes from the two *B. azoricus* symbionts were 100% identical and were 99.3% identical to BspSym. The core metabolic potential of the *Bathymodiolus* SOX symbionts is described in Appendix 1-Symbiont metabolism. A detailed description of the genomes is beyond the scope of this article and will be published elsewhere.

**General genome comparison**

*Bathymodiolus* symbiont genomes were more similar to the free-living SUP05 than to the clam symbionts in terms of size and GC content (Table 1). Analysis of codon usage showed that all three *Bathymodiolus* SOX symbiont genomes had a greater proportion of genes that may have been acquired through recent horizontal gene transfer (HGT) compared to the clam symbionts *Candidatus Vesicomyosocius okutanii*, and *Candidatus Ruthia magnifica* (Table 1). BspSym, the genome that assembled into the fewest contigs, lacked synteny compared to SUP05 and the clam symbionts, as shown by whole genome alignment (Figure 2—figure supplement 2). Genome alignment of BazSymA and BazSymB was not attempted because the assemblies were highly fragmented. The possibility of incorrect genome assembly for BspSym was ruled out for four regions by PCR amplification of sequences spanning the regions without synteny. For confirmation, one PCR product was Sanger-sequenced and found to be identical to the draft genome assembly of BspSym. These regions without genome synteny therefore most likely represent true genome reshuffling in *Bathymodiolus* symbionts.

The *Bathymodiolus* symbiont genomes had more mobile elements compared to their closest relatives (*Supplementary file 1A*). *Bathymodiolus* symbionts had between 13 and 23 transposases and three to five integrases. SUP05 had 14 transposases and one integrase. We did not find any
<table><caption>**Table 1.** Overview of the genomes compared in this study: SOX symbiont of B. sp. two individual SOX symbionts of B. azoricus, SOX symbiont Candidatus Vesicomyosocious okutanii, SOX symbiont of Calyptogena magnifica (Candidatus Ruthia magnifica), and free-living SUP05</caption><thead><tr><th>Genome</th><th>Collection site</th><th>Contigs</th><th>GC content (%)</th><th>Length/Span (Mb)<sup>a</sup></th><th>Number of CDSs</th><th>HGT</th><th>Estimated completeness<sup>‡</sup></th><th>Coverage<sup>§</sup></th><th>Separation method<sup>#</sup></th><th>References</th></tr></thead><tbody><tr><td>B. sp. symbiont (BspSym)</td><td>Lilliput</td><td>52</td><td>38.23</td><td>1.8/2.3</td><td>2225</td><td>33%</td><td>95.39%</td><td>22X</td><td>Filtration</td><td>Petersen et al., 2011, this study</td></tr><tr><td>B. azoricus symbiont (BazSymB)</td><td>Menez Gwen</td><td>239</td><td>38.20</td><td>1.5/1.7</td><td>1802</td><td>30%</td><td>90.60%</td><td>8X</td><td>Gradient centrifugation/ binning</td><td>This study</td></tr><tr><td>B. azoricus symbiont (BazSymA)*</td><td>Menez Gwen</td><td>506</td><td>37.58</td><td>1.85/1.85</td><td>2008</td><td>35%</td><td>97.70%</td><td>59X</td><td>Binning</td><td>This study</td></tr><tr><td>Ca. V. okutanii</td><td>Sagami Bay</td><td>1</td><td>31.59</td><td>1.0/1.0</td><td>980</td><td>26%</td><td>93.58%</td><td>-</td><td>Whole genome assembly</td><td>Kuwahara et al., 2007</td></tr><tr><td>Ca. R. magnifica</td><td>East Pacific Rise, 9°N</td><td>1</td><td>34.03</td><td>1.2/1.2</td><td>1210</td><td>23%</td><td>94.84%</td><td>-</td><td>Whole genome assembly</td><td>Newton et al., 2007</td></tr><tr><td>SUP05</td><td>Saanich Inlet</td><td>97</td><td>39.29</td><td>1.4/2.5</td><td>1586</td><td>30%</td><td>85.76%</td><td>-</td><td>Binning</td><td>Walsh et al., 2009</td></tr></tbody></table>

SOX, sulfur-oxidizing.

*SOX symbiont sequences recovered from metagenome of adductor muscle.

HGT = Genes that potentially originated from horizontal gene transfer.

†Length is the total length of sequence information on contigs without Ns, and span is the entire length of scaffold assembly including Ns.

‡The completeness of the genome was estimated with CheckM using a set of lineage-specific genes for proteobacteria (Parks et al., 2015).

§Median coverage.

\#Separation method indicates the experimental separation of symbionts from host tissue and co-occurring symbionts (filtration or gradient centrifugation), or the in silico separation of genomic information from hosts and co-occurring bacteria (binning).

DOI: 10.7554/eLife.07966.004
<table><caption>**Table 1.** Overview of the genomes compared in this study: SOX symbiont of B. sp. two individual SOX symbionts of B. azoricus, SOX symbiont Candidatus Vesicomyosocious okutanii, SOX symbiont of Calyptogena magnifica (Candidatus Ruthia magnifica), and free-living SUP05</caption><thead><tr><th>Genome</th><th>Collection site</th><th>Contigs</th><th>GC content (%)</th><th>Length/Span (Mb)<sup>a</sup></th><th>Number of CDSs</th><th>HGT</th><th>Estimated completeness<sup>‡</sup></th><th>Coverage<sup>§</sup></th><th>Separation method<sup>#</sup></th><th>References</th></tr></thead><tbody><tr><td>B. sp. symbiont (BspSym)</td><td>Lilliput</td><td>52</td><td>38.23</td><td>1.8/2.3</td><td>2225</td><td>33%</td><td>95.39%</td><td>22X</td><td>Filtration</td><td>Petersen et al., 2011, this study</td></tr><tr><td>B. azoricus symbiont (BazSymB)</td><td>Menez Gwen</td><td>239</td><td>38.20</td><td>1.5/1.7</td><td>1802</td><td>30%</td><td>90.60%</td><td>8X</td><td>Gradient centrifugation/binning</td><td>This study</td></tr><tr><td>B. azoricus symbiont (BazSymA)*</td><td>Menez Gwen</td><td>506</td><td>37.58</td><td>1.85/1.85</td><td>2008</td><td>35%</td><td>97.70%</td><td>59X</td><td>Binning</td><td>This study</td></tr><tr><td>Ca. V. okutanii</td><td>Sagami Bay</td><td>1</td><td>31.59</td><td>1.0/1.0</td><td>980</td><td>26%</td><td>93.58%</td><td>-</td><td>Whole genome assembly</td><td>Kuwahara et al., 2007</td></tr><tr><td>Ca. R. magnifica</td><td>East Pacific Rise, 9°N</td><td>1</td><td>34.03</td><td>1.2/1.2</td><td>1210</td><td>23%</td><td>94.84%</td><td>-</td><td>Whole genome assembly</td><td>Newton et al., 2007</td></tr><tr><td>SUP05</td><td>Saanich Inlet</td><td>97</td><td>39.29</td><td>1.4/2.5</td><td>1586</td><td>30%</td><td>85.76%</td><td>-</td><td>Binning</td><td>Walsh et al., 2009</td></tr></tbody></table>

SOX, sulfur-oxidizing.

*SOX symbiont sequences recovered from metagenome of adductor muscle.

HGT = Genes that potentially originated from horizontal gene transfer.

†Length is the total length of sequence information on contigs without Ns, and span is the entire length of scaffold assembly including Ns.

‡The completeness of the genome was estimated with CheckM using a set of lineage-specific genes for proteobacteria (Parks et al., 2015).

§Median coverage.

\#Separation method indicates the experimental separation of symbionts from host tissue and co-occurring symbionts (filtration or gradient centrifugation), or the in silico separation of genomic information from hosts and co-occurring bacteria (binning).

DOI: 10.7554/eLife.07966.004
transposases or integrases in the clam symbiont genomes. The *Bathymodiolus* symbionts were highly enriched in restriction-modification system genes (between 10 and 22 genes), whereas SUP05 only had one, and the clam symbionts had none. This large difference raises the possibility that restriction-modification systems are involved in genome reshuffling in the *Bathymodiolus* symbionts.

## Gene-based comparison reveals toxin-related genes specific to *Bathymodiolus* symbionts

Between 2.3 and 7.6% of the genes found only in the *Bathymodiolus* SOX symbionts but not in the clam symbionts and SUP05, genomes were annotated as toxin or virulence genes (Figure 2). Most were related to genes from one of three toxin classes: (1) the RTX (repeats in toxins) toxins, (2) MARTX (multifunctional autoprocessing RTX toxins), a sub-group of RTX toxins, and (3) YD repeat toxins (also called *rhs* genes as they were initially described as 'recombination hotspots'). Representatives from all three toxin-related genes (TRGs) classes were found in each of the three-draft genomes, except for RTX, which were not found in BazSymA. The number of genes from each class varied between the three genomes. We found the largest number of TRGs in the genome with the fewest contigs, BspSym, which had at least 33 YD repeat genes, eight RTX genes, and 19 MARTX-like genes. In the BazSymB genome, 14 YD repeat genes, two RTX genes, and up to 10 MARTX genes were found. BazSymA had 16 YD repeat genes, and one MARTX (Figure 2, Supplementary file 1B). This indicates that these toxin-related classes are common to the SOX symbionts of both B. sp. and B. azoricus. In the BspSym genome, which assembled into the largest contigs, 22 out of 88 TRGs were found directly upstream or downstream of mobile elements.

Toxin genes are known to have unusually high substitution rates due to an 'evolutionary arms race' with their targets (Jackson et al., 2009; Linhartová et al., 2010). Accordingly, many of the *Bathymodiolus* symbiont TRGs were highly variable between the symbionts of *B. azoricus* and *B. sp.*, but also between the symbionts from the two *B. azoricus* individuals, and even between different copies within one genome (see 'Variability of TRGs within *Bathymodiolus* SOX symbiont populations'). We therefore searched for homologs of the TRGs in the clam symbionts and SUP05 with a lower BSR of up to 0.25, but even with this reduced stringency, no hits were found (see 'Materials and methods').

We ruled out the possibility that the TRGs were not found in SUP05 draft genome because of its incompleteness (~89% complete), by searching for homologs of the *Bathymodiolus* symbionts TRGs in unbinned metagenomes and metatranscriptomes from hydrothermal plumes and OMZs that are enriched in SUP05. If free-living SUP05 also encoded these TRGs, we would expect to find them regularly in SUP05-enriched metagenomes and metatranscriptomes. Instead, no hits were found in four out of these six data sets. In a metagenome from the Lost City hydrothermal vent, we found one weak hit to a YD repeat gene (31% similarity). In a metagenome from the Guaymas Basin hydrothermal vent, we found one weak hit to an RTX gene (34% similarity). However, both of these metagenomes were from sites colonized by either *Bathymodiolus* mussels (Lost City), or *Riftia pachyptila* tubeworms (Guaymas Basin), whose symbionts also encode a hemolysin gene of the RTX class (Gardebrecht et al., 2012). These rare hits might therefore come from contamination by symbionts in the environment (Harmer et al., 2008). Considering the almost complete absence of genes similar to the TRGs of *Bathymodiolus* in SUP05-enriched next-generation sequence data sets, we conclude that these genes are specific to the *Bathymodiolus* SOX symbionts and are not found in their close symbiotic or free-living relatives.

## Relationships to other toxins

Most symbiont TRGs were so divergent that they could not be confidently aligned. One exception was the YD repeat genes, a few of which contained a conserved repeat region. We reconstructed the phylogeny of this conserved region. YD repeat sequences from the symbionts of both *Bathymodiolus* species formed a distinct cluster, distant from all other sequences in public databases (Figure 3, Figure 3—figure supplement 1). The *Bathymodiolus* symbiont sequences did not cluster according to their host species, but instead were intermixed, suggesting gene duplication events prior to the divergence of the *B. azoricus* and *B. sp.* symbiont lineages. The *Bathymodiolus* symbiont sequences fell into a cluster that contained mostly pathogenic bacteria such as *Yersinia pestis* and *Burkholderia pseudomallei*. This cluster was well supported by Bayesian analysis (0.91 posterior probability). This cluster also contained a number of beneficial symbionts such as *B. rhizoxinica*, which is an intracellular
transposases or integrases in the clam symbiont genomes. The *Bathymodiolus* symbionts were highly enriched in restriction-modification system genes (between 10 and 22 genes), whereas SUP05 only had one, and the clam symbionts had none. This large difference raises the possibility that restriction-modification systems are involved in genome reshuffling in the *Bathymodiolus* symbionts.

## Gene-based comparison reveals toxin-related genes specific to *Bathymodiolus* symbionts

Between 2.3 and 7.6% of the genes found only in the *Bathymodiolus* SOX symbionts but not in the clam symbionts and SUP05, genomes were annotated as toxin or virulence genes (Figure 2). Most were related to genes from one of three toxin classes: (1) the RTX (repeats in toxins) toxins, (2) MARTX (multifunctional autoprocessing RTX toxins), a sub-group of RTX toxins, and (3) YD repeat toxins (also called *rhs* genes as they were initially described as 'recombination hotspots'). Representatives from all three toxin-related genes (TRGs) classes were found in each of the three-draft genomes, except for RTX, which were not found in BazSymA. The number of genes from each class varied between the three genomes. We found the largest number of TRGs in the genome with the fewest contigs, BspSym, which had at least 33 YD repeat genes, eight RTX genes, and 19 MARTX-like genes. In the BazSymB genome, 14 YD repeat genes, two RTX genes, and up to 10 MARTX genes were found. BazSymA had 16 YD repeat genes, and one MARTX (Figure 2, Supplementary file 1B). This indicates that these toxin-related classes are common to the SOX symbionts of both B. sp. and B. azoricus. In the BspSym genome, which assembled into the largest contigs, 22 out of 88 TRGs were found directly upstream or downstream of mobile elements.

Toxin genes are known to have unusually high substitution rates due to an 'evolutionary arms race' with their targets (Jackson et al., 2009; Linhartová et al., 2010). Accordingly, many of the *Bathymodiolus* symbiont TRGs were highly variable between the symbionts of *B. azoricus* and *B. sp.*, but also between the symbionts from the two *B. azoricus* individuals, and even between different copies within one genome (see 'Variability of TRGs within *Bathymodiolus* SOX symbiont populations'). We therefore searched for homologs of the TRGs in the clam symbionts and SUP05 with a lower BSR of up to 0.25, but even with this reduced stringency, no hits were found (see 'Materials and methods').

We ruled out the possibility that the TRGs were not found in SUP05 draft genome because of its incompleteness (~89% complete), by searching for homologs of the *Bathymodiolus* symbionts TRGs in unbinned metagenomes and metatranscriptomes from hydrothermal plumes and OMZs that are enriched in SUP05. If free-living SUP05 also encoded these TRGs, we would expect to find them regularly in SUP05-enriched metagenomes and metatranscriptomes. Instead, no hits were found in four out of these six data sets. In a metagenome from the Lost City hydrothermal vent, we found one weak hit to a YD repeat gene (31% similarity). In a metagenome from the Guaymas Basin hydrothermal vent, we found one weak hit to an RTX gene (34% similarity). However, both of these metagenomes were from sites colonized by either *Bathymodiolus* mussels (Lost City), or *Riftia pachyptila* tubeworms (Guaymas Basin), whose symbionts also encode a hemolysin gene of the RTX class (Gardebrecht et al., 2012). These rare hits might therefore come from contamination by symbionts in the environment (Harmer et al., 2008). Considering the almost complete absence of genes similar to the TRGs of *Bathymodiolus* in SUP05-enriched next-generation sequence data sets, we conclude that these genes are specific to the *Bathymodiolus* SOX symbionts and are not found in their close symbiotic or free-living relatives.

## Relationships to other toxins

Most symbiont TRGs were so divergent that they could not be confidently aligned. One exception was the YD repeat genes, a few of which contained a conserved repeat region. We reconstructed the phylogeny of this conserved region. YD repeat sequences from the symbionts of both *Bathymodiolus* species formed a distinct cluster, distant from all other sequences in public databases (Figure 3, Figure 3—figure supplement 1). The *Bathymodiolus* symbiont sequences did not cluster according to their host species, but instead were intermixed, suggesting gene duplication events prior to the divergence of the *B. azoricus* and *B. sp.* symbiont lineages. The *Bathymodiolus* symbiont sequences fell into a cluster that contained mostly pathogenic bacteria such as *Yersinia pestis* and *Burkholderia pseudomallei*. This cluster was well supported by Bayesian analysis (0.91 posterior probability). This cluster also contained a number of beneficial symbionts such as *B. rhizoxinica*, which is an intracellular
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Genes shared between the *Bathymodiolus* and vesicomyid SOX symbionts and free-living SUP05. Protein-coding sequences from the *Bathymodiolus* sulfur-oxidizing (SOX) symbiont were compared to the clam symbiont genomes and to the SUP05 metagenome from *Walsh et al.* (2009) with BLAST score ratios (BSR). (A) Venn diagram of the shared and unique gene content in the clam symbionts, mussel symbionts, and SUP05 bacteria. Predicted protein sequences of each mussel symbiont were compared to a combined data set of the clam symbionts (Rma and Vok) and SUP05. Similarly, protein sequences of each clam symbiont were compared to a combined data set of mussel symbionts (BspSym, BazSymB, and BazSymA). Depending on the reference genome, the number of shared genes varies slightly and possibly reflects the presence of paralogous genes and redundant sequence information in these draft genomes. Abbreviations are explained in detail in Table 1. The BLAST score ratio (BSR) threshold was 0.4. (B) Venn diagram of mussel symbiont toxin-related genes (TRGs), calculated with a BSR threshold of 0.2.

DOI: 10.7554/eLife.07966.005</figcaption></figure>

The following figure supplements are available for figure 2:

**Figure supplement 1.** Maximum likelihood 16S rRNA phylogeny of the close relatives of the *Bathymodiolus* SOX symbionts.

DOI: 10.7554/eLife.07966.006

**Figure supplement 2.** Whole genome alignment.

DOI: 10.7554/eLife.07966.007

Figure 2. continued on next page
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Genes shared between the *Bathymodiolus* and vesicomyid SOX symbionts and free-living SUP05. Protein-coding sequences from the *Bathymodiolus* sulfur-oxidizing (SOX) symbiont were compared to the clam symbiont genomes and to the SUP05 metagenome from Walsh et al. (2009) with BLAST score ratios (BSR). (A) Venn diagram of the shared and unique gene content in the clam symbionts, mussel symbionts, and SUP05 bacteria. Predicted protein sequences of each mussel symbiont were compared to a combined data set of the clam symbionts (Rma and Vok) and SUP05. Similarly, protein sequences of each clam symbiont were compared to a combined data set of mussel symbionts (BspSym, BazSymB, and BazSymA). Depending on the reference genome, the number of shared genes varies slightly and possibly reflects the presence of paralogous genes and redundant sequence information in these draft genomes. Abbreviations are explained in detail in Table 1. The BLAST score ratio (BSR) threshold was 0.4. (B) Venn diagram of mussel symbiont toxin-related genes (TRGs), calculated with a BSR threshold of 0.2.

DOI: 10.7554/eLife.07966.005</figcaption></figure>

The following figure supplements are available for figure 2:

**Figure supplement 1.** Maximum likelihood 16S rRNA phylogeny of the close relatives of the *Bathymodiolus* SOX symbionts.

DOI: 10.7554/eLife.07966.006

**Figure supplement 2.** Whole genome alignment.

DOI: 10.7554/eLife.07966.007

Figure 2. continued on next page
Figure 2. Continued

**Figure supplement 3.** Metabolic reconstruction of the Bathymodiolus symbiont.

DOI: 10.7554/eLife.07966.008

symbiont of the fungus *Rhizopus*, and the *Photorhabdus* and *Xenorhabdus* symbionts of soil nematodes (Waterfield et al., 2001; Goodrich-Blair and Clarke, 2007; Moebius et al., 2014).

To overcome the difficulties in aligning these highly divergent TRGs, we constructed gene sequence similarity networks based on BLAST to depict relationships among and between the symbiont TRGs, and those in public databases. This analysis revealed distinct sequence clusters that contained genes with >25% similarity over at least half of the length of the gene (Figure 4). If a cluster

<figure><img src="image_6.png" /><figcaption>**Figure 3.** Phylogeny of YD repeat-containing proteins. The tree is a consensus of bayesian and maximum likelihood analyses, result of an alignment of 536 amino acids. Black circles represent branches with posterior probability >0.8 and bootstrap value >80. White circles represent branches with either posterior probability >0.8 or bootstrap value >80. The number of sequences per collapsed group is shown next to the gray bloks. Purple: organism found in intestinal microflora or in close association with another organism; green: free-living; red: pathogen.</figcaption></figure>

DOI: 10.7554/eLife.07966.009

The following figure supplement is available for figure 3:

**Figure supplement 1.** Consensus of bayesian and maximum likelihood phylogeny of YD proteins with identifiers.

DOI: 10.7554/eLife.07966.010
Figure 2. Continued

**Figure supplement 3.** Metabolic reconstruction of the Bathymodiolus symbiont.

DOI: 10.7554/eLife.07966.008

symbiont of the fungus *Rhizopus*, and the *Photorhabdus* and *Xenorhabdus* symbionts of soil nematodes (Waterfield *et al.*, 2001; Goodrich-Blair and Clarke, 2007; Moebius *et al.*, 2014).

To overcome the difficulties in aligning these highly divergent TRGs, we constructed gene sequence similarity networks based on BLAST to depict relationships among and between the symbiont TRGs, and those in public databases. This analysis revealed distinct sequence clusters that contained genes with >25% similarity over at least half of the length of the gene (Figure 4). If a cluster

<figure><img src="image_6.png" /><figcaption>**Figure 3.** Phylogeny of YD repeat-containing proteins. The tree is a consensus of bayesian and maximum likelihood analyses, result of an alignment of 536 amino acids. Black circles represent branches with posterior probability >0.8 and bootstrap value >80. White circles represent branches with either posterior probability >0.8 or bootstrap value >80. The number of sequences per collapsed group is shown next to the gray bloks. Purple: organism found in intestinal microflora or in close association with another organism; green: free-living; red: pathogen.</figcaption></figure>

DOI: 10.7554/eLife.07966.009

The following figure supplement is available for figure 3:

**Figure supplement 1.** Consensus of bayesian and maximum likelihood phylogeny of YD proteins with identifiers.

DOI: 10.7554/eLife.07966.010
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Protein similarity network of toxin-related proteins in the *Bathymodiolus* symbionts. Each node corresponds to a protein sequence and the links between nodes represent BLAST hits. The length of the edges is inversely proportional to the sequence similarity. Protein clusters containing RTX or multifunctional autoprocessing RTX (MARTX) proteins are shown in the red panel on the left, and sequence clusters containing YD repeats are shown in the gray panel on the right. Arrowheads are proteins from *B. azoricus* symbionts, and triangles are proteins from *B. sp.* symbionts. The symbols are colored in green if they were identified in the *Bathymodiolus* symbionts as YD repeat-containing genes, red if they were identified as RTX genes, and purple for MARTX genes. Some protein sequences were similar to the TRGs but not annotated as such as these are partial genes that did not have any conserved domain. If the clusters contained mostly genes with a particular annotation, we named the clusters after these annotations, for example, cluster 'TcB/TcC' contained proteins annotated as TcB or TcC.
DOI: 10.7554/eLife.07966.011</figcaption></figure>

Figure 4. continued on next page
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Protein similarity network of toxin-related proteins in the *Bathymodiolus* symbionts. Each node corresponds to a protein sequence and the links between nodes represent BLAST hits. The length of the edges is inversely proportional to the sequence similarity. Protein clusters containing RTX or multifunctional autoprocessing RTX (MARTX) proteins are shown in the red panel on the left, and sequence clusters containing YD repeats are shown in the gray panel on the right. Arrowheads are proteins from *B. azoricus* symbionts, and triangles are proteins from *B. sp.* symbionts. The symbols are colored in green if they were identified in the *Bathymodiolus* symbionts as YD repeat-containing genes, red if they were identified as RTX genes, and purple for MARTX genes. Some protein sequences were similar to the TRGs but not annotated as such as these are partial genes that did not have any conserved domain. If the clusters contained mostly genes with a particular annotation, we named the clusters after these annotations, for example, cluster 'TcB/TcC' contained proteins annotated as TcB or TcC.
DOI: 10.7554/eLife.07966.011</figcaption></figure>

Figure 4. continued on next page
Figure 4. Continued

The following figure supplements are available for figure 4:

**Figure supplement 1.** Network of toxin-related proteins in the *Bathymodiolus* symbionts with BLAST hits from Vibrio, Photorhabdus, Xenorhabdus, and Pseudomonas highlighted.

DOI: 10.7554/eLife.07966.012

**Figure supplement 2.** Genomic architecture of MARTX regions.

DOI: 10.7554/eLife.07966.013

contained at least one gene that was similar to at least one other gene in another cluster (similarity cut-offs as above), then these clusters were joined to create a larger network. They were also joined if both clusters contained genes that had similarity to another gene in the database (i.e., if they could be joined by at most two steps). This allowed us to identify distinct sub-groups within the three toxin-related classes, and to identify toxin sequences from public databases that were most similar to the *Bathymodiolus* symbiont genes.

TRGs from the *Bathymodiolus* symbionts clustered together with toxin and TRGs from phylogenetically diverse organisms including characterized toxins of gammaproteobacterial Vibrio and Pseudomonas, and TRGs of the gammaproteobacteria Shewanella, the actinobacterial Rhodococcus, the cyanobacterium Trichodesmium, and the firmicute Caldicellulosiruptor (see e.g., **Figure 4—figure supplement 1**). The RTX genes clustered into two separate networks, one that had similarity to RTX secretion and activation genes, and one that had similarity to RTX toxins. Five distinct networks contained MARTX genes. One of these included genes from the symbiont MARTX1 cluster, and genes from other organisms that were annotated as MARTX or filamentous hemagglutinin. One network contained some genes that we classified as RTX and some we classified as MARTX, reflecting their shared features such as the RTX repeats. Eight MARTX genes had no significant hits to any other gene in public databases.

The YD repeat genes formed five distinct networks. Sequences from the first three had structural similarity to TcB and TcC, two subunits of the ABC toxins of Photorhabdus and Xenorhabdus, the beneficial symbionts of entomopathogenic nematodes. The B and C subunits form a cage-like structure that encapsulates the toxic domain (an adenosine diphosphate (ADP) ribosylation domain, located at the C terminus of the C subunit). The A subunit forms a syringe-like structure, which delivers the toxin to the insect cell (*Meusch et al., 2014*). The genes in the fourth YD network had structural similarity to TcA genes that encode the syringe-like A subunit. The fifth YD network had similarity to genes annotated as RhsB, which was shown to play a role in bacteria-bacteria competition in *Escherichia coli* (*Poole et al., 2011*).

## MARTX and YD repeat genes are enriched in the genomes of host-associated bacteria

The *Bathymodiolus* symbiont genomes encoded more YD repeat and MARTX genes than any other genome that we compared them to (Figure 5, Figure 5—figure supplements 1–3). This is remarkable considering the relatively small size of their genomes, and the fact that they are still incomplete. A few published genomes encoded more RTX genes, but these were much larger (>5 Mbp) (Figure 5—figure supplement 3).

The vast majority of RTX, MARTX, and YD repeat proteins have not been functionally characterized. The few proteins whose function has been studied in detail are from bacteria that are known pathogens or cultured strains that can form biofilms. Because of this, it is generally assumed that RTX, MARTX, and YD repeat proteins function in host-microbe interactions, in microbe-microbe antagonism, or in biofilm formation, but this has not been extensively tested. To further investigate the functional role of the TRGs encoded by the *Bathymodiolus* symbionts, we tested whether similar genes are more likely to be found in bacteria that live in a particular niche (extracellular host-associated, intracellular host-associated, or free-living), or that express a particular phenotype (pathogenesis or biofilm formation).

First, we used the Kruskal-Wallis one-way analysis of variance to determine whether the distribution of the three TRGs classes differed significantly (I) between biofilm-forming vs non-biofilm-forming bacteria, (II) between pathogenic and non-pathogenic bacteria, and (III) between free-
Figure 4. Continued

The following figure supplements are available for figure 4:

**Figure supplement 1.** Network of toxin-related proteins in the *Bathymodiolus* symbionts with BLAST hits from Vibrio, Photorhabdus, Xenorhabdus, and Pseudomonas highlighted.

DOI: 10.7554/eLife.07966.012

**Figure supplement 2.** Genomic architecture of MARTX regions.

DOI: 10.7554/eLife.07966.013

contained at least one gene that was similar to at least one other gene in another cluster (similarity cut-offs as above), then these clusters were joined to create a larger network. They were also joined if both clusters contained genes that had similarity to another gene in the database (i.e., if they could be joined by at most two steps). This allowed us to identify distinct sub-groups within the three toxin-related classes, and to identify toxin sequences from public databases that were most similar to the *Bathymodiolus* symbiont genes.

TRGs from the *Bathymodiolus* symbionts clustered together with toxin and TRGs from phylogenetically diverse organisms including characterized toxins of gammaproteobacterial Vibrio and Pseudomonas, and TRGs of the gammaproteobacteria Shewanella, the actinobacterial Rhodococcus, the cyanobacterium Trichodesmium, and the firmicute Caldicellulosiruptor (see e.g., **Figure 4—figure supplement 1**). The RTX genes clustered into two separate networks, one that had similarity to RTX secretion and activation genes, and one that had similarity to RTX toxins. Five distinct networks contained MARTX genes. One of these included genes from the symbiont MARTX1 cluster, and genes from other organisms that were annotated as MARTX or filamentous hemagglutinin. One network contained some genes that we classified as RTX and some we classified as MARTX, reflecting their shared features such as the RTX repeats. Eight MARTX genes had no significant hits to any other gene in public databases.

The YD repeat genes formed five distinct networks. Sequences from the first three had structural similarity to TcB and TcC, two subunits of the ABC toxins of Photorhabdus and Xenorhabdus, the beneficial symbionts of entomopathogenic nematodes. The B and C subunits form a cage-like structure that encapsulates the toxic domain (an adenosine diphosphate (ADP) ribosylation domain, located at the C terminus of the C subunit). The A subunit forms a syringe-like structure, which delivers the toxin to the insect cell (*Meusch et al., 2014*). The genes in the fourth YD network had structural similarity to TcA genes that encode the syringe-like A subunit. The fifth YD network had similarity to genes annotated as RhsB, which was shown to play a role in bacteria–bacteria competition in *Escherichia coli* (*Poole et al., 2011*).

## MARTX and YD repeat genes are enriched in the genomes of host-associated bacteria

The *Bathymodiolus* symbiont genomes encoded more YD repeat and MARTX genes than any other genome that we compared them to (Figure 5, Figure 5—figure supplements 1–3). This is remarkable considering the relatively small size of their genomes, and the fact that they are still incomplete. A few published genomes encoded more RTX genes, but these were much larger (>5 Mbp) (Figure 5—figure supplement 3).

The vast majority of RTX, MARTX, and YD repeat proteins have not been functionally characterized. The few proteins whose function has been studied in detail are from bacteria that are known pathogens or cultured strains that can form biofilms. Because of this, it is generally assumed that RTX, MARTX, and YD repeat proteins function in host-microbe interactions, in microbe–microbe antagonism, or in biofilm formation, but this has not been extensively tested. To further investigate the functional role of the TRGs encoded by the *Bathymodiolus* symbionts, we tested whether similar genes are more likely to be found in bacteria that live in a particular niche (extracellular host-associated, intracellular host-associated, or free-living), or that express a particular phenotype (pathogenesis or biofilm formation).

First, we used the Kruskal–Wallis one-way analysis of variance to determine whether the distribution of the three TRGs classes differed significantly (I) between biofilm-forming vs non-biofilm-forming bacteria, (II) between pathogenic and non-pathogenic bacteria, and (III) between free-
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Distribution of the three major TRGs classes according to lifestyle. Each dot represents one sequenced genome. The sum of TRGs is on the Y axis, and the total number of genes predicted in each genome is on the X axis. Free-living bacteria are shown in red, host-associated bacteria that live outside of host cells are in green, and host-associated bacteria that can live inside host cells are shown in blue. The positions of the Bathymodiolus SOX symbionts are indicated. A detailed overview of all organisms that had similar TRGs to the SOX symbiont with the number of TRGs is shown in *Supplementary file 1B*.</figcaption></figure>

DOI: 10.7554/eLife.07966.014

The following figure supplements are available for figure 5:

**Figure supplement 1.** YD genes per genome, normalized to the total gene count.

DOI: 10.7554/eLife.07966.015

Figure 5. continued on next page
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Distribution of the three major TRGs classes according to lifestyle. Each dot represents one sequenced genome. The sum of TRGs is on the Y axis, and the total number of genes predicted in each genome is on the X axis. Free-living bacteria are shown in red, host-associated bacteria that live outside of host cells are in green, and host-associated bacteria that can live inside host cells are shown in blue. The positions of the Bathymodiolus SOX symbionts are indicated. A detailed overview of all organisms that had similar TRGs to the SOX symbiont with the number of TRGs is shown in *Supplementary file 1B*.</figcaption></figure>

DOI: 10.7554/eLife.07966.014

The following figure supplements are available for figure 5:

**Figure supplement 1.** YD genes per genome, normalized to the total gene count.

DOI: 10.7554/eLife.07966.015

Figure 5. continued on next page
Figure 5. Continued

**Figure supplement 2.** MARTX genes per genome, normalized to the total gene count.

DOI: 10.7554/eLife.07966.016

**Figure supplement 3.** RTX genes per genome, normalized to the total gene count.

DOI: 10.7554/eLife.07966.017

living bacteria, host-associated intracellular bacteria, and host-associated extracellular bacteria, without considering whether the bacteria were pathogenic. There was no significant enrichment of any TRG category in bacteria known to form biofilms vs those that do not (Table 2). One class, MARTX was significantly enriched in the genomes of pathogenic vs non-pathogenic bacteria (p-value = 0.007, Kruskal-Wallis test). There was also a significant bias in the distribution of genes encoding MARTX and YD repeat genes in bacteria according to their lifestyle (extracellular host-associated, intracellular host-associated, or free-living).

When three categories are tested, such as in (III) above, the Kruskal-Wallis test does not identify which category the bias is associated with. To tease apart which of these three niche categories was most enriched in TRGs, we did Mann-Whitney-Wilcoxon tests (Table 3). These showed that both YD repeat and MARTX genes were enriched in the genomes of host-associated microbes (YD repeat: p-value = 0.026, MARTX: p-values = 2.125e-6, 1.618e-6, Mann-Whitney-Wilcoxon test). While MARTX genes were enriched in host-associated bacteria regardless of their location, YD repeat genes were only significantly enriched in intracellular bacteria. In contrast to YD repeat and MARTX genes, RTX did not show any enrichment in the three defined categories. RTX are therefore widely distributed among bacteria and are just as likely to be found in free-living and host-associated bacteria (Appendix 2).

Bacteria that are closely related often have similar genomic and physiological features. However, toxin genes are commonly gained through HGT, which may weaken the phylogenetic signal in their distribution patterns (reviewed by Dobrindt et al., 2004; Gogarten and Townsend, 2005). To tease apart the possible phylogenetic influence on the TRGs distribution, we used Permanova to test whether any of the three classes was enriched in particular phylogenetic groups at the class, order, and family levels. Only RTX genes were significantly enriched, and only at the order level (p-value = 0.0159) (Supplementary file 1C). Therefore, phylogeny is not the main driver in the toxin-related distribution of YD repeats genes and MARTX.

## Variability of TRGs within *Bathymodiolus* SOX symbiont populations

Toxin genes often have unusually high substitution rates, making them highly variable compared to non-toxin genes (e.g., Ohno et al., 1997; Davies et al., 2002). We compared the substitution rates of all genes in the *Bathymodiolus* SOX symbiont genomes within the population of symbionts associated with each mussel species. This was done for each of the two species, *B. azoricus* and *B. sp.* by mapping transcriptome reads from three *Bathymodiolus* individuals to the draft genomes of their respective symbionts (see below for transcriptomes). We calculated the number of single nucleotide polymorphisms (SNPs) per kb per gene. The number of SNPs in most of the TRGs was not significantly higher than the genome-wide average (Figure 6). However, we found 22 TRGs that did have

<table><caption>**Table 2.** p-values obtained with Kruskal-Wallis rank sum test</caption><thead><tr><th></th><th>B/NB df = 1</th><th>P/NP df = 1</th><th>Ext/Int/FL df = 2</th></tr></thead><tbody><tr><td>YD</td><td>0.097</td><td>0.5217</td><td>0.010*</td></tr><tr><td>RTX</td><td>0.715</td><td>0.793</td><td>0.308</td></tr><tr><td>MARTX</td><td>0.773</td><td>0.007*</td><td>3.21e-06*</td></tr></tbody></table>

The three main lifestyle categories were tested against each toxin-related class. Number of TRGs per genome was normalized to the total gene count.

FL = free-living, Ext = extracellular host-associated, Int = intracellular host-associated, P = pathogen, NP = non-pathogen, B = found in biofilms, NB = not found in biofilms, df = degrees of freedom, TRG, toxin-related gene, MARTX, multifunctional autoprocessing RTX.

*p-value was considered to be significant (p < 0.05).

DOI: 10.7554/eLife.07966.018
Figure 5. Continued

**Figure supplement 2.** MARTX genes per genome, normalized to the total gene count.

DOI: 10.7554/eLife.07966.016

**Figure supplement 3.** RTX genes per genome, normalized to the total gene count.

DOI: 10.7554/eLife.07966.017

living bacteria, host-associated intracellular bacteria, and host-associated extracellular bacteria, without considering whether the bacteria were pathogenic. There was no significant enrichment of any TRG category in bacteria known to form biofilms vs those that do not (Table 2). One class, MARTX was significantly enriched in the genomes of pathogenic vs non-pathogenic bacteria (p-value = 0.007, Kruskal-Wallis test). There was also a significant bias in the distribution of genes encoding MARTX and YD repeat genes in bacteria according to their lifestyle (extracellular host-associated, intracellular host-associated, or free-living).

When three categories are tested, such as in (III) above, the Kruskal-Wallis test does not identify which category the bias is associated with. To tease apart which of these three niche categories was most enriched in TRGs, we did Mann-Whitney-Wilcoxon tests (Table 3). These showed that both YD repeat and MARTX genes were enriched in the genomes of host-associated microbes (YD repeat: p-value = 0.026, MARTX: p-values = 2.125e-6, 1.618e-6, Mann-Whitney-Wilcoxon test). While MARTX genes were enriched in host-associated bacteria regardless of their location, YD repeat genes were only significantly enriched in intracellular bacteria. In contrast to YD repeat and MARTX genes, RTX did not show any enrichment in the three defined categories. RTX are therefore widely distributed among bacteria and are just as likely to be found in free-living and host-associated bacteria (Appendix 2).

Bacteria that are closely related often have similar genomic and physiological features. However, toxin genes are commonly gained through HGT, which may weaken the phylogenetic signal in their distribution patterns (reviewed by Dobrindt et al., 2004; Gogarten and Townsend, 2005). To tease apart the possible phylogenetic influence on the TRGs distribution, we used Permanova to test whether any of the three classes was enriched in particular phylogenetic groups at the class, order, and family levels. Only RTX genes were significantly enriched, and only at the order level (p-value = 0.0159) (Supplementary file 1C). Therefore, phylogeny is not the main driver in the toxin-related distribution of YD repeats genes and MARTX.

## Variability of TRGs within *Bathymodiolus* SOX symbiont populations

Toxin genes often have unusually high substitution rates, making them highly variable compared to non-toxin genes (e.g., Ohno et al., 1997; Davies et al., 2002). We compared the substitution rates of all genes in the *Bathymodiolus* SOX symbiont genomes within the population of symbionts associated with each mussel species. This was done for each of the two species, *B. azoricus* and *B. sp.* by mapping transcriptome reads from three *Bathymodiolus* individuals to the draft genomes of their respective symbionts (see below for transcriptomes). We calculated the number of single nucleotide polymorphisms (SNPs) per kb per gene. The number of SNPs in most of the TRGs was not significantly higher than the genome-wide average (Figure 6). However, we found 22 TRGs that did have

<table><caption>**Table 2.** p-values obtained with Kruskal-Wallis rank sum test</caption><thead><tr><th></th><th>B/NB df = 1</th><th>P/NP df = 1</th><th>Ext/Int/FL df = 2</th></tr></thead><tbody><tr><td>YD</td><td>0.097</td><td>0.5217</td><td>0.010*</td></tr><tr><td>RTX</td><td>0.715</td><td>0.793</td><td>0.308</td></tr><tr><td>MARTX</td><td>0.773</td><td>0.007*</td><td>3.21e-06*</td></tr></tbody></table>

The three main lifestyle categories were tested against each toxin-related class. Number of TRGs per genome was normalized to the total gene count.

FL = free-living, Ext = extracellular host-associated, Int = intracellular host-associated, P = pathogen, NP = non-pathogen, B = found in biofilms, NB = not found in biofilms, df = degrees of freedom, TRG, toxin-related gene, MARTX, multifunctional autoprocessing RTX.

*p-value was considered to be significant (p < 0.05).

DOI: 10.7554/eLife.07966.018
<table><caption>**Table 3.** p-values obtained with Mann-Whitney-Wilcoxon test for enrichment of YD and MARTX genes similar to those from the SOX symbiont</caption><thead><tr><th></th><th>FL/Ext</th><th>FL/Int</th><th>Ext/Int</th></tr></thead><tbody><tr><td>YD</td><td>0.129</td><td>0.026</td><td>0.006*</td></tr><tr><td>MARTX</td><td>2.125e<sup>-06*</sup></td><td>1.618e<sup>-06*</sup></td><td>0.751</td></tr></tbody></table>

FL = free-living, Ext = extracellular host-associated, Int = intracellular host-associated, MARTX, multifunctional autoprocessing RTX, SOX, sulfur-oxidizing.

*p-value was considered to be significant (p < 0.05).

DOI: 10.7554/eLife.07966.019

significantly more SNPs than most of the other genes in the genomes. Among the 22 highly variable genes, we found representatives of each TRG class: YD repeats, RTX, and MARTX (Figure 6).

## Expression of TRGs

Transcriptome sequencing revealed that all predicted TRGs of the SOX symbionts in *B. azoricus* and *B. sp. gills* were expressed. Reads mapping to TRGs accounted for 0.67–1.71% of mRNA in *B. azoricus* symbionts and 0.58–3.14% in *B. sp. symbionts*. All TRGs were found in the transcriptomes of at least one of the three

individuals that we sequenced per species (*Supplementary file 1D*). The expression levels of some genes from the RTX, MARTX, and YD repeats classes were in some cases higher than the expression of the essential Calvin cycle gene ribulose bisphosphate carboxylase/oxidase (RuBisCO), which

<figure><img src="image_9.png" /><figcaption>**Figure 6.** Single nucleotide polymorphisms per gene. The number of single nucleotide polymorphisms (SNPs) per gene was normalized according to the length minus regions of unknown sequence for genes containing N's. Genes smaller than 150 bp were not considered. The dotted line represents the median plus one standard deviation of the number of SNPs per gene per kb.</figcaption></figure>

DOI: 10.7554/eLife.07966.020

The following source data is available for figure 6:

**Source data 1.** Variability in TRGs encoded by the *Bathymodiolus* SOX symbionts.

DOI: 10.7554/eLife.07966.021
<table><caption>**Table 3.** p-values obtained with Mann-Whitney-Wilcoxon test for enrichment of YD and MARTX genes similar to those from the SOX symbiont</caption><thead><tr><th></th><th>FL/Ext</th><th>FL/Int</th><th>Ext/Int</th></tr></thead><tbody><tr><td>YD</td><td>0.129</td><td>0.026</td><td>0.006*</td></tr><tr><td>MARTX</td><td>2.125e<sup>-06*</sup></td><td>1.618e<sup>-06*</sup></td><td>0.751</td></tr></tbody></table>

FL = free-living, Ext = extracellular host-associated, Int = intracellular host-associated, MARTX, multifunctional autoprocessing RTX, SOX, sulfur-oxidizing.

*p-value was considered to be significant (p < 0.05).

DOI: 10.7554/eLife.07966.019

significantly more SNPs than most of the other genes in the genomes. Among the 22 highly variable genes, we found representatives of each TRG class: YD repeats, RTX, and MARTX (Figure 6).

## Expression of TRGs

Transcriptome sequencing revealed that all predicted TRGs of the SOX symbionts in *B. azoricus* and *B. sp. gills* were expressed. Reads mapping to TRGs accounted for 0.67–1.71% of mRNA in *B. azoricus* symbionts and 0.58–3.14% in *B. sp. symbionts*. All TRGs were found in the transcriptomes of at least one of the three

individuals that we sequenced per species (*Supplementary file 1D*). The expression levels of some genes from the RTX, MARTX, and YD repeats classes were in some cases higher than the expression of the essential Calvin cycle gene ribulose bisphosphate carboxylase/oxidase (RuBisCO), which

<figure><img src="image_9.png" /><figcaption>**Figure 6.** Single nucleotide polymorphisms per gene. The number of single nucleotide polymorphisms (SNPs) per gene was normalized according to the length minus regions of unknown sequence for genes containing N's. Genes smaller than 150 bp were not considered. The dotted line represents the median plus one standard deviation of the number of SNPs per gene per kb.</figcaption></figure>

DOI: 10.7554/eLife.07966.020

The following source data is available for figure 6:

**Source data 1.** Variability in TRGs encoded by the *Bathymodiolus* SOX symbionts.

DOI: 10.7554/eLife.07966.021
accounted for 0.03–0.5% of mRNA in *B. azoricus* symbionts and 0.11–1.04% in *B. sp.* symbionts. We could also confirm expression of some TRGs at the protein level in *B. azoricus* mussels (samples of *B. sp.* were not available for proteomics). We analyzed (I) gradient centrifugation fractions that were enriched in the *Bathymodiolus* SOX symbiont and (II) whole gill tissue. Some of the characterized toxin proteins are associated with membranes. We therefore also analyzed a fraction enriched in membrane proteins. In these proteomes, we found 12 YD repeat proteins and one RTX. We also found a number of toxin-associated proteins such as one RTX activator and two RTX toxin transporters (Table 4). Nine of the 12 identified YD repeat proteins were present in the membrane proteome. Our method allowed us to identify symbiont-encoded proteins that were found in higher relative abundance in the whole gill fraction compared to the symbiont-enriched fraction (see Appendix 4). These proteins are potentially exported to the host tissue, indicating that they play a role in direct host-symbiont interactions. One RTX activator and one YD repeat protein were significantly enriched in the host tissue (Table 4).

Discussion

Origin of TRGs in *Bathymodiolus* symbionts

Two scenarios could explain the origin of the large complement of TRGs in the *Bathymodiolus* symbionts and not in their close relatives: firstly, the TRGs could have been in the genome of their last

<table><caption>**Table 4.** Toxin-related proteins found in the proteome of the SOX symbiont from *B. azoricus*</caption><thead><tr><th>Proteome</th><th>Identifier</th><th>Category</th><th>Annotation</th><th>Molecular weight (kDa)</th><th>Max. number of unique peptides*</th></tr></thead><tbody><tr><td>SM</td><td>Host_EST_000107</td><td>YD</td><td>IPR006530IYD repeat</td><td>43</td><td>11</td></tr><tr><td>SM</td><td>Host_EST_000115</td><td>YD</td><td>IPR006530IYD repeat</td><td>42</td><td>12</td></tr><tr><td>N</td><td>Host_EST_000248</td><td>YD</td><td>IPR006530IYD repeat</td><td>37</td><td>7</td></tr><tr><td>M</td><td>Host_EST_002123</td><td>YD</td><td>IPR006530IYD repeat</td><td>24</td><td>5</td></tr><tr><td>SM†</td><td>Thio_BAZ_1943_contig360420_0</td><td>RTX (activator)</td><td>Hemolysin-activating lysine-acyltransferase (Hemolysin C)</td><td>19</td><td>3</td></tr><tr><td>SM</td><td>Tox_BAZ_119_contig00027_0</td><td>YD</td><td>RHS repeat-associated core domain-containing protein</td><td>202</td><td>17</td></tr><tr><td>SM</td><td>Tox_BAZ_120_contig00027_1</td><td>YD</td><td>Virulence plasmid 28.1 kDa A protein</td><td>62</td><td>11</td></tr><tr><td>SM</td><td>Tox_BAZ_1734_contig02141_2</td><td>RTX (transporter)</td><td>Secretion protein HlyD family protein</td><td>43</td><td>10</td></tr><tr><td>SM</td><td>Tox_BAZ_2494_contig00030_0</td><td>YD</td><td>Virulence plasmid 28.1 kDa A protein</td><td>183</td><td>33</td></tr><tr><td>M</td><td>Tox_BAZ_3202_scaffold00038_7</td><td>RTX</td><td>Hemolysins and related proteins containing CBS domains</td><td>35</td><td>2</td></tr><tr><td>SM</td><td>Tox_BAZ_525_contig104979_0</td><td>YD</td><td>Virulence plasmid 28.1 kDa A protein</td><td>52</td><td>2</td></tr><tr><td>S†</td><td>ToxAzor_892893</td><td>YD</td><td>Rhs</td><td>114</td><td>2</td></tr><tr><td>SM</td><td>ToxSMAR_1260BAT01109</td><td>YD</td><td>[weak similarity to] Toxin complex/plasmid virulence protein</td><td>321</td><td>8</td></tr><tr><td>N</td><td>ToxSMAR_2052BAT01788, Thio_BAZ_1733_contig02141_1 or Thio_BAZ_2580_scaffold00010_8</td><td>RTX (transporter)</td><td>Toxin secretion ATP-binding protein</td><td>79</td><td>5</td></tr><tr><td>S</td><td>ToxSMAR893-894</td><td>YD</td><td>Rhs family protein</td><td>103</td><td>1</td></tr><tr><td>SM</td><td>ToxAzor_890891</td><td>YD</td><td>Rhs family protein</td><td>67</td><td>6</td></tr></tbody></table>

S = soluble proteome, M = membrane-enriched proteome, SM = found in both proteomes, SOX, sulfur-oxidizing.

*The highest number of unique peptides detected in one sample.

†Proteins that are potentially exported by the symbiont to the host gill tissue.

DOI: 10.7554/eLife.07966.022
accounted for 0.03–0.5% of mRNA in *B. azoricus* symbionts and 0.11–1.04% in *B. sp.* symbionts. We could also confirm expression of some TRGs at the protein level in *B. azoricus* mussels (samples of *B. sp.* were not available for proteomics). We analyzed (I) gradient centrifugation fractions that were enriched in the *Bathymodiolus* SOX symbiont and (II) whole gill tissue. Some of the characterized toxin proteins are associated with membranes. We therefore also analyzed a fraction enriched in membrane proteins. In these proteomes, we found 12 YD repeat proteins and one RTX. We also found a number of toxin-associated proteins such as one RTX activator and two RTX toxin transporters (Table 4). Nine of the 12 identified YD repeat proteins were present in the membrane proteome. Our method allowed us to identify symbiont-encoded proteins that were found in higher relative abundance in the whole gill fraction compared to the symbiont-enriched fraction (see Appendix 4). These proteins are potentially exported to the host tissue, indicating that they play a role in direct host-symbiont interactions. One RTX activator and one YD repeat protein were significantly enriched in the host tissue (Table 4).

Discussion

Origin of TRGs in *Bathymodiolus* symbionts

Two scenarios could explain the origin of the large complement of TRGs in the *Bathymodiolus* symbionts and not in their close relatives: firstly, the TRGs could have been in the genome of their last

<table><caption>**Table 4.** Toxin-related proteins found in the proteome of the SOX symbiont from *B. azoricus*</caption><thead><tr><th>Proteome</th><th>Identifier</th><th>Category</th><th>Annotation</th><th>Molecular weight (kDa)</th><th>Max. number of unique peptides*</th></tr></thead><tbody><tr><td>SM</td><td>Host_EST_000107</td><td>YD</td><td>IPR006530IYD repeat</td><td>43</td><td>11</td></tr><tr><td>SM</td><td>Host_EST_000115</td><td>YD</td><td>IPR006530IYD repeat</td><td>42</td><td>12</td></tr><tr><td>N</td><td>Host_EST_000248</td><td>YD</td><td>IPR006530IYD repeat</td><td>37</td><td>7</td></tr><tr><td>M</td><td>Host_EST_002123</td><td>YD</td><td>IPR006530IYD repeat</td><td>24</td><td>5</td></tr><tr><td>SM†</td><td>Thio_BAZ_1943_contig360420_0</td><td>RTX (activator)</td><td>Hemolysin-activating lysine-acyltransferase (Hemolysin C)</td><td>19</td><td>3</td></tr><tr><td>SM</td><td>Tox_BAZ_119_contig00027_0</td><td>YD</td><td>RHS repeat-associated core domain-containing protein</td><td>202</td><td>17</td></tr><tr><td>SM</td><td>Tox_BAZ_120_contig00027_1</td><td>YD</td><td>Virulence plasmid 28.1 kDa A protein</td><td>62</td><td>11</td></tr><tr><td>SM</td><td>Tox_BAZ_1734_contig02141_2</td><td>RTX (transporter)</td><td>Secretion protein HlyD family protein</td><td>43</td><td>10</td></tr><tr><td>SM</td><td>Tox_BAZ_2494_contig00030_0</td><td>YD</td><td>Virulence plasmid 28.1 kDa A protein</td><td>183</td><td>33</td></tr><tr><td>M</td><td>Tox_BAZ_3202_scaffold00038_7</td><td>RTX</td><td>Hemolysins and related proteins containing CBS domains</td><td>35</td><td>2</td></tr><tr><td>SM</td><td>Tox_BAZ_525_contig104979_0</td><td>YD</td><td>Virulence plasmid 28.1 kDa A protein</td><td>52</td><td>2</td></tr><tr><td>S†</td><td>ToxAzor_892893</td><td>YD</td><td>Rhs</td><td>114</td><td>2</td></tr><tr><td>SM</td><td>ToxSMAR_1260BAT01109</td><td>YD</td><td>[weak similarity to] Toxin complex/plasmid virulence protein</td><td>321</td><td>8</td></tr><tr><td>N</td><td>ToxSMAR_2052BAT01788, Thio_BAZ_1733_contig02141_1 or Thio_BAZ_2580_scaffold00010_8</td><td>RTX (transporter)</td><td>Toxin secretion ATP-binding protein</td><td>79</td><td>5</td></tr><tr><td>S</td><td>ToxSMAR893-894</td><td>YD</td><td>Rhs family protein</td><td>103</td><td>1</td></tr><tr><td>SM</td><td>ToxAzor_890891</td><td>YD</td><td>Rhs family protein</td><td>67</td><td>6</td></tr></tbody></table>

S = soluble proteome, M = membrane-enriched proteome, SM = found in both proteomes, SOX, sulfur-oxidizing.

*The highest number of unique peptides detected in one sample.

†Proteins that are potentially exported by the symbiont to the host gill tissue.

DOI: 10.7554/eLife.07966.022
common ancestor but were all subsequently lost in both the clam symbionts and SUP05. Alternatively, the *Bathymodiolus* symbionts could have acquired these genes via HGT after their divergence from the clam symbionts and SUP05. Toxins are often found in 'genomic islands' that have been acquired through HGT (reviewed by Lindsay et al., 1998; Ochman et al., 2000; Dobrindt et al., 2004; Soucy et al., 2015). Several observations point towards HGT of the TRGs into the *Bathymodiolus* symbionts rather than their loss by the clam symbionts and SUP05. Firstly, 63-68% of the TRGs could be identified as potentially horizontally acquired based on codon usage analysis, in contrast to 30-35% predicted for all coding sequences (Table 1, Figure 6—source data 1). This means that their transfer was relatively recent, as the codon usage of these genes has not yet adapted to one typical of the symbiont genomes. Secondly, the content of the *Bathymodiolus* symbiont genomes attest to the major role that HGT has played in their evolution. They are enriched in mobile elements such as transposases and restriction-modification systems compared to their closest relatives (Supplementary file 1A). The lack of synteny we observed in the symbiont genomes is consistent with the presence of mobile elements and major HGT events (Kobayashi, 2001; Rocha and Danchin, 2002; Achaz et al., 2003). Thirdly, the TRGs from the *Bathymodiolus* symbionts are similar to genes from distantly related bacteria. Finally, mobile elements were regularly found directly upstream or downstream of the TRGs in the *Bathymodiolus* symbiont genomes. The linkage of mobile elements with some of these genes could explain the mechanism of their transfer into the *Bathymodiolus* symbiont genomes. This could also explain why each genome contained multiple copies of TRGs, as mobile elements are also prone to duplication (Reams and Neidle, 2004). Considering these observations, and the absence of these TRGs in their close relatives, we consider it most likely that they were acquired through HGT.

**Are some of the TRGs used for direct beneficial interactions between hosts and symbionts?**

Genes from the second major class of TRGs, MARTX, that were similar to those from *Bathymodiolus* symbiont, were significantly enriched in beneficial and pathogenic host-associated bacteria. One of the regions encoding MARTX genes from the *Bathymodiolus* SOX symbiont has a domain structure similar to the filamentous hemagglutinin FhaB from *B. pertussis* and *Bordetella bronchiseptica* (MARTX1, see Appendix 3). FhaB is involved in attachment of *Bordetella* to their human host and suppression of the immune response (reviewed by Melvin et al., 2014). *B. bronchiseptica* has two distinct phenotypic stages: an infective stage, where *fhaB* is upregulated, and a non-infective stage, where *fhaB* is downregulated. The non-infective stage is necessary for its survival in the environment outside of the host. The lifestyle of the *Bathymodiolus* symbiont has striking similarities to *B. bronchiseptica*. The symbionts must also survive in the environment to be transmitted from one host generation to the next. Our transcriptomes showed that these MARTX genes are expressed by *Bathymodiolus* symbionts within the host tissue. Unfortunately, we do not have samples to test whether they are downregulated in environmental symbiont stages. If so, it may have a similar function to its homologs in pathogenic *Bordetella*. MARTX-like genes also mediate cell-cell attachment in the symbiotic bacterial consortium '*Chlorochromatium aggregatum*' (Vogl et al., 2008; Liu et al., 2013).

**Are some of the YD repeat proteins used for antagonistic bacteria-bacteria interactions?**

Like the MARTX, YD repeat proteins were also significantly more enriched in host-associated compared to free-living bacteria. Some characterized YD repeat proteins function in competition between closely related bacterial strains (Waterfield et al., 2001; Kung et al., 2012; Koskiniemi et al., 2013). In its intracellular niche, why would the *Bathymodiolus* symbiont need to compete against other bacteria? Multiple strains of SOX symbionts can co-occur in *Bathymodiolus* mussels (Won et al., 2003a; DeChaine et al., 2006; Duperron et al., 2008). These may compete with each other for nutrients and energy, or for space within host cells. Bacteria that express toxins to inhibit their close relatives need an immunity protein to protect them from each toxin (Zhang et al., 2012; Benz and Meinhart, 2014). These immunity proteins are encoded immediately downstream or upstream of the toxin protein. The toxin-immunity pair is usually linked to genes encoding a type VI secretion system, which is the mechanism of toxin delivery for all so far described YD proteins involved in bacteria-bacteria competition (Zhang et al., 2012; Benz and Meinhart, 2014). None of the YD repeat proteins in the *Bathymodiolus* SOX symbiont genome was found in an operon with an
common ancestor but were all subsequently lost in both the clam symbionts and SUP05. Alternatively, the *Bathymodiolus* symbionts could have acquired these genes via HGT after their divergence from the clam symbionts and SUP05. Toxins are often found in 'genomic islands' that have been acquired through HGT (reviewed by Lindsay et al., 1998; Ochman et al., 2000; Dobrindt et al., 2004; Soucy et al., 2015). Several observations point towards HGT of the TRGs into the *Bathymodiolus* symbionts rather than their loss by the clam symbionts and SUP05. Firstly, 63-68% of the TRGs could be identified as potentially horizontally acquired based on codon usage analysis, in contrast to 30-35% predicted for all coding sequences (Table 1, Figure 6—source data 1). This means that their transfer was relatively recent, as the codon usage of these genes has not yet adapted to one typical of the symbiont genomes. Secondly, the content of the *Bathymodiolus* symbiont genomes attest to the major role that HGT has played in their evolution. They are enriched in mobile elements such as transposases and restriction-modification systems compared to their closest relatives (Supplementary file 1A). The lack of synteny we observed in the symbiont genomes is consistent with the presence of mobile elements and major HGT events (Kobayashi, 2001; Rocha and Danchin, 2002; Achaz et al., 2003). Thirdly, the TRGs from the *Bathymodiolus* symbionts are similar to genes from distantly related bacteria. Finally, mobile elements were regularly found directly upstream or downstream of the TRGs in the *Bathymodiolus* symbiont genomes. The linkage of mobile elements with some of these genes could explain the mechanism of their transfer into the *Bathymodiolus* symbiont genomes. This could also explain why each genome contained multiple copies of TRGs, as mobile elements are also prone to duplication (Reams and Neidle, 2004). Considering these observations, and the absence of these TRGs in their close relatives, we consider it most likely that they were acquired through HGT.

**Are some of the TRGs used for direct beneficial interactions between hosts and symbionts?**

Genes from the second major class of TRGs, MARTX, that were similar to those from *Bathymodiolus* symbiont, were significantly enriched in beneficial and pathogenic host-associated bacteria. One of the regions encoding MARTX genes from the *Bathymodiolus* SOX symbiont has a domain structure similar to the filamentous hemagglutinin FhaB from *B. pertussis* and *Bordetella bronchiseptica* (MARTX1, see Appendix 3). FhaB is involved in attachment of *Bordetella* to their human host and suppression of the immune response (reviewed by Melvin et al., 2014). *B. bronchiseptica* has two distinct phenotypic stages: an infective stage, where *fhaB* is upregulated, and a non-infective stage, where *fhaB* is downregulated. The non-infective stage is necessary for its survival in the environment outside of the host. The lifestyle of the *Bathymodiolus* symbiont has striking similarities to *B. bronchiseptica*. The symbionts must also survive in the environment to be transmitted from one host generation to the next. Our transcriptomes showed that these MARTX genes are expressed by *Bathymodiolus* symbionts within the host tissue. Unfortunately, we do not have samples to test whether they are downregulated in environmental symbiont stages. If so, it may have a similar function to its homologs in pathogenic *Bordetella*. MARTX-like genes also mediate cell-cell attachment in the symbiotic bacterial consortium '*Chlorochromatium aggregatum*' (Vogl et al., 2008; Liu et al., 2013).

**Are some of the YD repeat proteins used for antagonistic bacteria-bacteria interactions?**

Like the MARTX, YD repeat proteins were also significantly more enriched in host-associated compared to free-living bacteria. Some characterized YD repeat proteins function in competition between closely related bacterial strains (Waterfield et al., 2001; Kung et al., 2012; Koskiniemi et al., 2013). In its intracellular niche, why would the *Bathymodiolus* symbiont need to compete against other bacteria? Multiple strains of SOX symbionts can co-occur in *Bathymodiolus* mussels (Won et al., 2003a; DeChaine et al., 2006; Duperron et al., 2008). These may compete with each other for nutrients and energy, or for space within host cells. Bacteria that express toxins to inhibit their close relatives need an immunity protein to protect them from each toxin (Zhang et al., 2012; Benz and Meinhart, 2014). These immunity proteins are encoded immediately downstream or upstream of the toxin protein. The toxin-immunity pair is usually linked to genes encoding a type VI secretion system, which is the mechanism of toxin delivery for all so far described YD proteins involved in bacteria-bacteria competition (Zhang et al., 2012; Benz and Meinhart, 2014). None of the YD repeat proteins in the *Bathymodiolus* SOX symbiont genome was found in an operon with an
identifiable immunity protein, and no type VI secretion system gene was found in any of our draft genomes. The arrangement of YD repeat genes in the *Bathymodiolus* symbiont genomes and the lack of genes encoding type VI secretion systems are therefore inconsistent with a role in competition between closely related symbiont strains.

*Bathymodiolus* mussels are also infected by bacterial intranuclear pathogens called *Candidatus Endonucleobacter bathymodioli*, which are related to the genus *Endozoicomonas* (*Zielinski et al., 2009*). Ca. E. bathymodioli invades the mussel cell nuclei where it multiplies, eventually bursting the infected cell. Intranuclear bacteria are never found in the nuclei of symbiont-containing cells, which led *Zielinski et al.* (2009) to hypothesize that the symbionts can protect their host cells against infection. Consistent with this hypothesis, growth inhibition assays showed that *B. azoricus* gill tissue homogenates inhibit the growth of a broader spectrum of pathogens compared to the symbiont-free foot tissues (*Bettencourt et al., 2007*). The mechanisms of protection against cultured bacterial pathogens and Ca. E. bathymodioli are unknown. There is no evidence from their genomes that the SOX symbionts of *B. sp.* and *B. azoricus* produce antibiotics, but they do have genes for bacteriocin production that were expressed in the six transcriptomes analyzed in this study. Expression of some of the TRGs discovered here, for example, those related to toxins involved in bacteria-bacteria competition, or the production of bacteriocins by symbionts could explain the absence of intranuclear bacteria from symbiont-hosting cells.

**Are some of the TRGs used for protection against eukaryotic parasites?**

Symbionts of other marine invertebrate hosts are able to recognize, enter, and survive within host cells without a large number of TRGs. For example, the genomes of the SOX symbionts of hydrothermal vent *Riftia* tubeworms and the heterotrophic symbionts of whale-fall *Osedax* worms are virtually complete, but they only contain one or a few RTX genes, and no YD repeat genes (*Robidart et al., 2008; Gardebrecht et al., 2012; Goffredi et al., 2014*). There is overwhelming evidence that SOX symbionts are beneficial for their *Bathymodiolus* mussel hosts (Appendix 1). It is therefore highly unlikely that the *Bathymodiolus* SOX symbionts are pathogens that have been mistaken for beneficial symbionts. The remarkably large number of TRGs in the *Bathymodiolus* SOX symbionts bears striking similarity to the arsenal of toxins encoded by *Candidatus* Hamiltonella defensa, which is a facultative symbiont of aphids (*Oliver et al., 2003, 2005*). Both symbionts encode multiple copies of RTX and YD repeat proteins (*Degnan and Moran, 2008*). Ca. Hamiltonella defensa is a defensive symbiont that protects its host from attack by parasitic wasps (*Oliver et al., 2003*). Its protective effect is linked to its complement of RTX and YD repeat toxins (*Degnan and Moran, 2008; Oliver et al., 2010*). Based on our phylogeny of YD repeats, those from the *Bathymodiolus* symbionts cluster together with sequences from both pathogens and beneficial symbionts. Of the beneficial symbionts in this cluster, all except the *Bathymodiolus* symbionts have been shown to produce exotoxins that damage the organisms their host parasitizes, either a plant in the case of *B. rhizoxinica* or an insect in the case of *Photorhabdus/Xenorhabdus* (*Partida-Martinez et al., 2007*). This raises the intriguing possibility that some of the TRGs in the SOX symbiont genomes might function in protecting the mussel hosts against eukaryotic parasites.

Compared to our knowledge of parasitism in shallow-water bivalves, little is known about parasitism in deep-sea *Bathymodiolus* mussels. This is surprising considering that these incredibly dense communities would be ideal habitats for parasites (*Moreira and López-García, 2003*). Two studies have investigated parasitism in *Bathymodiolus* mussels based on the microscopic identification of unusual 'inclusions' in mussel tissues (*Powell et al., 1999; Ward et al., 2004*). The most abundant parasites resembled *Bucephalus*-like trematodes of the phylum Platyhelminthes, which are common in shallow-water mussels (e.g., *Wardle, 1988; Lauenstein et al., 1993; da Silva et al., 2002; Minguez et al., 2012*). *Bucephalus* trematodes infect the gonads of their mussel hosts, which often results in sterilization (*Hopkins, 1957; Coustau et al., 1991*). Like their shallow-water relatives, the *Bucephalus*-like trematodes were abundant in the gonads of *Bathymodiolus childressi* from cold seeps in the Gulf of Mexico (*Powell et al., 1999*). *Powell et al.* (1999) estimated that due to this heavy infection, up to 40% of *B. childressi* populations are reproductively compromised.

The distribution of these trematode parasites has not yet been systematically investigated in *Bathymodiolus*. However, of the three species so far studied, only *B. childressi* was infected by trematodes (*Powell et al., 1999; Ward et al., 2004*). *B. childressi* is one of the few *Bathymodiolus*
identifiable immunity protein, and no type VI secretion system gene was found in any of our draft genomes. The arrangement of YD repeat genes in the *Bathymodiolus* symbiont genomes and the lack of genes encoding type VI secretion systems are therefore inconsistent with a role in competition between closely related symbiont strains.

*Bathymodiolus* mussels are also infected by bacterial intranuclear pathogens called *Candidatus Endonucleobacter bathymodioli*, which are related to the genus *Endozoicomonas* (*Zielinski et al., 2009*). Ca. E. bathymodioli invades the mussel cell nuclei where it multiplies, eventually bursting the infected cell. Intranuclear bacteria are never found in the nuclei of symbiont-containing cells, which led *Zielinski et al.* (2009) to hypothesize that the symbionts can protect their host cells against infection. Consistent with this hypothesis, growth inhibition assays showed that *B. azoricus* gill tissue homogenates inhibit the growth of a broader spectrum of pathogens compared to the symbiont-free foot tissues (*Bettencourt et al., 2007*). The mechanisms of protection against cultured bacterial pathogens and Ca. E. bathymodioli are unknown. There is no evidence from their genomes that the SOX symbionts of *B. sp.* and *B. azoricus* produce antibiotics, but they do have genes for bacteriocin production that were expressed in the six transcriptomes analyzed in this study. Expression of some of the TRGs discovered here, for example, those related to toxins involved in bacteria-bacteria competition, or the production of bacteriocins by symbionts could explain the absence of intranuclear bacteria from symbiont-hosting cells.

**Are some of the TRGs used for protection against eukaryotic parasites?**

Symbionts of other marine invertebrate hosts are able to recognize, enter, and survive within host cells without a large number of TRGs. For example, the genomes of the SOX symbionts of hydrothermal vent *Riftia* tubeworms and the heterotrophic symbionts of whale-fall *Osedax* worms are virtually complete, but they only contain one or a few RTX genes, and no YD repeat genes (*Robidart et al., 2008; Gardebrecht et al., 2012; Goffredi et al., 2014*). There is overwhelming evidence that SOX symbionts are beneficial for their *Bathymodiolus* mussel hosts (Appendix 1). It is therefore highly unlikely that the *Bathymodiolus* SOX symbionts are pathogens that have been mistaken for beneficial symbionts. The remarkably large number of TRGs in the *Bathymodiolus* SOX symbionts bears striking similarity to the arsenal of toxins encoded by *Candidatus* Hamiltonella defensa, which is a facultative symbiont of aphids (*Oliver et al., 2003, 2005*). Both symbionts encode multiple copies of RTX and YD repeat proteins (*Degnan and Moran, 2008*). Ca. Hamiltonella defensa is a defensive symbiont that protects its host from attack by parasitic wasps (*Oliver et al., 2003*). Its protective effect is linked to its complement of RTX and YD repeat toxins (*Degnan and Moran, 2008; Oliver et al., 2010*). Based on our phylogeny of YD repeats, those from the *Bathymodiolus* symbionts cluster together with sequences from both pathogens and beneficial symbionts. Of the beneficial symbionts in this cluster, all except the *Bathymodiolus* symbionts have been shown to produce exotoxins that damage the organisms their host parasitizes, either a plant in the case of *B. rhizoxinica* or an insect in the case of *Photorhabdus/Xenorhabdus* (*Partida-Martinez et al., 2007*). This raises the intriguing possibility that some of the TRGs in the SOX symbiont genomes might function in protecting the mussel hosts against eukaryotic parasites.

Compared to our knowledge of parasitism in shallow-water bivalves, little is known about parasitism in deep-sea *Bathymodiolus* mussels. This is surprising considering that these incredibly dense communities would be ideal habitats for parasites (*Moreira and López-García, 2003*). Two studies have investigated parasitism in *Bathymodiolus* mussels based on the microscopic identification of unusual 'inclusions' in mussel tissues (*Powell et al., 1999; Ward et al., 2004*). The most abundant parasites resembled *Bucephalus*-like trematodes of the phylum Platyhelminthes, which are common in shallow-water mussels (e.g., *Wardle, 1988; Lauenstein et al., 1993; da Silva et al., 2002; Minguez et al., 2012*). *Bucephalus* trematodes infect the gonads of their mussel hosts, which often results in sterilization (*Hopkins, 1957; Coustau et al., 1991*). Like their shallow-water relatives, the *Bucephalus*-like trematodes were abundant in the gonads of *Bathymodiolus childressi* from cold seeps in the Gulf of Mexico (*Powell et al., 1999*). *Powell et al.* (1999) estimated that due to this heavy infection, up to 40% of *B. childressi* populations are reproductively compromised.

The distribution of these trematode parasites has not yet been systematically investigated in *Bathymodiolus*. However, of the three species so far studied, only *B. childressi* was infected by trematodes (*Powell et al., 1999; Ward et al., 2004*). *B. childressi* is one of the few *Bathymodiolus*
species that only associates with MOX symbionts, but not with SOX symbionts. If many of the TRGs encoded by the *Bathymodiolus* SOX symbiont are being used to defend its host against parasites, as is hypothesized for Ca. H. defensa, then this could help to explain why *B. childressi* is so heavily infected by trematodes. The MOX symbionts of *B. azoricus* and *B. childressi* do not encode the abundant TRGs of the *Bathymodiolus* SOX symbiont (Antony CP, personal communication, May 2015).

Our SNP analysis provides further support for the hypothesis that some TRGs may be used for direct beneficial interactions, and some may be used for indirect interactions such as protection against parasites. Genes involved in direct host-symbiont interactions such as recognition and communication are expected to be conserved within the symbiont population of one host species (Jiggins et al., 2002; Bailly et al., 2006). Consistent with this, eight out of nine genes in the MARTX1 region, which we hypothesize may be involved in attachment to the host, do not have a significantly larger number of SNPs per kb compared to the rest of the genome (Figure 6—source data 1). In the ninth gene, annotated as a hypothetical protein, SNPs per kb were slightly above average. In contrast to genes involved in direct host-symbiont interactions, those involved in indirect interactions such as defense against parasites are expected to be highly diverse (see Appendix 2). The large sequence variability in 22 of the TRGs is therefore consistent with a role for these genes in protecting the host against parasites.

## Conclusions

The genomes of the uncultured *Bathymodiolus* SOX symbionts encode a unique arsenal of TRGs, unexpected for a beneficial, nutritional symbiont. We hypothesize that the *Bathymodiolus* SOX symbiont has 'tamed' these genes for use in beneficial interactions with their host. Some of the TRGs may benefit the symbiosis by protecting the symbionts and their hosts from their natural enemies. In most cases, symbionts are either nutritional, that is, their primary role is to provide their host with most or all of its nutrition, or they are defensive (Douglas, 2014). The *Bathymodiolus* symbiont is therefore unusual, as it may play an essential role in both nutrition and defense. The TRGs were most likely acquired by HGT, and this may be the mechanism by which its free-living ancestors acquired the ability to form an intimate relationship with marine animals.

Remarkably, the *Bathymodiolus* SOX symbionts encode a larger complement of these TRGs than any so far sequenced pathogen, suggesting that these 'toxins', although initially discovered in pathogens, may in fact belong to larger protein families that function in both beneficial and pathogenic host-microbe interactions. An alternative to the hypothesis that toxins may be tamed for use in beneficial interactions would be that 'symbiosis factors' may be commandeered for use in harmful interactions. Given our recent recognition of the ubiquity and vast natural diversity of mutualistic interactions between bacteria and eukaryotes (McFall-Ngai et al., 2013), it is possible that many of the genes that are currently annotated as toxins may have first evolved through beneficial host-microbe associations.

## Materials and methods

### Sampling and processing of *Bathymodiolus* mussels

We collected *Bathymodiolus* mussels in Lilliput, Menez Gwen, and Lucky Strike vent sites on the MAR. *Bathymodiolus* sp. from Lilliput on the SMAR were sampled and processed for genome sequencing as in Petersen et al. (2011). For transcriptomics, we sampled mussels from the SMAR at 09°32.85′S, 13°12.64′W. Specimens of *B. azoricus* were collected in three cruise expeditions to Menez Gwen at (i) 37°45.5777′N, 31°38.2611′W during the MOMARETO cruise, (ii) 37°50.68′N, 31°31.17′W during the RV Meteor cruise M82-3, and (iii) Lucky Strike at 37°16′58.5′′N, 32°16′32.2′′W during the Biobaz Cruise. The adductor muscle was dissected from samples of the MOMARETO cruise, while the gill tissue was dissected from mussels collected during the RV Meteor cruise. For samples from the RV Meteor cruise, we used a combination of differential and rate-zonal centrifugation to enrich *Bathymodiolus* SOX symbiont from gill tissue for genomic and proteomic analyses. Samples for transcriptomics were fixed on board in RNAlater (Sigma, Germany) according to the manufacturer instructions and stored at -80°C. An overview of the samples used in this study is shown in *Supplementary file 1E*.
species that only associates with MOX symbionts, but not with SOX symbionts. If many of the TRGs encoded by the *Bathymodiolus* SOX symbiont are being used to defend its host against parasites, as is hypothesized for Ca. H. defensa, then this could help to explain why *B. childressi* is so heavily infected by trematodes. The MOX symbionts of *B. azoricus* and *B. childressi* do not encode the abundant TRGs of the *Bathymodiolus* SOX symbiont (Antony CP, personal communication, May 2015).

Our SNP analysis provides further support for the hypothesis that some TRGs may be used for direct beneficial interactions, and some may be used for indirect interactions such as protection against parasites. Genes involved in direct host-symbiont interactions such as recognition and communication are expected to be conserved within the symbiont population of one host species (Jiggins et al., 2002; Bailly et al., 2006). Consistent with this, eight out of nine genes in the MARTX1 region, which we hypothesize may be involved in attachment to the host, do not have a significantly larger number of SNPs per kb compared to the rest of the genome (Figure 6—source data 1). In the ninth gene, annotated as a hypothetical protein, SNPs per kb were slightly above average. In contrast to genes involved in direct host-symbiont interactions, those involved in indirect interactions such as defense against parasites are expected to be highly diverse (see Appendix 2). The large sequence variability in 22 of the TRGs is therefore consistent with a role for these genes in protecting the host against parasites.

## Conclusions

The genomes of the uncultured *Bathymodiolus* SOX symbionts encode a unique arsenal of TRGs, unexpected for a beneficial, nutritional symbiont. We hypothesize that the *Bathymodiolus* SOX symbiont has 'tamed' these genes for use in beneficial interactions with their host. Some of the TRGs may benefit the symbiosis by protecting the symbionts and their hosts from their natural enemies. In most cases, symbionts are either nutritional, that is, their primary role is to provide their host with most or all of its nutrition, or they are defensive (Douglas, 2014). The *Bathymodiolus* symbiont is therefore unusual, as it may play an essential role in both nutrition and defense. The TRGs were most likely acquired by HGT, and this may be the mechanism by which its free-living ancestors acquired the ability to form an intimate relationship with marine animals.

Remarkably, the *Bathymodiolus* SOX symbionts encode a larger complement of these TRGs than any so far sequenced pathogen, suggesting that these 'toxins', although initially discovered in pathogens, may in fact belong to larger protein families that function in both beneficial and pathogenic host-microbe interactions. An alternative to the hypothesis that toxins may be tamed for use in beneficial interactions would be that 'symbiosis factors' may be commandeered for use in harmful interactions. Given our recent recognition of the ubiquity and vast natural diversity of mutualistic interactions between bacteria and eukaryotes (McFall-Ngai et al., 2013), it is possible that many of the genes that are currently annotated as toxins may have first evolved through beneficial host-microbe associations.

## Materials and methods

### Sampling and processing of *Bathymodiolus* mussels

We collected *Bathymodiolus* mussels in Lilliput, Menez Gwen, and Lucky Strike vent sites on the MAR. *Bathymodiolus* sp. from Lilliput on the SMAR were sampled and processed for genome sequencing as in Petersen et al. (2011). For transcriptomics, we sampled mussels from the SMAR at 09°32.85′S, 13°12.64′W. Specimens of *B. azoricus* were collected in three cruise expeditions to Menez Gwen at (i) 37°45.5777′N, 31°38.2611′W during the MOMARETO cruise, (ii) 37°50.68′N, 31°31.17′W during the RV Meteor cruise M82-3, and (iii) Lucky Strike at 37°16′58.5′′N, 32°16′32.2′′W during the Biobaz Cruise. The adductor muscle was dissected from samples of the MOMARETO cruise, while the gill tissue was dissected from mussels collected during the RV Meteor cruise. For samples from the RV Meteor cruise, we used a combination of differential and rate-zonal centrifugation to enrich *Bathymodiolus* SOX symbiont from gill tissue for genomic and proteomic analyses. Samples for transcriptomics were fixed on board in RNAlater (Sigma, Germany) according to the manufacturer instructions and stored at -80°C. An overview of the samples used in this study is shown in *Supplementary file 1E*.
## DNA extraction, sequencing, genome assembly and binning

DNA extraction and genome sequencing of the B. sp. SMAR SOX symbiont was described in *Petersen et al. (2011)*. Briefly, the gill tissue of a single individual was ground in a glass tissue homogenizer and frozen until further processing. In the home laboratory, the homogenate was diluted in phosphate-buffered saline (PBS)1x, centrifuged at 400×g for one minute, and the supernatant filtered through a 12-µm GTTP filter (Millipore, Germany). Centrifugation and filtration was repeated 20 times. The filtrates were passed through GTTP filters of 8 µm, 5 µm, 3 µm, and 2 µm. Cells collected on the 0.2-µm filter were used for DNA extraction after *Zhou et al. (1996)* with an initial incubation overnight at 37°C in extraction buffer and proteinase K. 6 kb mate-paired reads were sequenced with 454-Titanium and 36 bp Illumina reads. 454 reads were assembled with Newbler v2.3 (454 Life Sciences Corporation) and 569 pyrosequencing errors were corrected using the Illumina reads.

DNA from the *B. azoricus* SOX symbiont enrichments from three individuals (gradient pellets, see Appendix 4) was extracted according to *Zhou et al. (1996)*. Genomic DNA was extracted from adductor muscle using a CTAB/PVP extraction procedure (2% CTAB, 1% PVP, 1.4 M NaCl, 0.2% beta-mercaptoethanol, 100 mM Tris HCl pH 8, 0.1 mg ml⁻¹ proteinase K). After complete digestion of tissues (1 hr at 60°C), the mixture was incubated with 1 µl of RNase for 30 min at 37°C. An equal volume of chloroform-isoamyl alcohol (24:1) was added and tubes were slowly mixed by inversion for 3 min before a 10 min centrifugation at 14,000 rpm and 4°C. The supernatant was collected in a fresh tube, and DNA was precipitated with 2/3 volume of cold isopropanol (1 hr at -20°C). The DNA pellet was recovered by centrifugation (14,000 rpm at 4°C for 20 min), washed with 75% cold ethanol, air-dried, and suspended in 100 µl of sterile water. 454 sequencing was done by Genoscope to sequence the gradient pellet from gill tissue, and by OIST to sequence the adductor muscle of *B. azoricus*. For the gradient pellet, a 3 kb insert 454 library was prepared according to manufacturer protocols for mate-pair sequencing. 630752 reads were generated on a Titanium FLX sequencing machine and assembled using Newbler software (version 08172012). 1310 contigs larger than 500 bp were obtained, forming 130 scaffolds of a total length of 1668565 bp. The assembly from the adductor muscle was done with Newbler v. 2.7 (454 Life Sciences Corporation) as described in *Takeuchi et al. (2012)* resulting in 644000 contigs of a total length of 510449434 bp. The adductor muscle and gradient pellet metagenomes of *B. azoricus* were binned to separate the SOX symbiont from the MOX symbiont and host genomes with Metawatt V. 1.7, which uses tetranucleotide frequencies, coverage, GC content, and taxonomic information for binning (*Strous et al., 2012*). Only sequences longer than 800 bp were considered for further analyses. Since we could only recover an 829 bp fragment of the 16S rRNA from BazSymB, the same DNA that was used for 454 sequencing was used as template for PCR amplification with the universal primers GM3f/GM4r (*Muyzer et al., 1995*). The PCR product was directly sequenced with Sanger and assembled using Geneious V7 (*Kearse et al., 2012*). The 16S rRNA of BazSymB can be found under the accession number (LN871183).

## Genome annotation

We annotated the genomes of the *Bathymodiolus* symbionts (BspSym: PRJNA65421, BazSymA: PRJEB8263, and BazSymB: PRJEB8264), *Candidatus* V. okutanii (NC_009465), *Candidatus* R. magnifica (CP000488), and SUP05 metagenome (ACSG01000000; GQ351266 to GQ351269 and GQ369726) with the following workflow: we used Glimmer (*Delcher et al., 2007*) for open reading frame (ORF) prediction. Ribosomal RNA genes were detected with RNAmmer (*Lagesen et al., 2007*) and tRNAs with tRNAscan-SE (*Lowe and Eddy, 1997*). Annotation was done with GenDB 2.4 (*Meyer, 2003*) and supplemented by JCoast 1.7 (*Richter et al., 2008*) to integrate the results of BLASTp (cut-off e-value of 10.0) against sequence databases NCBI-nr (*Altschul et al., 1997*), SwissProt (*Boeckmann, 2003*), KEGG (*Kanehisa et al., 2011*), COG (*Tatusov, 2000*), Pfam (*Bateman et al., 2004*), and InterPro (*Hunter et al., 2009*). TMHMM (*Krogh et al., 2001*) was used for transmembrane helix analysis and SignalP (*Emanuelsson et al., 2007*) for signal peptide predictions. Sequences of the cytochrome c oxidase subunit 1 of the *Bathymodiolus* mussels were submitted to the European Nucleotide Archive when available (*Supplementary file 1E*).

## PCR amplification of regions with lack of synteny

We designed primers to amplify four regions covering a region with lack of synteny. The primer sequences and annealing temperatures are listed in *Supplementary file 1F*. The PCR program
## DNA extraction, sequencing, genome assembly and binning

DNA extraction and genome sequencing of the B. sp. SMAR SOX symbiont was described in *Petersen et al. (2011)*. Briefly, the gill tissue of a single individual was ground in a glass tissue homogenizer and frozen until further processing. In the home laboratory, the homogenate was diluted in phosphate-buffered saline (PBS)1x, centrifuged at 400×g for one minute, and the supernatant filtered through a 12-µm GTTP filter (Millipore, Germany). Centrifugation and filtration was repeated 20 times. The filtrates were passed through GTTP filters of 8 µm, 5 µm, 3 µm, and 2 µm. Cells collected on the 0.2-µm filter were used for DNA extraction after *Zhou et al. (1996)* with an initial incubation overnight at 37°C in extraction buffer and proteinase K. 6 kb mate-paired reads were sequenced with 454-Titanium and 36 bp Illumina reads. 454 reads were assembled with Newbler v2.3 (454 Life Sciences Corporation) and 569 pyrosequencing errors were corrected using the Illumina reads.

DNA from the *B. azoricus* SOX symbiont enrichments from three individuals (gradient pellets, see Appendix 4) was extracted according to *Zhou et al. (1996)*. Genomic DNA was extracted from adductor muscle using a CTAB/PVP extraction procedure (2% CTAB, 1% PVP, 1.4 M NaCl, 0.2% beta-mercaptoethanol, 100 mM Tris HCl pH 8, 0.1 mg ml⁻¹ proteinase K). After complete digestion of tissues (1 hr at 60°C), the mixture was incubated with 1 µl of RNase for 30 min at 37°C. An equal volume of chloroform-isoamyl alcohol (24:1) was added and tubes were slowly mixed by inversion for 3 min before a 10 min centrifugation at 14,000 rpm and 4°C. The supernatant was collected in a fresh tube, and DNA was precipitated with 2/3 volume of cold isopropanol (1 hr at -20°C). The DNA pellet was recovered by centrifugation (14,000 rpm at 4°C for 20 min), washed with 75% cold ethanol, air-dried, and suspended in 100 µl of sterile water. 454 sequencing was done by Genoscope to sequence the gradient pellet from gill tissue, and by OIST to sequence the adductor muscle of *B. azoricus*. For the gradient pellet, a 3 kb insert 454 library was prepared according to manufacturer protocols for mate-pair sequencing. 630752 reads were generated on a Titanium FLX sequencing machine and assembled using Newbler software (version 08172012). 1310 contigs larger than 500 bp were obtained, forming 130 scaffolds of a total length of 1668565 bp. The assembly from the adductor muscle was done with Newbler v. 2.7 (454 Life Sciences Corporation) as described in *Takeuchi et al. (2012)* resulting in 644000 contigs of a total length of 510449434 bp. The adductor muscle and gradient pellet metagenomes of *B. azoricus* were binned to separate the SOX symbiont from the MOX symbiont and host genomes with Metawatt V. 1.7, which uses tetranucleotide frequencies, coverage, GC content, and taxonomic information for binning (*Strous et al., 2012*). Only sequences longer than 800 bp were considered for further analyses. Since we could only recover an 829 bp fragment of the 16S rRNA from BazSymB, the same DNA that was used for 454 sequencing was used as template for PCR amplification with the universal primers GM3f/GM4r (*Muyzer et al., 1995*). The PCR product was directly sequenced with Sanger and assembled using Geneious V7 (*Kearse et al., 2012*). The 16S rRNA of BazSymB can be found under the accession number (LN871183).

## Genome annotation

We annotated the genomes of the *Bathymodiolus* symbionts (BspSym: PRJNA65421, BazSymA: PRJEB8263, and BazSymB: PRJEB8264), *Candidatus* V. okutanii (NC_009465), *Candidatus* R. magnifica (CP000488), and SUP05 metagenome (ACSG01000000; GQ351266 to GQ351269 and GQ369726) with the following workflow: we used Glimmer (*Delcher et al., 2007*) for open reading frame (ORF) prediction. Ribosomal RNA genes were detected with RNAmmer (*Lagesen et al., 2007*) and tRNAs with tRNAscan-SE (*Lowe and Eddy, 1997*). Annotation was done with GenDB 2.4 (*Meyer, 2003*) and supplemented by JCoast 1.7 (*Richter et al., 2008*) to integrate the results of BLASTp (cut-off e-value of 10.0) against sequence databases NCBI-nr (*Altschul et al., 1997*), SwissProt (*Boeckmann, 2003*), KEGG (*Kanehisa et al., 2011*), COG (*Tatusov, 2000*), Pfam (*Bateman et al., 2004*), and InterPro (*Hunter et al., 2009*). TMHMM (*Krogh et al., 2001*) was used for transmembrane helix analysis and SignalP (*Emanuelsson et al., 2007*) for signal peptide predictions. Sequences of the cytochrome c oxidase subunit 1 of the *Bathymodiolus* mussels were submitted to the European Nucleotide Archive when available (*Supplementary file 1E*).

## PCR amplification of regions with lack of synteny

We designed primers to amplify four regions covering a region with lack of synteny. The primer sequences and annealing temperatures are listed in *Supplementary file 1F*. The PCR program
consisted of an initial denaturation step of 98°C for 30 s, followed by 35 cycles at 98°C for 10 s, specific annealing temperature for 30 s, 72°C for 2 min, and a final extension at 72°C for 10 min. We obtained a PCR product of the expected size based on our assembly of all four targeted regions. We sequenced one of these by Sanger sequencing using ABI BigDye v3.1 and the ABI PRISM 3100 genetic analyzer (Applied Biosystems, Foster City, CA).

## Genome analysis and comparison to close relatives

We used CheckM to evaluate the completeness of our draft genomes with a set of single-copy marker genes that are specific to proteobacteria (lineage-specific marker set of CheckM p_Proteobacteria, UID3880) (Parks et al., 2015). We estimated the similarity of BazSymB and BazSymA draft genomes with the mean and standard deviation of genes with bi-directional best BLAST hits. The initial comparison of gene content between the clam symbionts, mussel symbionts draft genomes, and SUP05 metagenome was done by BLAST score ratio (BSR) with a BSR cut-off of 0.4 (Rasko et al., 2005). Since toxin genes are expected to have a higher mutation rate, we compared the toxin distribution among the three mussel symbiont draft genomes and their closest relatives with a BSR cut-off of 0.2.

Because Ca. R. magnifica is the largest vesicomyid genome available, a whole genome comparison with SUP05 was done using a Dotplot produced by Ugene (http://ugene.unipro.ru) with a minimum length of 50 bp and 90% similarity. Since a high gene synteny could be observed, Ca. R. magnifica was used as scaffold to reorder the contigs from SUP05 and the B. sp. symbiont using Mauve (Rissman et al., 2009). We did not do this analysis for the B. azoricus symbiont genomes because they were highly fragmented. To compare the influence of HGT on genome evolution, we used the method of Davis and Olsen (2010) to identify genes with different codon usage patterns, indicating that they may have been recently acquired via HGT.

## Identification, classification, and structural analysis of TRGs

The initial genome annotation contained many genes annotated as toxins (see Appendix 4 for details of the automatic annotation). Most were related to toxins from three broad previously defined classes, the YD repeats, RTX, and MARTX toxins. Most of the YD repeat and RTX genes could be identified accurately by automatic annotation due to the presence of signature repeat regions. Some of the predicted YD genes appeared to be truncated. Alignment with the closest hits was used to extend the partial ORFs that were not correctly predicted. To curate the annotation of MARTX genes, we built Hidden Markov Model (HMM) profiles for the characterized MARTX domains in Vibrio cholerae: actin cross-linking domain, Rho GTPase inactivation, and cysteine protease domain (Satchell, 2011). The profiles were used to scan the SOX symbiont genomes with HMMER (Finn et al., 2011). We searched for functional domains in the proteins identified as TRGs with SMART (Ponting et al., 1999). To analyze the diversity and redundancy of the TRGs in the SOX genomes, we constructed a protein similarity network. All protein sequences of the SOX symbionts were searched against Uniprot with BLAST (coverage >50%, similarity >25%, e-value < 10⁻⁵). We also used BLAST of all against all sequences of the Bathymodiolus SOX symbiont to recover partial TRGs that could not be identified by automatic annotation. The two searches were combined to produce a sequence similarity network based on transitivity clustering using the plugin Blast2SimilarityGraph in Cytoscape (Srinivasan and Moon, 1999; Shannon et al., 2003; Wittkop et al., 2010). Only the sub-networks that were connected in a maximum of two steps to the TRGs were considered. All TRGs of the SOX symbionts were submitted to Phyre2 (Kelley and Sternberg, 2009) to predict the secondary structure of the protein. We looked for clusters in the protein similarity network that could be associated with a subunit of a toxin complex or an active domain.

## Search for TRGs in SUP05

None of the three toxin-related classes has characteristic patterns, profiles, or domains that can be searched with standard tools. Often, their only shared feature is a repeat region, such as the RTX repeat in RTX and MARTX, which is a calcium-binding site containing G and N residues. However, the number of G and N residues and the length of the repeat is highly variable, which makes profile searches impossible. Moreover, these repeat regions can be shared by other calcium-binding proteins such as integrins and fibronectin, which do not act as toxins. The sequence and length of the YD
consisted of an initial denaturation step of 98°C for 30 s, followed by 35 cycles at 98°C for 10 s, specific annealing temperature for 30 s, 72°C for 2 min, and a final extension at 72°C for 10 min. We obtained a PCR product of the expected size based on our assembly of all four targeted regions. We sequenced one of these by Sanger sequencing using ABI BigDye v3.1 and the ABI PRISM 3100 genetic analyzer (Applied Biosystems, Foster City, CA).

## Genome analysis and comparison to close relatives

We used CheckM to evaluate the completeness of our draft genomes with a set of single-copy marker genes that are specific to proteobacteria (lineage-specific marker set of CheckM p_Proteobacteria, UID3880) (Parks et al., 2015). We estimated the similarity of BazSymB and BazSymA draft genomes with the mean and standard deviation of genes with bi-directional best BLAST hits. The initial comparison of gene content between the clam symbionts, mussel symbionts draft genomes, and SUP05 metagenome was done by BLAST score ratio (BSR) with a BSR cut-off of 0.4 (Rasko et al., 2005). Since toxin genes are expected to have a higher mutation rate, we compared the toxin distribution among the three mussel symbiont draft genomes and their closest relatives with a BSR cut-off of 0.2.

Because Ca. R. magnifica is the largest vesicomyid genome available, a whole genome comparison with SUP05 was done using a Dotplot produced by Ugene (http://ugene.unipro.ru) with a minimum length of 50 bp and 90% similarity. Since a high gene synteny could be observed, Ca. R. magnifica was used as scaffold to reorder the contigs from SUP05 and the B. sp. symbiont using Mauve (Rissman et al., 2009). We did not do this analysis for the B. azoricus symbiont genomes because they were highly fragmented. To compare the influence of HGT on genome evolution, we used the method of Davis and Olsen (2010) to identify genes with different codon usage patterns, indicating that they may have been recently acquired via HGT.

## Identification, classification, and structural analysis of TRGs

The initial genome annotation contained many genes annotated as toxins (see Appendix 4 for details of the automatic annotation). Most were related to toxins from three broad previously defined classes, the YD repeats, RTX, and MARTX toxins. Most of the YD repeat and RTX genes could be identified accurately by automatic annotation due to the presence of signature repeat regions. Some of the predicted YD genes appeared to be truncated. Alignment with the closest hits was used to extend the partial ORFs that were not correctly predicted. To curate the annotation of MARTX genes, we built Hidden Markov Model (HMM) profiles for the characterized MARTX domains in Vibrio cholerae: actin cross-linking domain, Rho GTPase inactivation, and cysteine protease domain (Satchell, 2011). The profiles were used to scan the SOX symbiont genomes with HMMER (Finn et al., 2011). We searched for functional domains in the proteins identified as TRGs with SMART (Ponting et al., 1999). To analyze the diversity and redundancy of the TRGs in the SOX genomes, we constructed a protein similarity network. All protein sequences of the SOX symbionts were searched against Uniprot with BLAST (coverage >50%, similarity >25%, e-value < 10⁻⁵). We also used BLAST of all against all sequences of the Bathymodiolus SOX symbiont to recover partial TRGs that could not be identified by automatic annotation. The two searches were combined to produce a sequence similarity network based on transitivity clustering using the plugin Blast2SimilarityGraph in Cytoscape (Srinivasan and Moon, 1999; Shannon et al., 2003; Wittkop et al., 2010). Only the sub-networks that were connected in a maximum of two steps to the TRGs were considered. All TRGs of the SOX symbionts were submitted to Phyre2 (Kelley and Sternberg, 2009) to predict the secondary structure of the protein. We looked for clusters in the protein similarity network that could be associated with a subunit of a toxin complex or an active domain.

## Search for TRGs in SUP05

None of the three toxin-related classes has characteristic patterns, profiles, or domains that can be searched with standard tools. Often, their only shared feature is a repeat region, such as the RTX repeat in RTX and MARTX, which is a calcium-binding site containing G and N residues. However, the number of G and N residues and the length of the repeat is highly variable, which makes profile searches impossible. Moreover, these repeat regions can be shared by other calcium-binding proteins such as integrins and fibronectin, which do not act as toxins. The sequence and length of the YD
repeat is similarly variable. To identify homologs of the *Bathymodiolus* symbiont TRGs in published metagenomes and metatranscriptomes, we therefore used BLASTp (coverage >50%, similarity >25%, e-value < 0.001) (*Supplementary file 1G*).

Best hits were blasted against the SOX symbiont proteins to search for signatures of the TRGs. We looked for bacterial genomes that had similar TRGs to the *Bathymodiolus* SOX symbiont, using the Integrated Microbial Genomes (IMG) database (Markowitz et al., 2011) with BLASTp (>50% coverage, > 25% similarity, e-value 0.001). 122 genomes were retrieved and manually curated for the potential of pathogenicity and biofilm formation, as well as the lifestyle categories extracellular host-associated, intracellular, and free-living bacteria. Kruskall–Wallis and Mann–Whitney–Wilcoxon tests were used to estimate if the increased number of a class of TR is biased in a certain lifestyle. All statistical analyses were done in R.

## Phylogenetic analyses

The 16S rRNA gene sequences of the *Bathymodiolus* symbionts were imported into the Silva database (release Ref 119) and initially aligned with the SINA aligner (Pruesse et al., 2012). The final alignment was refined with MAFFT (Katoh et al., 2002). A maximum likelihood tree was estimated from the alignment of 1653 nucleotide positions using RaxML with 100 bootstrap replicates.

To construct the YD phylogeny, the data set from BspSym was used to obtain related sequences from GenBank and *B. azoricus*. Because two metagenomes of *B. azoricus* were analyzed, we used CD-hit to remove redundancy at 100% similarity (Li and Godzik, 2006). YD repeat proteins are very variable in the C and N terminus. Therefore, the selection criterion we used was the presence of the conserved *rhs* domain. Only this 'core' region was used for phylogeny (Jackson et al., 2009). YD proteins were aligned with MAFFT with BLOSUM30 (Katoh et al., 2002). Phylogenetic analyses were done in ARB with maximum likelihood and bayesian reconstructions using a filter of 10% similarity, which resulted in 536 amino acid positions (Ronquist and Huelsenbeck, 2003; Ludwig et al., 2004; Stamatakis, 2006). Bootstrap support was calculated with 100 replicates in maximum likelihood. MrBayes was run for 9 million generations and two independent runs of four heated chains. A consensus tree of both methods was constructed. Polytomies were introduced when both methods did not agree.

## Statistical analyses of TRGs content in microbial genomes

Genomes with curated metadata are available through IMG (Markowitz et al., 2011). Nevertheless, most entries contained no information in categories we were interested in such as biofilm formation and intracellular lifestyle. IMG was used as starting point to search for curated genomes that had genes with similarity to the TRGs of the mussel symbionts (>10% similarity, e-value < e⁻⁵). The protein sequences of these genomes were retrieved from NCBI based on the genome name (modified Python script from Sixing Huang, Max Planck Institute Bremen). BLASTp was used to retrieve protein sequences related to the SOX TRGs (>50% coverage, > 25% similarity, e-value of 0.001). 122 genomes had at least one gene that was similar to at least one TRG. We considered the genomes of different strains as independent events for statistical analysis, as even closely related strains can have different lifestyles. These genomes were manually curated based on a literature search and classified according to the following categories: (1) lifestyle: divided into three sub-categories (a) associated: if the organism is at any stage host-associated (based on IMG metadata), (b) intracellular: includes obligate and facultative symbionts (searched in Google scholar with the keywords 'intracellular bacteria' and references read for more details when the abstract was not sufficient), (c) free-living: bacteria that are not host-associated and not intracellular; (2) pathogen: bacteria that can produce disease (information obtained from IMG); (3) biofilm: bacteria found in biofilms (Google scholar keywords 'biofilm' and 'microbial mat' with literature analyses when not clear).

The data were formatted and merged with self-written Perl scripts and *R Development Core Team* (2011). Genomes with similar toxins or TRGs to the *Bathymodiolus* symbionts are shown in *Supplementary file 1B*. The sum of the number of TRGs belonging to YD, RTX, and MARTX was normalized with the total gene count for each genome and multiplied by a factor of 1000. To compare the number of TRGs against the lifestyle categories, we used Kruskall–Wallis test. A post hoc analysis was carried out on significant p-values for the associated category with the Mann–Whitney–Wilcoxon test. Statistical analyses were done in R. We tested whether the bacteria were enriched in any of the
repeat is similarly variable. To identify homologs of the *Bathymodiolus* symbiont TRGs in published metagenomes and metatranscriptomes, we therefore used BLASTp (coverage >50%, similarity >25%, e-value < 0.001) (*Supplementary file 1G*).

Best hits were blasted against the SOX symbiont proteins to search for signatures of the TRGs. We looked for bacterial genomes that had similar TRGs to the *Bathymodiolus* SOX symbiont, using the Integrated Microbial Genomes (IMG) database (Markowitz et al., 2011) with BLASTp (>50% coverage, > 25% similarity, e-value 0.001). 122 genomes were retrieved and manually curated for the potential of pathogenicity and biofilm formation, as well as the lifestyle categories extracellular host-associated, intracellular, and free-living bacteria. Kruskall–Wallis and Mann–Whitney–Wilcoxon tests were used to estimate if the increased number of a class of TR is biased in a certain lifestyle. All statistical analyses were done in R.

## Phylogenetic analyses

The 16S rRNA gene sequences of the *Bathymodiolus* symbionts were imported into the Silva database (release Ref 119) and initially aligned with the SINA aligner (Pruesse et al., 2012). The final alignment was refined with MAFFT (Katoh et al., 2002). A maximum likelihood tree was estimated from the alignment of 1653 nucleotide positions using RaxML with 100 bootstrap replicates.

To construct the YD phylogeny, the data set from BspSym was used to obtain related sequences from GenBank and *B. azoricus*. Because two metagenomes of *B. azoricus* were analyzed, we used CD-hit to remove redundancy at 100% similarity (Li and Godzik, 2006). YD repeat proteins are very variable in the C and N terminus. Therefore, the selection criterion we used was the presence of the conserved *rhs* domain. Only this 'core' region was used for phylogeny (Jackson et al., 2009). YD proteins were aligned with MAFFT with BLOSUM30 (Katoh et al., 2002). Phylogenetic analyses were done in ARB with maximum likelihood and bayesian reconstructions using a filter of 10% similarity, which resulted in 536 amino acid positions (Ronquist and Huelsenbeck, 2003; Ludwig et al., 2004; Stamatakis, 2006). Bootstrap support was calculated with 100 replicates in maximum likelihood. MrBayes was run for 9 million generations and two independent runs of four heated chains. A consensus tree of both methods was constructed. Polytomies were introduced when both methods did not agree.

## Statistical analyses of TRGs content in microbial genomes

Genomes with curated metadata are available through IMG (Markowitz et al., 2011). Nevertheless, most entries contained no information in categories we were interested in such as biofilm formation and intracellular lifestyle. IMG was used as starting point to search for curated genomes that had genes with similarity to the TRGs of the mussel symbionts (>10% similarity, e-value < e⁻⁵). The protein sequences of these genomes were retrieved from NCBI based on the genome name (modified Python script from Sixing Huang, Max Planck Institute Bremen). BLASTp was used to retrieve protein sequences related to the SOX TRGs (>50% coverage, > 25% similarity, e-value of 0.001). 122 genomes had at least one gene that was similar to at least one TRG. We considered the genomes of different strains as independent events for statistical analysis, as even closely related strains can have different lifestyles. These genomes were manually curated based on a literature search and classified according to the following categories: (1) lifestyle: divided into three sub-categories (a) associated: if the organism is at any stage host-associated (based on IMG metadata), (b) intracellular: includes obligate and facultative symbionts (searched in Google scholar with the keywords 'intracellular bacteria' and references read for more details when the abstract was not sufficient), (c) free-living: bacteria that are not host-associated and not intracellular; (2) pathogen: bacteria that can produce disease (information obtained from IMG); (3) biofilm: bacteria found in biofilms (Google scholar keywords 'biofilm' and 'microbial mat' with literature analyses when not clear).

The data were formatted and merged with self-written Perl scripts and R Development Core Team (2011). Genomes with similar toxins or TRGs to the *Bathymodiolus* symbionts are shown in *Supplementary file 1B*. The sum of the number of TRGs belonging to YD, RTX, and MARTX was normalized with the total gene count for each genome and multiplied by a factor of 1000. To compare the number of TRGs against the lifestyle categories, we used Kruskall–Wallis test. A post hoc analysis was carried out on significant p-values for the associated category with the Mann–Whitney–Wilcoxon test. Statistical analyses were done in R. We tested whether the bacteria were enriched in any of the
three toxin-related classes at the class, order, or family by using one-way Permanova with 9999 permutations and Euclidean distances. p-values were corrected with Bonferroni correction for multiple testing.

## Transcriptomics

To extract the total RNA of three individuals of *B. azoricus*, the gill tissue was incubated overnight in RNAlater (Sigma) at 4°C. A fragment of the gill was dissected and homogenized. RNA was extracted with RNeasy Plus MicroKit (Qiagen, Hilden, Germany) according to the manufacturer's instructions. To remove cell debris and to improve RNA yield, we used QIAshredder Mini Spin Columns (Qiagen). The quality of the RNA was assessed with Agilent 2100 Bioanalyzer. The RNA was used for cDNA synthesis with the Ovation RNA-Seq System V2 (NuGEN, San Carlos, CA). To extract the total RNA of three individuals of *B. sp.*, the gill tissue was placed separately on liquid nitrogen, homogenized, and stored overnight at 4°C in self-made RNAlater (10 mM ethylenediaminetetraacetic acid (EDTA), 25 mM tri-sodiumcitrate-dihydrate, 5.3 M ammonium sulfate, adjusted to pH 5.2). After removal of RNAlater, samples were incubated in 600 µl RLT-β-mercaptoethanol buffer (1:100) for 10 min and homogenized on QIAshredder columns (Qiagen). Total RNA was extracted with RNeasy Mini Kit (Qiagen). We applied the DNA-free DNAse Treatment and Removal Kit according to the manufacturer's instructions (Invitrogen, Carlsbad, CA/Ambion, Austin, TX). RNA quality was checked on the Experion Automated Electrophoresis Station using the RNA StdSens Analysis protocol (BioRad, Hercules, CA).

Libraries of *B. spp.* were generated with the Illumina TruSeq RNA Sample Preparation Kit and sequenced 2 x 100 paired-end on an Illumina HiSeq 2000 platform at the Institute of Clinical Molecular Biology (Kiel). A total of 32.9, 38.2, and 38.4 million reads were sequenced per individual of *B. azoricus*. Libraries of *B. azoricus* were generated with DNA library prep kit for Illumina (BioLABS, Frankfurt am Main, Germany) and sequenced as single 100-bp reads on an Illumina HiSeq 2500 platform at the Max Planck Genome Centre (Cologne). A total of 4.3, 4.8, and 6.9 million reads were sequenced per individual of *B. azoricus*. Adaptor removal and quality trimming was done with Nesoni (http://www.vicbioinformatics.com/software.nesoni.shtml) using a quality threshold of 20. To remove ribosomal sequences from the data, we mapped the reads against the SILVA 115 SSU database with Bowtie2 and kept those reads that failed to align (*Langmead and Salzberg, 2012; Quast et al., 2013*). The abundance of the transcripts per gene was estimated with Rockhopper that uses upper quartile normalization (*McClure et al., 2013*). The expression values of the TRGs were normalized to the expression of RubisCO (*Supplementary file 1D*). Transcriptome reads that mapped to the SOX symbionts with Bbmap (http://bbmap.sourceforge.net/) were deposited in the European Nucleotide Archive under the accession numbers PRJEB7941 for *B. azoricus* and PRJEB7943 for *B. sp*.

## SNP analysis

We calculated SNPs to compare substitution rates of genes in the SOX symbiont genomes using three transcriptomes of *B. sp* and three of *B. azoricus* (*Supplementary file 1E*). Reads were normalized to a coverage of maximum 200 and minimum of five with BBNorm v33 (http://sourceforge.net/projects/bbmap/). Transcriptomes of *Bathymodiolus* sp. were mapped against the draft genome BspSym and transcriptomes of *B. azoricus* were mapped to the draft genomes of BazSymA and BazSymB with BBMap v33. We only considered those reads that mapped to a single position in the reference genome and that had higher than 90% identity alignments. SNPs were called independently for the draft genomes BspSym, BazSymA, and BazSymB with the Genome Analysis ToolKit as described by *De Wit et al. (2012)* with some modifications (*McKenna et al., 2010*). In summary, regions needing realignment were identified and realigned over intervals. SNPs and insertions or deletions (InDels) were called with the haplotype caller with a minimum confidence of 20. A filter was applied around InDels with a mask extension of 5. SNPs per gene were obtained with BEDTools (*Quinlan and Hall, 2010*). The number of SNPs per gene was normalized according to the gene length minus regions of unknown sequence for genes containing Ns. We did not consider genes shorter than 150 bp nucleotides or outlier genes at the ends of scaffolds that had an unusual number of SNPs.

## Proteomics

Soluble proteins were extracted from SOX symbiont enrichments, host cytosolic fractions, and whole gill and foot tissue in biological duplicates (Appendix 4). All proteome samples were obtained from the M82-3 cruise. Membrane proteins were extracted from the SOX symbiont enrichments and the
three toxin-related classes at the class, order, or family by using one-way Permanova with 9999 permutations and Euclidean distances. p-values were corrected with Bonferroni correction for multiple testing.

## Transcriptomics

To extract the total RNA of three individuals of *B. azoricus*, the gill tissue was incubated overnight in RNAlater (Sigma) at 4°C. A fragment of the gill was dissected and homogenized. RNA was extracted with RNeasy Plus MicroKit (Qiagen, Hilden, Germany) according to the manufacturer's instructions. To remove cell debris and to improve RNA yield, we used QIAshredder Mini Spin Columns (Qiagen). The quality of the RNA was assessed with Agilent 2100 Bioanalyzer. The RNA was used for cDNA synthesis with the Ovation RNA-Seq System V2 (NuGEN, San Carlos, CA). To extract the total RNA of three individuals of *B. sp.*, the gill tissue was placed separately on liquid nitrogen, homogenized, and stored overnight at 4°C in self-made RNAlater (10 mM ethylenediaminetetraacetic acid (EDTA), 25 mM tri-sodiumcitrate-dihydrate, 5.3 M ammonium sulfate, adjusted to pH 5.2). After removal of RNAlater, samples were incubated in 600 µl RLT-β-mercaptoethanol buffer (1:100) for 10 min and homogenized on QIAshredder columns (Qiagen). Total RNA was extracted with RNeasy Mini Kit (Qiagen). We applied the DNA-free DNAse Treatment and Removal Kit according to the manufacturer's instructions (Invitrogen, Carlsbad, CA/Ambion, Austin, TX). RNA quality was checked on the Experion Automated Electrophoresis Station using the RNA StdSens Analysis protocol (BioRad, Hercules, CA).

Libraries of *B. spp.* were generated with the Illumina TruSeq RNA Sample Preparation Kit and sequenced 2 x 100 paired-end on an Illumina HiSeq 2000 platform at the Institute of Clinical Molecular Biology (Kiel). A total of 32.9, 38.2, and 38.4 million reads were sequenced per individual of *B. azoricus*. Libraries of *B. azoricus* were generated with DNA library prep kit for Illumina (BioLABS, Frankfurt am Main, Germany) and sequenced as single 100-bp reads on an Illumina HiSeq 2500 platform at the Max Planck Genome Centre (Cologne). A total of 4.3, 4.8, and 6.9 million reads were sequenced per individual of *B. azoricus*. Adaptor removal and quality trimming was done with Nesoni (http://www.vicbioinformatics.com/software.nesoni.shtml) using a quality threshold of 20. To remove ribosomal sequences from the data, we mapped the reads against the SILVA 115 SSU database with Bowtie2 and kept those reads that failed to align (*Langmead and Salzberg, 2012; Quast et al., 2013*). The abundance of the transcripts per gene was estimated with Rockhopper that uses upper quartile normalization (*McClure et al., 2013*). The expression values of the TRGs were normalized to the expression of RubisCO (*Supplementary file 1D*). Transcriptome reads that mapped to the SOX symbionts with Bbmap (http://bbmap.sourceforge.net/) were deposited in the European Nucleotide Archive under the accession numbers PRJEB7941 for *B. azoricus* and PRJEB7943 for *B. sp*.

## SNP analysis

We calculated SNPs to compare substitution rates of genes in the SOX symbiont genomes using three transcriptomes of *B. sp* and three of *B. azoricus* (*Supplementary file 1E*). Reads were normalized to a coverage of maximum 200 and minimum of five with BBNorm v33 (http://sourceforge.net/projects/bbmap/). Transcriptomes of *Bathymodiolus* sp. were mapped against the draft genome BspSym and transcriptomes of *B. azoricus* were mapped to the draft genomes of BazSymA and BazSymB with BBMap v33. We only considered those reads that mapped to a single position in the reference genome and that had higher than 90% identity alignments. SNPs were called independently for the draft genomes BspSym, BazSymA, and BazSymB with the Genome Analysis ToolKit as described by *De Wit et al. (2012)* with some modifications (*McKenna et al., 2010*). In summary, regions needing realignment were identified and realigned over intervals. SNPs and insertions or deletions (InDels) were called with the haplotype caller with a minimum confidence of 20. A filter was applied around InDels with a mask extension of 5. SNPs per gene were obtained with BEDTools (*Quinlan and Hall, 2010*). The number of SNPs per gene was normalized according to the gene length minus regions of unknown sequence for genes containing Ns. We did not consider genes shorter than 150 bp nucleotides or outlier genes at the ends of scaffolds that had an unusual number of SNPs.

## Proteomics

Soluble proteins were extracted from SOX symbiont enrichments, host cytosolic fractions, and whole gill and foot tissue in biological duplicates (Appendix 4). All proteome samples were obtained from the M82-3 cruise. Membrane proteins were extracted from the SOX symbiont enrichments and the
whole gill tissue samples. We used 1D-PAGE followed by liquid chromatography (1D-PAGE-LC) to separate proteins and peptides as described previously with minor modifications (Heinz et al., 2012). MS/MS spectra were acquired with a LTQ Orbitrap Velos mass spectrometer (Thermo Fisher, Bremen, Germany) for soluble proteins and a LTQ Orbitrap Classic (Thermo Fisher Scientific Inc., Waltham, MA) for membrane proteins (Appendix 4). MS/MS data were searched against two databases using the SEQUEST algorithm (Eng et al., 1994). The first database, designated 'reduced' database, contained protein sequences from the SOX and MOX symbionts from *Bathymodiolus*, as well as from the host (see Appendix 4). The second database contained in addition sequences from host-related bivalves and symbiont-related bacteria (*Supplementary file 1H*). False discovery rates were determined with the target-decoy search strategy as described by Elias and Gygi (Elias and Gygi, 2007; Kleiner et al., 2012a).

## Acknowledgements

We thank the captains and crews of the research vessels and ROVs involved in the sampling effort, including Francois Lallier and Arnaud Tanguy who provided samples for sequencing from the MoMARETO cruise (Sarrazin et al., 2006) with the N/O Pourquoi Pas? and ROV Victor 6000. This project was funded by the Max Planck Society, an ERC Advanced Grant and a Gordon and Betty Moore Foundation Marine Microbial Initiative Investigator Award to ND, the DFG Cluster of Excellence 'The Ocean in the Earth System' at MARUM, Bremen, a DAAD scholarship to LS, OIST funding to NS, and the Marie Curie Initial Training Network 'Symbiomics' (project no. 264774) for the proteomics. We thank Harald Gruber-Vodicka, Hanno Teeling, and Brandon Seah for useful discussions, Michael Richter for his help with GenDB, and Sixing Huang for providing a Python script. We thank the three reviewers whose excellent suggestions helped us to improve the manuscript.

## Additional information

### Funding

<table><thead><tr><td>Funder</td><td>Grant reference</td><td>Author</td></tr></thead><tbody><tr><td>Max-Planck- Gesellschaft (MPG)</td><td></td><td>Lizbeth Sayavedra, Manuel Kleiner, Silke Wetzel, Dennis Fink, Nicole Dubilier, Jillian M Petersen</td></tr><tr><td>Gordon and Betty Moore Foundation</td><td>MMI Investigator Grant 3811</td><td>Nicole Dubilier</td></tr><tr><td>European Research Council (ERC)</td><td>Bathybiome</td><td>Lizbeth Sayavedra, Manuel Kleiner, Silke Wetzel, Dennis Fink, Nicole Dubilier, Jillian M Petersen</td></tr><tr><td>Deutsche Forschungsgemeinschaft (DFG)</td><td></td><td>Lizbeth Sayavedra, Manuel Kleiner, Silke Wetzel, Nicole Dubilier, Jillian M Petersen</td></tr><tr><td>German Academic Exchange Service (DAAD)</td><td></td><td>Lizbeth Sayavedra</td></tr><tr><td>European Research Council (ERC)</td><td>Symbiomics ITN 264774</td><td>Ruby Ponnudurai, Thomas Schweder, Stephanie Markert</td></tr><tr><td>Okinawa Institute of Science and Technology Graduate University (OIST)</td><td></td><td>Nori Satoh</td></tr></tbody></table>

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

### Author contributions

LS, Wrote the paper; analyzed genomic and transcriptomic data, Conception and design, Acquisition of data, Analysis and interpretation of data, Drafting or revising the article; MK, RP, Performed proteomic analyses and processed mass spectrometry data, Drafting or revising the article; SW, Developed symbiont extraction method, Drafting or revising the article; EP, VB, Sequenced and assembled samples sent to Genoscope, Drafting or revising the article; NS, ES, Sequenced and
whole gill tissue samples. We used 1D-PAGE followed by liquid chromatography (1D-PAGE-LC) to separate proteins and peptides as described previously with minor modifications (Heinz et al., 2012). MS/MS spectra were acquired with a LTQ Orbitrap Velos mass spectrometer (Thermo Fisher, Bremen, Germany) for soluble proteins and a LTQ Orbitrap Classic (Thermo Fisher Scientific Inc., Waltham, MA) for membrane proteins (Appendix 4). MS/MS data were searched against two databases using the SEQUEST algorithm (Eng et al., 1994). The first database, designated 'reduced' database, contained protein sequences from the SOX and MOX symbionts from *Bathymodiolus*, as well as from the host (see Appendix 4). The second database contained in addition sequences from host-related bivalves and symbiont-related bacteria (*Supplementary file 1H*). False discovery rates were determined with the target-decoy search strategy as described by Elias and Gygi (Elias and Gygi, 2007; Kleiner et al., 2012a).

## Acknowledgements

We thank the captains and crews of the research vessels and ROVs involved in the sampling effort, including Francois Lallier and Arnaud Tanguy who provided samples for sequencing from the MoMARETO cruise (Sarrazin et al., 2006) with the N/O Pourquoi Pas? and ROV Victor 6000. This project was funded by the Max Planck Society, an ERC Advanced Grant and a Gordon and Betty Moore Foundation Marine Microbial Initiative Investigator Award to ND, the DFG Cluster of Excellence 'The Ocean in the Earth System' at MARUM, Bremen, a DAAD scholarship to LS, OIST funding to NS, and the Marie Curie Initial Training Network 'Symbiomics' (project no. 264774) for the proteomics. We thank Harald Gruber-Vodicka, Hanno Teeling, and Brandon Seah for useful discussions, Michael Richter for his help with GenDB, and Sixing Huang for providing a Python script. We thank the three reviewers whose excellent suggestions helped us to improve the manuscript.

## Additional information

### Funding

<table><thead><tr><td>Funder</td><td>Grant reference</td><td>Author</td></tr></thead><tbody><tr><td>Max-Planck- Gesellschaft (MPG)</td><td></td><td>Lizbeth Sayavedra, Manuel Kleiner, Silke Wetzel, Dennis Fink, Nicole Dubilier, Jillian M Petersen</td></tr><tr><td>Gordon and Betty Moore Foundation</td><td>MMI Investigator Grant 3811</td><td>Nicole Dubilier</td></tr><tr><td>European Research Council (ERC)</td><td>Bathybiome</td><td>Lizbeth Sayavedra, Manuel Kleiner, Silke Wetzel, Dennis Fink, Nicole Dubilier, Jillian M Petersen</td></tr><tr><td>Deutsche Forschungsgemeinschaft (DFG)</td><td></td><td>Lizbeth Sayavedra, Manuel Kleiner, Silke Wetzel, Nicole Dubilier, Jillian M Petersen</td></tr><tr><td>German Academic Exchange Service (DAAD)</td><td></td><td>Lizbeth Sayavedra</td></tr><tr><td>European Research Council (ERC)</td><td>Symbiomics ITN 264774</td><td>Ruby Ponnudurai, Thomas Schweder, Stephanie Markert</td></tr><tr><td>Okinawa Institute of Science and Technology Graduate University (OIST)</td><td></td><td>Nori Satoh</td></tr></tbody></table>

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

### Author contributions

LS, Wrote the paper; analyzed genomic and transcriptomic data, Conception and design, Acquisition of data, Analysis and interpretation of data, Drafting or revising the article; MK, RP, Performed proteomic analyses and processed mass spectrometry data, Drafting or revising the article; SW, Developed symbiont extraction method, Drafting or revising the article; EP, VB, Sequenced and assembled samples sent to Genoscope, Drafting or revising the article; NS, ES, Sequenced and
assembled samples in OIST, Drafting or revising the article; DF, DF performed gradient centrifugation during the M82-3 cruise, Drafting or revising the article; CB, TBHR, PR, MBS, Sequenced transcriptomes from B. sp., Drafting or revising the article; DB, TS, SM, Performed proteomic analyses and processed mass spectrometry data, Drafting or revising the article; ND, Conception and design, Drafting or revising the article; JMP, Wrote the paper, Conception and design, Drafting or revising the article

## Additional files

### Supplementary file

* Supplementary file 1. (A) Number of mobile elements in the genomes compared in this study. (B) Genomes with toxin-related genes (TRGs) similar to those of the sulfur-oxidizing (SOX) symbionts of *Bathymodiolus*. The number of genes per TRGs class is shown. (C) p-values obtained with one-way Permanova were corrected with Bonferroni correction for multiple testing. Number of TRGs per genome was normalized to the total gene count. (D) Transcriptome counts of three individuals from *B. azoricus* and three individuals from *B. sp.* were mapped to their respective reference genomes with Rockhopper. Expression values of TRGs were normalized to the expression of RubisCO. (E) Samples used in this study. (F) Primer sequences and annealing temperatures used to detect genome rearrangements. (G) Metagenomes and metatranscriptomes enriched in SUP05 from oxygen minimum zones (OMZ) or hydrothermal vents. (H) Amino acid sequences from the following genomes were used in the reference database for proteomic analysis (IncDB). The genomes belong to relatives of the SOX and methane-oxidizing (MOX) symbionts of *B. azoricus*, as well as the mussel host. (I) Details of expressed toxin-related proteins identified with proteomics. The values are given in % normalized spectral abundance factor (NSAF), which is a normalized spectral abundance factor that gives the relative abundance of a protein in a sample in %.

DOI: 10.7554/eLife.07966.023

### Major datasets

The following datasets were generated:

<table><thead><tr><td>Author(s)</td><td>Year</td><td>Dataset title</td><td>Dataset ID and/or URL</td><td>Database, license, and accessibility information</td></tr></thead><tbody><tr><td>Sayavedra L, Kleiner M, Ponnudurai R, Wetzel S, Pelletier E, Barbe V, Satoh N, Shoguchi E, Fink D, Breusing C, Reusch TBH, Rosenstiel P, Schilhabel M, Becher D, Schweder T, Markert S, Dubilier N, Petersen JM</td><td>2015</td><td>Endosymbiont of Bathymodiolus sp.</td><td>http://www.ebi.ac.uk/ena/data/view/PRJNA65421</td><td>Publicly available at the EBI European Nucleotide Archive (Accession no: PRJNA65421).</td></tr><tr><td>Sayavedra L, Kleiner M, Ponnudurai R, Wetzel S, Pelletier E, Barbe V, Satoh N, Shoguchi E, Fink D, Breusing C, Reusch TBH, Rosenstiel P, Schilhabel M, Becher D, Schweder T, Markert S, Dubilier N, Petersen JM</td><td>2015</td><td>Genome of Sulfur-oxidizer endosymbiont of Bathymodiolus azoricus from Menez Gwen (BazSymA)</td><td>http://www.ebi.ac.uk/ena/data/view/PRJEB8263</td><td>Publicly available at the EBI European Nucleotide Archive (Accession no: PRJEB8263).</td></tr><tr><td>Sayavedra L, Kleiner M, Ponnudurai R, Wetzel S, Pelletier E, Barbe V, Satoh N, Shoguchi E, Fink D, Breusing C, Reusch TBH, Rosenstiel P, Schilhabel M, Becher D, Schweder T, Markert S, Dubilier N, Petersen JM</td><td>2015</td><td>Genome of Sulfur-oxidizer endosymbiont of Bathymodiolus azoricus from Menez Gwen (BazSymB)</td><td>http://www.ebi.ac.uk/ena/data/view/PRJEB8264</td><td>Publicly available at the EBI European Nucleotide Archive (Accession no: PRJEB8264).</td></tr></tbody></table>
assembled samples in OIST, Drafting or revising the article; DF, DF performed gradient centrifugation during the M82-3 cruise, Drafting or revising the article; CB, TBHR, PR, MBS, Sequenced transcriptomes from B. sp., Drafting or revising the article; DB, TS, SM, Performed proteomic analyses and processed mass spectrometry data, Drafting or revising the article; ND, Conception and design, Drafting or revising the article; JMP, Wrote the paper, Conception and design, Drafting or revising the article

## Additional files

### Supplementary file

* Supplementary file 1. (A) Number of mobile elements in the genomes compared in this study. (B) Genomes with toxin-related genes (TRGs) similar to those of the sulfur-oxidizing (SOX) symbionts of *Bathymodiolus*. The number of genes per TRGs class is shown. (C) p-values obtained with one-way Permanova were corrected with Bonferroni correction for multiple testing. Number of TRGs per genome was normalized to the total gene count. (D) Transcriptome counts of three individuals from *B. azoricus* and three individuals from *B. sp.* were mapped to their respective reference genomes with Rockhopper. Expression values of TRGs were normalized to the expression of RubisCO. (E) Samples used in this study. (F) Primer sequences and annealing temperatures used to detect genome rearrangements. (G) Metagenomes and metatranscriptomes enriched in SUP05 from oxygen minimum zones (OMZ) or hydrothermal vents. (H) Amino acid sequences from the following genomes were used in the reference database for proteomic analysis (IncDB). The genomes belong to relatives of the SOX and methane-oxidizing (MOX) symbionts of *B. azoricus*, as well as the mussel host. (I) Details of expressed toxin-related proteins identified with proteomics. The values are given in % normalized spectral abundance factor (NSAF), which is a normalized spectral abundance factor that gives the relative abundance of a protein in a sample in %.

DOI: 10.7554/eLife.07966.023

### Major datasets

The following datasets were generated:

<table><thead><tr><td>Author(s)</td><td>Year</td><td>Dataset title</td><td>Dataset ID and/or URL</td><td>Database, license, and accessibility information</td></tr></thead><tbody><tr><td>Sayavedra L, Kleiner M, Ponnudurai R, Wetzel S, Pelletier E, Barbe V, Satoh N, Shoguchi E, Fink D, Breusing C, Reusch TBH, Rosenstiel P, Schilhabel M, Becher D, Schweder T, Markert S, Dubilier N, Petersen JM</td><td>2015</td><td>Endosymbiont of Bathymodiolus sp.</td><td>http://www.ebi.ac.uk/ena/data/view/PRJNA65421</td><td>Publicly available at the EBI European Nucleotide Archive (Accession no: PRJNA65421).</td></tr><tr><td>Sayavedra L, Kleiner M, Ponnudurai R, Wetzel S, Pelletier E, Barbe V, Satoh N, Shoguchi E, Fink D, Breusing C, Reusch TBH, Rosenstiel P, Schilhabel M, Becher D, Schweder T, Markert S, Dubilier N, Petersen JM</td><td>2015</td><td>Genome of Sulfur-oxidizer endosymbiont of Bathymodiolus azoricus from Menez Gwen (BazSymA)</td><td>http://www.ebi.ac.uk/ena/data/view/PRJEB8263</td><td>Publicly available at the EBI European Nucleotide Archive (Accession no: PRJEB8263).</td></tr><tr><td>Sayavedra L, Kleiner M, Ponnudurai R, Wetzel S, Pelletier E, Barbe V, Satoh N, Shoguchi E, Fink D, Breusing C, Reusch TBH, Rosenstiel P, Schilhabel M, Becher D, Schweder T, Markert S, Dubilier N, Petersen JM</td><td>2015</td><td>Genome of Sulfur-oxidizer endosymbiont of Bathymodiolus azoricus from Menez Gwen (BazSymB)</td><td>http://www.ebi.ac.uk/ena/data/view/PRJEB8264</td><td>Publicly available at the EBI European Nucleotide Archive (Accession no: PRJEB8264).</td></tr></tbody></table>
<table><thead><tr><td>Author(s)</td><td>Year</td><td>Dataset title</td><td>Dataset ID and/or URL</td><td>Database, license, and accessibility information</td></tr></thead><tbody><tr><td>Sayavedra L, Kleiner M, Ponnudurai R, Wetzel S, Pelletier E, Barbe V, Satoh N, Shoguchi E, Fink D, Breusing C, Reusch TBH, Rosenstiel P, Schilhabel M, Becher D, Schweder T, Markert S, Dubilier N, Petersen JM</td><td>2015</td><td>Mitochondrial COI gene for cytochrome oxidase subunit 1</td><td>http://www.ebi.ac.uk/ena/data/view/LN833433-LN833440</td><td>Publicly available at the EBI European Nucleotide Archive (Accession no: LN833433-LN833440).</td></tr><tr><td>Sayavedra L, Kleiner M, Ponnudurai R, Wetzel S, Pelletier E, Barbe V, Satoh N, Shoguchi E, Fink D, Breusing C, Reusch TBH, Rosenstiel P, Schilhabel M, Becher D, Schweder T, Markert S, Dubilier N, Petersen JM</td><td>2015</td><td>Transcriptome of Sulfur-oxidizer endosymbiont of Bathymodiolus azoricus</td><td>http://www.ebi.ac.uk/ena/data/view/PRJEB7941</td><td>Publicly available at the EBI European Nucleotide Archive (Accession no: PRJEB7941).</td></tr><tr><td>Sayavedra L, Kleiner M, Ponnudurai R, Wetzel S, Pelletier E, Barbe V, Satoh N, Shoguchi E, Fink D, Breusing C, Reusch TBH, Rosenstiel P, Schilhabel M, Becher D, Schweder T, Markert S, Dubilier N, Petersen JM</td><td>2015</td><td>Transcriptome of Sulfur-oxidizer endosymbiont of Bathymodiolus sp. 9 South</td><td>http://www.ebi.ac.uk/ena/data/view/PRJEB7943</td><td>Publicly available at the EBI European Nucleotide Archive (Accession no: PRJEB7943).</td></tr><tr><td>Sayavedra L, Kleiner M, Ponnudurai R, Wetzel S, Pelletier E, Barbe V, Satoh N, Shoguchi E, Fink D, Breusing C, Reusch TBH, Rosenstiel P, Schilhabel M, Becher D, Schweder T, Markert S, Dubilier N, Petersen JM</td><td>2015</td><td>Bathymodiolus azoricus thioautotrophic gill symbiont partial 16S rRNA gene, isolate BazSymB</td><td>http://www.ebi.ac.uk/ena/data/view/LN871183</td><td>Publicly available at the EBI European Nucleotide Archive (Accession no: LN871183).</td></tr></tbody></table>

## References

* Achaz G, Coissac E, Netter P, Rocha EPC. 2003. Associations between inverted repeats and the structural evolution of bacterial genomes. *Genetics* **164**:1279-1289.

* Altschul SF, Madden TL, Schäffer AA, Zhang J, Zhang Z, Miller W, Lipman DJ. 1997. Gapped BLAST and PSI-BLAST: a new generation of protein database search programs. *Nucleic Acids Research* **25**:3389-3402. doi: 10.1093/nar/25.17.3389.

* Anantharaman K, Breier JA, Sheik CS, Dick GJ. 2012. Evidence for hydrogen oxidation and metabolic plasticity in widespread deep-sea sulfur-oxidizing bacteria. *Proceedings of the National Academy of Sciences of USA* **110**:330-335. doi: 10.1073/pnas.1215340110.

* Backert S, Meyer TF. 2006. Type IV secretion systems and their effectors in bacterial pathogenesis. *Current Opinion in Microbiology* **9**:207-217. doi: 10.1016/j.mib.2006.02.008.

* Badger MR, Bek EJ. 2008. Multiple RuBisCO forms in proteobacteria: their functional significance in relation to CO₂ acquisition by the CBB cycle. *Journal of Experimental Botany* **59**:1525-1541. doi: 10.1093/jxb/erm297.

* Bailly X, Olivieri I, De Mita S, Cleyet-Marel J-C, Béna G. 2006. Recombination and selection shape the molecular diversity pattern of nitrogen-fixing *Sinorhizobium* sp. associated to *Medicago*. *Molecular Ecology* **15**:2719-2734. doi: 10.1111/j.1365-294X.2006.02969.x.

* Bateman A, Coin L, Durbin R, Finn RD, Hollich V, Griffiths-Jones S, Khanna A, Marshall M, Moxon S, Sonnhammer EL, Studholme DJ, Yeats C, Eddy SR. 2004. The Pfam protein families database. *Nucleic Acids Research* **32**:D138-D141. doi: 10.1093/nar/gkh121.

* Benz J, Meinhart A. 2014. Antibacterial effector/immunity systems: it's just the tip of the iceberg. *Current Opinion in Microbiology* **17**:1-10. doi: 10.1016/j.mib.2013.11.002.
<table><thead><tr><td>Author(s)</td><td>Year</td><td>Dataset title</td><td>Dataset ID and/or URL</td><td>Database, license, and accessibility information</td></tr></thead><tbody><tr><td>Sayavedra L, Kleiner M, Ponnudurai R, Wetzel S, Pelletier E, Barbe V, Satoh N, Shoguchi E, Fink D, Breusing C, Reusch TBH, Rosenstiel P, Schilhabel M, Becher D, Schweder T, Markert S, Dubilier N, Petersen JM</td><td>2015</td><td>Mitochondrial COI gene for cytochrome oxidase subunit 1</td><td>http://www.ebi.ac.uk/ena/data/view/LN833433-LN833440</td><td>Publicly available at the EBI European Nucleotide Archive (Accession no: LN833433-LN833440).</td></tr><tr><td>Sayavedra L, Kleiner M, Ponnudurai R, Wetzel S, Pelletier E, Barbe V, Satoh N, Shoguchi E, Fink D, Breusing C, Reusch TBH, Rosenstiel P, Schilhabel M, Becher D, Schweder T, Markert S, Dubilier N, Petersen JM</td><td>2015</td><td>Transcriptome of Sulfur-oxidizer endosymbiont of Bathymodiolus azoricus</td><td>http://www.ebi.ac.uk/ena/data/view/PRJEB7941</td><td>Publicly available at the EBI European Nucleotide Archive (Accession no: PRJEB7941).</td></tr><tr><td>Sayavedra L, Kleiner M, Ponnudurai R, Wetzel S, Pelletier E, Barbe V, Satoh N, Shoguchi E, Fink D, Breusing C, Reusch TBH, Rosenstiel P, Schilhabel M, Becher D, Schweder T, Markert S, Dubilier N, Petersen JM</td><td>2015</td><td>Transcriptome of Sulfur-oxidizer endosymbiont of Bathymodiolus sp. 9 South</td><td>http://www.ebi.ac.uk/ena/data/view/PRJEB7943</td><td>Publicly available at the EBI European Nucleotide Archive (Accession no: PRJEB7943).</td></tr><tr><td>Sayavedra L, Kleiner M, Ponnudurai R, Wetzel S, Pelletier E, Barbe V, Satoh N, Shoguchi E, Fink D, Breusing C, Reusch TBH, Rosenstiel P, Schilhabel M, Becher D, Schweder T, Markert S, Dubilier N, Petersen JM</td><td>2015</td><td>Bathymodiolus azoricus thioautotrophic gill symbiont partial 16S rRNA gene, isolate BazSymB</td><td>http://www.ebi.ac.uk/ena/data/view/LN871183</td><td>Publicly available at the EBI European Nucleotide Archive (Accession no: LN871183).</td></tr></tbody></table>

## References

* Achaz G, Coissac E, Netter P, Rocha EPC. 2003. Associations between inverted repeats and the structural evolution of bacterial genomes. *Genetics* **164**:1279-1289.

* Altschul SF, Madden TL, Schäffer AA, Zhang J, Zhang Z, Miller W, Lipman DJ. 1997. Gapped BLAST and PSI-BLAST: a new generation of protein database search programs. *Nucleic Acids Research* **25**:3389-3402. doi: 10.1093/nar/25.17.3389.

* Anantharaman K, Breier JA, Sheik CS, Dick GJ. 2012. Evidence for hydrogen oxidation and metabolic plasticity in widespread deep-sea sulfur-oxidizing bacteria. *Proceedings of the National Academy of Sciences of USA* **110**:330-335. doi: 10.1073/pnas.1215340110.

* Backert S, Meyer TF. 2006. Type IV secretion systems and their effectors in bacterial pathogenesis. *Current Opinion in Microbiology* **9**:207-217. doi: 10.1016/j.mib.2006.02.008.

* Badger MR, Bek EJ. 2008. Multiple RuBisCO forms in proteobacteria: their functional significance in relation to CO₂ acquisition by the CBB cycle. *Journal of Experimental Botany* **59**:1525-1541. doi: 10.1093/jxb/erm297.

* Bailly X, Olivieri I, De Mita S, Cleyet-Marel J-C, Béna G. 2006. Recombination and selection shape the molecular diversity pattern of nitrogen-fixing *Sinorhizobium* sp. associated to *Medicago*. *Molecular Ecology* **15**:2719-2734. doi: 10.1111/j.1365-294X.2006.02969.x.

* Bateman A, Coin L, Durbin R, Finn RD, Hollich V, Griffiths-Jones S, Khanna A, Marshall M, Moxon S, Sonnhammer EL, Studholme DJ, Yeats C, Eddy SR. 2004. The Pfam protein families database. *Nucleic Acids Research* **32**:D138-D141. doi: 10.1093/nar/gkh121.

* Benz J, Meinhart A. 2014. Antibacterial effector/immunity systems: it's just the tip of the iceberg. *Current Opinion in Microbiology* **17**:1-10. doi: 10.1016/j.mib.2013.11.002.
Bettencourt R, Roch P, Stefanni S, Rosa D, Colaço A, Serrão Santos R. 2007. Deep sea immunity: unveiling immune constituents from the hydrothermal vent mussel *Bathymodiolus azoricus*. *Marine Environmental Research* **64**:108–127. doi: 10.1016/j.marenvres.2006.12.010.

Bhavsar AP, Guttman JA, Finlay BB. 2007. Manipulation of host-cell pathways by bacterial pathogens. *Nature* **449**:827–834. doi: 10.1038/nature06247.

Boeckmann B. 2003. The SWISS-PROT protein knowledgebase and its supplement TrEMBL in 2003. *Nucleic Acids Research* **31**:365–370. doi: 10.1093/nar/gkg095.

Bradford MM. 1976. A rapid and sensitive method for the quantitation of microgram quantities of protein utilizing the principle of protein-dye binding. *Analytical Biochemistry* **72**:248–254.

Brinkhoff T, Giebel H-A, Simon M. 2008. Diversity, ecology, and genomics of the Roseobacter clade: a short overview. *Archives of Microbiology* **189**:531–539. doi: 10.1007/s00203-008-0353-y.

Buchan A, González JM, Moran MA. 2005. Overview of the marine roseobacter lineage. *Applied and Environment Microbiology* **71**:5665–5677. doi: 10.1128/AEM.71.10.5665-5677.2005.

Caspi R, Altman T, Billington R, Dreher K, Foerster H, Fulcher CA, Holland TA, Keseler IM, Kothari A, Kubo A, Krummenacker M, Latendresse M, Mueller LA, Ong Q, Paley S, Subhraveti P, Weaver DS, Weerasinghe D, Zhang P, Karp PD. 2014. The MetaCyc database of metabolic pathways and enzymes and the BioCyc collection of pathway/genome databases. *Nucleic Acids Research* **42**:D459–D471. doi: 10.1093/nar/gkt1103.

Cavanaugh CM, Robinson JJ. 1996. CO₂ fixation in chemoautotroph-invertebrate symbioses: expression of form I and form II RuBisCO. In: Lidstrom ME, Tabita FR, editors. *Microbial growth on C1 compounds*: Springer Netherlands. p. 285–292.

Cavanaugh CM, McKiness ZP, Newton ILG, Stewart FJ. 2006. Marine chemosynthetic symbioses. In: Dworkin M, Falkow S, Rosenberg E, Schleifer K-H, Stackebrandt E, editors. *The Prokaryotes*: Springer New York. p. 475–507.

Chambers MC, Maclean B, Burke R, Amodei D, Ruderman DL, Neumann S, Gatto L, Fischer B, Pratt B, Egertson J, Hoff K, Kessner D, Tasman N, Shulman N, Frewen B, Baker TA, Brusniak MY, Paulse C, Creasy D, Flashner L, Kani K, Moulding C, Seymour SL, Nuwaysir LM, Lefebvre B, Kuhlmann F, Roark J, Rainer P, Detlev S, Hemenway T, Huhmer A, Langridge J, Connolly B, Chadick T, Holly K, Eckels J, Deutsch EW, Moritz RL, Katz JE, Agus DB, MacCoss M, Tabb DL, Mallick P. 2012. A cross-platform toolkit for mass spectrometry and proteomics. *Nature Biotechnology* **30**:918–920. doi: 10.1038/nbt.2377.

Christie-Oleza JA, Piña-Villalonga JM, Bosch R, Nogales B, Armengaud J. 2012. Comparative proteogenomics of twelve Roseobacter exoproteomes reveals different adaptive strategies among these marine bacteria. *Molecular & Cellular Proteomics* **11**:M111.013110. doi: 10.1074/mcp.M111.013110.

Clayton AL, Oakeson KF, Gutin M, Pontes A, Dunn DM, von Niederhausern AC, Weiss RB, Fisher M, Dale C. 2012. A novel human-infection-derived bacterium provides insights into the evolutionary origins of mutualistic insect-bacterial symbioses. *PLOS Genetics* **8**:e1002990. doi: 10.1371/journal.pgen.1002990.

Cossart P, Sansonetti PJ. 2004. Bacterial invasion: the paradigms of enteroinvasive pathogens. *Science* **304**:242–248. doi: 10.1126/science.1090124.

Coustau C, Renaud F, Maillard C, Pasteur N, Delay B. 1991. Differential susceptibility to a trematode parasite among genotypes of the *Mytilus edulis/galloprovincialis* complex. *Genetics Research* **57**:207–212. doi: 10.1017/S0016672300029359.

da Silva PM, Magalhães AR, Barracco MA. 2002. Effects of *Bucephalus* sp. (Trematoda: Bucephalidae) on *Perna perna* mussels from a culture station in Ratones Grande Island, Brazil. *Journal of Invertebrate Pathology* **79**:154–162. doi: 10.1016/S0022-2011(02)00026-5.

Dale C, Moran NA. 2006. Molecular interactions between bacterial symbionts and their hosts. *Cell* **126**:453–465. doi: 10.1016/j.cell.2006.07.014.

Dale C, Young SA, Haydon DT, Welburn SC. 2001. The insect endosymbiont *Sodalis glossinidius* utilizes a type III secretion system for cell invasion. *Proceedings of the National Academy of Sciences of USA* **98**:1883–1888. doi: 10.1073/pnas.021450998.

Davies RL, Campbell S, Whittam TS. 2002. Mosaic structure and molecular evolution of the leukotoxin operon (IktCABD) in Mannheimia (Pasteurella) haemolytica, *Mannheimia glucosida*, and *Pasteurella trehalosi*. *Journal of Bacteriology* **184**:266–277. doi: 10.1128/JB.184.1.266-277.2002.

Davis JJ, Olsen GJ. 2010. Characterizing the native codon usages of a genome: an axis projection approach. *Molecular Biology and Evolution* **28**:211–221. doi: 10.1093/molbev/msq185.

De Wit P, Pespeni MH, Ladner JT, Barshis DJ, Seneca F, Jaris H, Therkildsen NO, Morikawa M, Palumbi SR. 2012. The simple fool's guide to population genomics via RNA-Seq: an introduction to high-throughput sequencing data analysis. *Molecular Ecology Resources* **12**:1058–1067. doi: 10.1111/1755-0998.12003.

DeChaine EG, Bates AE, Shank TM, Cavanaugh CM. 2006. Off-axis symbiosis found: characterization and biogeography of bacterial symbionts of *Bathymodiolus* mussels from Lost City hydrothermal vents. *Environmental Microbiology* **8**:1902–1912. doi: 10.1111/j.1462-2920.2005.01113.x.

Degnan PH, Moran NA. 2008. Diverse phage-encoded toxins in a protective insect endosymbiont. *Applied and Environmental Microbiology* **74**:6782–6791. doi: 10.1128/AEM.01285-08.

Delcher AL, Bratke KA, Powers EC, Salzberg SL. 2007. Identifying bacterial genes and endosymbiont DNA with Glimmer. *Bioinformatics* **23**:673–679. doi: 10.1093/bioinformatics/btm009.

Distel DL, Lee HK, Cavanaugh CM. 1995. Intracellular coexistence of methano- and thioautotrophic bacteria in a hydrothermal vent mussel. *Proceedings of the National Academy of Sciences of USA* **92**:9598–9602. doi: 10.1073/pnas.92.21.9598.

Dmytrenko O, Russell SL, Loo WT, Fontanez KM, Liao L, Roeselers G, Sharma R, Stewart FJ, Newton IL, Woyke T, Wu D, Lang JM, Eisen JA, Cavanaugh CM. 2014. The genome of the intracellular bacterium of the
Bettencourt R, Roch P, Stefanni S, Rosa D, Colaço A, Serrão Santos R. 2007. Deep sea immunity: unveiling immune constituents from the hydrothermal vent mussel *Bathymodiolus azoricus*. *Marine Environmental Research* **64**:108-127. doi: 10.1016/j.marenvres.2006.12.010.

Bhavsar AP, Guttman JA, Finlay BB. 2007. Manipulation of host-cell pathways by bacterial pathogens. *Nature* **449**:827-834. doi: 10.1038/nature06247.

Boeckmann B. 2003. The SWISS-PROT protein knowledgebase and its supplement TrEMBL in 2003. *Nucleic Acids Research* **31**:365–370. doi: 10.1093/nar/gkg095.

Bradford MM. 1976. A rapid and sensitive method for the quantitation of microgram quantities of protein utilizing the principle of protein-dye binding. *Analytical Biochemistry* **72**:248-254.

Brinkhoff T, Giebel H-A, Simon M. 2008. Diversity, ecology, and genomics of the Roseobacter clade: a short overview. *Archives of Microbiology* **189**:531–539. doi: 10.1007/s00203-008-0353-y.

Buchan A, González JM, Moran MA. 2005. Overview of the marine roseobacter lineage. *Applied and Environment Microbiology* **71**:5665-5677. doi: 10.1128/AEM.71.10.5665-5677.2005.

Caspi R, Altman T, Billington R, Dreher K, Foerster H, Fulcher CA, Holland TA, Keseler IM, Kothari A, Kubo A, Krummenacker M, Latendresse M, Mueller LA, Ong Q, Paley S, Subhraveti P, Weaver DS, Weerasinghe D, Zhang P, Karp PD. 2014. The MetaCyc database of metabolic pathways and enzymes and the BioCyc collection of pathway/genome databases. *Nucleic Acids Research* **42**:D459–D471. doi: 10.1093/nar/gkt1103.

Cavanaugh CM, Robinson JJ. 1996. CO₂ fixation in chemoautotroph-invertebrate symbioses: expression of form I and form II RuBisCO. In: Lidstrom ME, Tabita FR, editors. *Microbial growth on C1 compounds*: Springer Netherlands. p. 285-292.

Cavanaugh CM, McKiness ZP, Newton ILG, Stewart FJ. 2006. Marine chemosynthetic symbioses. In: Dworkin M, Falkow S, Rosenberg E, Schleifer K-H, Stackebrandt E, editors. *The Prokaryotes*: Springer New York. p. 475-507.

Chambers MC, Maclean B, Burke R, Amodei D, Ruderman DL, Neumann S, Gatto L, Fischer B, Pratt B, Egertson J, Hoff K, Kessner D, Tasman N, Shulman N, Frewen B, Baker TA, Brusniak MY, Paulse C, Creasy D, Flashner L, Kani K, Moulding C, Seymour SL, Nuwaysir LM, Lefebvre B, Kuhlmann F, Roark J, Rainer P, Detlev S, Hemenway T, Huhmer A, Langridge J, Connolly B, Chadick T, Holly K, Eckels J, Deutsch EW, Moritz RL, Katz JE, Agus DB, MacCoss M, Tabb DL, Mallick P. 2012. A cross-platform toolkit for mass spectrometry and proteomics. *Nature Biotechnology* **30**:918–920. doi: 10.1038/nbt.2377.

Christie-Oleza JA, Piña-Villalonga JM, Bosch R, Nogales B, Armengaud J. 2012. Comparative proteogenomics of twelve Roseobacter exoproteomes reveals different adaptive strategies among these marine bacteria. *Molecular & Cellular Proteomics* **11**:M111.013110. doi: 10.1074/mcp.M111.013110.

Clayton AL, Oakeson KF, Gutin M, Pontes A, Dunn DM, von Niederhausern AC, Weiss RB, Fisher M, Dale C. 2012. A novel human-infection-derived bacterium provides insights into the evolutionary origins of mutualistic insect-bacterial symbioses. *PLOS Genetics* **8**:e1002990. doi: 10.1371/journal.pgen.1002990.

Cossart P, Sansonetti PJ. 2004. Bacterial invasion: the paradigms of enteroinvasive pathogens. *Science* **304**:242–248. doi: 10.1126/science.1090124.

Coustau C, Renaud F, Maillard C, Pasteur N, Delay B. 1991. Differential susceptibility to a trematode parasite among genotypes of the *Mytilus edulis/galloprovincialis* complex. *Genetics Research* **57**:207-212. doi: 10.1017/S0016672300029359.

da Silva PM, Magalhães AR, Barracco MA. 2002. Effects of *Bucephalus* sp. (Trematoda: Bucephalidae) on *Perna perna* mussels from a culture station in Ratones Grande Island, Brazil. *Journal of Invertebrate Pathology* **79**:154-162. doi: 10.1016/S0022-2011(02)00026-5.

Dale C, Moran NA. 2006. Molecular interactions between bacterial symbionts and their hosts. *Cell* **126**:453-465. doi: 10.1016/j.cell.2006.07.014.

Dale C, Young SA, Haydon DT, Welburn SC. 2001. The insect endosymbiont *Sodalis glossinidius* utilizes a type III secretion system for cell invasion. *Proceedings of the National Academy of Sciences of USA* **98**:1883-1888. doi: 10.1073/pnas.021450998.

Davies RL, Campbell S, Whittam TS. 2002. Mosaic structure and molecular evolution of the leukotoxin operon (IktCABD) in Mannheimia (Pasteurella) haemolytica, *Mannheimia glucosida*, and *Pasteurella trehalosi*. *Journal of Bacteriology* **184**:266–277. doi: 10.1128/JB.184.1.266-277.2002.

Davis JJ, Olsen GJ. 2010. Characterizing the native codon usages of a genome: an axis projection approach. *Molecular Biology and Evolution* **28**:211–221. doi: 10.1093/molbev/msq185.

De Wit P, Pespeni MH, Ladner JT, Barshis DJ, Seneca F, Jaris H, Therkildsen NO, Morikawa M, Palumbi SR. 2012. The simple fool's guide to population genomics via RNA-Seq: an introduction to high-throughput sequencing data analysis. *Molecular Ecology Resources* **12**:1058-1067. doi: 10.1111/1755-0998.12003.

DeChaine EG, Bates AE, Shank TM, Cavanaugh CM. 2006. Off-axis symbiosis found: characterization and biogeography of bacterial symbionts of *Bathymodiolus* mussels from Lost City hydrothermal vents. *Environmental Microbiology* **8**:1902–1912. doi: 10.1111/j.1462-2920.2005.01113.x.

Degnan PH, Moran NA. 2008. Diverse phage-encoded toxins in a protective insect endosymbiont. *Applied and Environmental Microbiology* **74**:6782–6791. doi: 10.1128/AEM.01285-08.

Delcher AL, Bratke KA, Powers EC, Salzberg SL. 2007. Identifying bacterial genes and endosymbiont DNA with Glimmer. *Bioinformatics* **23**:673–679. doi: 10.1093/bioinformatics/btm009.

Distel DL, Lee HK, Cavanaugh CM. 1995. Intracellular coexistence of methano- and thioautotrophic bacteria in a hydrothermal vent mussel. *Proceedings of the National Academy of Sciences of USA* **92**:9598-9602. doi: 10.1073/pnas.92.21.9598.

Dmytrenko O, Russell SL, Loo WT, Fontanez KM, Liao L, Roeselers G, Sharma R, Stewart FJ, Newton IL, Woyke T, Wu D, Lang JM, Eisen JA, Cavanaugh CM. 2014. The genome of the intracellular bacterium of the
coastal bivalve, *Solemya velum*: a blueprint for thriving in and out of symbiosis. *BMC Genomics* **15**:924. doi: 10.1186/1471-2164-15-924.

Dobrindt U, Hochhut B, Hentschel U, Hacker J. 2004. Genomic islands in pathogenic and environmental microorganisms. *Nature Reviews Microbiology* **2**:414–424. doi: 10.1038/nrmicro884.

Douglas AE. 2014. *The symbiotic habit*.

Dubilier N, Bergin C, Lott C. 2008. Symbiotic diversity in marine animals: the art of harnessing chemosynthesis. *Nature Reviews Microbiology* **6**:725–740. doi: 10.1038/nrmicro1992.

Duperron S, Bergin C, Zielinski F, Blazejak A, Perthaler A, McKiness ZP, DeChaine E, Cavanaugh CM, Dubilier N. 2006. A dual symbiosis shared by two mussel species, *Bathymodiolus azoricus* and *Bathymodiolus puteoserpentis* (Bivalvia: Mytilidae), from hydrothermal vents along the northern Mid-Atlantic Ridge. *Environmental Microbiology* **8**:1441–1447.

Duperron S, Halary S, Lorion J, Sibuet M, Gaill F. 2008. Unexpected co-occurrence of six bacterial symbionts in the gills of the cold seep mussel *Idas* sp. (Bivalvia: Mytilidae). *Environmental Microbiology* **10**:433–445. doi: 10.1111/j.1462-2920.2007.01465.x.

Duperron S, Gaudron SM, Rodrigues CF, Cunha MR, Decker C, Olu K. 2013. An overview of chemosynthetic symbioses in bivalves from the North Atlantic and Mediterranean sea. *Biogeosciences* **10**:3241–3267. doi: 10.5194/bg-10-3241-2013.

Elias JE, Gygi SP. 2007. Target-decoy search strategy for increased confidence in large-scale protein identifications by mass spectrometry. *Nature Methods* **4**:207–214. doi: 10.1038/nmeth1019.

Emanuelsson O, Brunak S, von Heijne G, Nielsen H. 2007. Locating proteins in the cell using TargetP, SignalP and related tools. *Nature Protocols* **2**:953–971. doi: 10.1038/nprot.2007.131.

Eng JK, McCormack AL, Yates JR III. 1994. An approach to correlate tandem mass spectral data of peptides with amino acid sequences in a protein database. *Journal of the American Society for Mass Spectrometry* **5**:976–989. doi: 10.1016/1044-0305(94)80016-2.

Eymann C, Dreisbach A, Albrecht D, Bernhardt J, Becher D, Gentner S, Tam le T, Büttner K, Buurman G, Scharf C, Venz S, Völker U, Hecker M. 2004. A comprehensive proteome map of growing *Bacillus subtilis* cells. *Proteomics* **4**:2849–2876. doi: 10.1002/pmic.200400907.

Finn RD, Clements J, Eddy SR. 2011. HMMER web server: interactive sequence similarity searching. *Nucleic Acids Research* **39**:W29. doi: 10.1093/nar/gkr367.

Fisher CR, Brooks JM, Vodenichar JS, Zande JM, Childress JJ, Burke RA Jr. 1993. The co-occurrence of methanotrophic and chemoautotrophic sulfur-oxidizing bacterial symbionts in a deep-sea mussel. *Marine Ecology* **14**:277–289. doi: 10.1111/j.1439-0485.1993.tb00001.x.

Florens L, Carozza MJ, Swanson SK, Fournier M, Coleman MK, Workman JL, Washburn MP. 2006. Analyzing chromatin remodeling complexes using shotgun proteomics and normalized spectral abundance factors. *Methods* **40**:303–311. doi: 10.1016/j.ymeth.2006.07.028.

Fontanez KM, Cavanaugh CM. 2014. Evidence for horizontal transmission from multilocus phylogeny of deep-sea mussel (Mytilidae) symbionts: horizontal transmission of mussel symbionts. *Environmental Microbiology* **16**:3608–3621. doi: 10.1111/1462-2920.12379.

Galagan JE. 2014. Genomic insights into tuberculosis. *Nature Reviews Genetics* **15**:307–320. doi: 10.1038/nrg3664.

Gardebrecht A, Markert S, Sievert SM, Felbeck H, Thürmer A, Albrecht D, Wollherr A, Kabisch J, Le Bris N, Lehmann R, Daniel R, Liesegang H, Hecker M, Schweder T. 2012. Physiological homogeneity among the endosymbionts of *Riftia pachyptila* and *Tevnia jerichonana* revealed by proteogenomics. *The ISME Journal* **6**:766–776. doi: 10.1038/ismej.2011.137.

Gebruk AV, Chevaldonné P, Shank T, Lutz RA, Vrijenhoek RC. 2000. Deep-sea hydrothermal vent communities of the Logatchev area (14'45'N, Mid-Atlantic Ridge): diverse biotopes and high biomass. *Journal of the Marine Biological Association of the UK* **80**:383–393. doi: 10.1017/S0025315499002088.

Goffredi SK, Yi H, Zhang Q, Klann JE, Struve IA, Vrijenhoek RC, Brown CT. 2014. Genomic versatility and functional variation between two dominant heterotrophic symbionts of deep-sea *Osedax* worms. *The ISME Journal* **8**:908–924. doi: 10.1038/ismej.2013.201.

Gogarten JP, Townsend JP. 2005. Horizontal gene transfer, genome innovation and evolution. *Nature Reviews Microbiology* **3**:679–687. doi: 10.1038/nrmicro1204.

Goodrich-Blair H, Clarke DJ. 2007. Mutualism and pathogenesis in *Xenorhabdus* and *Photorhabdus*: two roads to the same destination. *Molecular Microbiology* **64**:260–268. doi: 10.1111/j.1365-2958.2007.05671.x.

Haglund CM, Welch MD. 2011. Pathogens and polymers: microbe-host interactions illuminate the cytoskeleton. *The Journal of Cell Biology* **195**:7–17. doi: 10.1083/jcb.201103148.

Halary S, Riou V, Gaill F, Boudier T, Duperron S. 2008. 3D FISH for the quantification of methane- and sulphur-oxidizing endosymbionts in bacteriocytes of the hydrothermal vent mussel *Bathymodiolus azoricus*. *The ISME Journal* **2**:284–292. doi: 10.1038/ismej.2008.3.

Harada M, Yoshida T, Kuwahara H, Shimamura S, Takaki Y, Kato C, Miwa T, Miyake H, Maruyama T. 2009. Expression of genes for sulfur oxidation in the intracellular chemoautotrophic symbiont of the deep-sea bivalve *Calyptogena okutanii*. *Extremophiles* **13**:895–903. doi: 10.1007/s00792-009-0277-8.

Harmer TL, Rotjan RD, Nussbaumer AD, Bright M, Ng AW, DeChaine EG, Cavanaugh CM. 2008. Free-living tube worm endosymbionts found at deep-sea vents. *Applied and Environmental Microbiology* **74**:3895–3898. doi: 10.1128/AEM.02470-07.

Hart MM, Forsythe J, Oshowski B, Bücking H, Jansa J, Kiers ET. 2013. Hiding in a crowd—does diversity facilitate persistence of a low-quality fungal partner in the mycorrhizal symbiosis? *Symbiosis* **59**:47–56. doi: 10.1007/s13199-012-0197-8.