The bivalve Thyasira cf. gouldi hosts chemoautotrophic symbiont populations with strain level diversity

Sample collection

Thyasirids were collected from Bonne Bay, Newfoundland, Canada on nine occasions between October 2009 and May 2012 (Table S1); permits for field sampling (NL 572 11 and NL 992 12) were obtained from Fisheries and Oceans Canada. Sediment was collected using a Peterson grab (radius = 10.5 cm, length = 30 cm, volume = 0.01 m³) from three sites within the fjord (Fig. 1): Neddy’s Harbour (17–30 m depth), Deer Arm (29–36 m depth) and South East Arm (20–35 m depth). Thyasirids were retrieved from sediments using a sieve with 1 mm mesh, and symbiotic individuals (Thyasira cf. gouldi OTU 1 and 2, distinguished by their shell shape; Batstone et al., 2014), were retained.

DNA extraction and gene sequencing

The gills of symbiotic Thyasira cf. gouldi specimens were dissected and immediately frozen or stored in 95% ethanol. Following the protocol for animal tissues, total DNA was extracted from gills using QIAGEN DNeasy^® Blood and Tissue kit spin columns and stored at −20 °C in the elusion buffer provided. The PCR amplification of 16S rRNA and RuBisCO gene sequences was conducted using 12.5 µl of Green Dream Master Mix, 1.5 µl of template DNA, 1 µl of forward primer, 1 µl of reverse primer, and 9 µl of water. A 1323 bp fragment of the 16S rRNA gene was amplified using primers 27F (5′AGAGTTTGGATCMTGGCTCAG 3′) and 1492R (5′CGGTTACCTTGTTACGACTT 3′) (Lane, 1991). Thermocycler settings were: 94 °C for 3 min, 35 cycles of (94 °C for 1 min, 50 °C for 30 sec, 72 °C for 1.5 min) and a final extension of 72 °C for 10 min. A 296 bp fragment of the RuBisCO Form II gene (previously identified within the symbionts of T. cf. gouldi; Dufour et al., 2014) was amplified using primers 663F (5′ATCATCAARCTSGGCCTGCGTCCC 3′) and 1033R (5′MGAGGTGACSGCRCCGTGRCCRGCMCRTG 3′) (Widmer et al., 1998); initial denaturation was at 95 °C for 2 min, followed by 30 cycles of 95 °C for 1 min, 62 °C for 1 min and 72 °C for 30 sec, and a final elongation at 72 °C for 5 min. PCR products were cleaned using Agencourt AMPure XP Beads (Beckman Coulter, Brea, CA, USA) following the manufacturer’s protocol, and sent to The Center for Advanced Genomics, Toronto, Canada for Sanger sequencing.

Phylogenetic analysis

Sequences were checked for quality, manually trimmed from both ends, and corresponding forward and reverse sequences from a single clam individual were combined into contiguous sequences (contigs) using SEQUENCHER^® 5.1 (Gene Codes Corp., Ann Arbor, MI, USA).

Contigs were then aligned in MEGA 7 (Kumar, Stecher & Tamura, 2016) using the ClustalW algorithm (Thompson, Higgins & Gibson, 1994). We paid particular attention to any sites with double peaks in the chromatographs (i.e., where IUPAC degenerate base symbols were assigned by the sequencing software) on corresponding forward and reverse sequences: see Fig. S1 for examples. We examined the identity and position of double bases within the alignment, as they could signal the presence of more than one strain or phylotype within a particular host bivalve (as in the vesicomyid tubeworms examined by Vrijenhoek, Duhaime & Jones, 2007).

Maximum likelihood trees were constructed for both 16S rRNA and RuBisCO genes obtained herein, with the tree for the former gene including additional sequences from Batstone & Dufour (2016). Sequences with degenerate bases were not included in the RuBisCO tree as the high number of heterogeneous sequences considerably reduced bootstrap numbers. Appropriate models were identified using MEGA7. Three distinct clusters were identified in each tree, and representative sequences of the highest quality (and without ambiguities) were selected for comparison to sequences of chemosynthetic symbionts and free-living bacteria for which both 16S rRNA and RuBisCO Form II sequences were available in Genbank (Table S2).

The degree of divergence within and between phylotype groups was calculated in MEGA 7, using the Jukes-Cantor model for 16S rRNA and the Tamura-Nei model (Tamura & Nei, 1993) with 5 discrete gamma distributions for RuBisCO. Variance from the model was estimated with 1000 bootstrap replicates. Percent similarity between pairs of representative sequences from each phylotype was calculated using BLAST.

Statistical analysis of gene pairings and of spatial, temporal, and host size effects

We used chi-square tests to examine whether 16S rRNA and RuBisCO phylotypes identified unambiguously from a series of host specimens were independent of each other. A similar approach was used to assess independence between 16S and RuBisCO phylotypes and (1) sampling site, (2) sampling date, and (3) host size (set as the following categories: <3 mm, 3–4 mm, and >4 mm). Statistical analyses were performed using R (Ihaka & Gentleman, 1996).

Description of sequences and evolutionary patterns

From the bivalves examined in this study, we obtained 16S rRNA and RuBisCO Form II sequences from 54 and 50 host individuals, respectively (Table S1). A phylogenetic tree using all sequences from this study was constructed for each gene; only sequences with no degenerate bases were used for the RuBisCO tree (Figs. S2 and S3). The 16S rRNA sequences formed 3 phylotypes, in agreement with (Batstone & Dufour, 2016); these are hereafter referred to as 16S A, 16S B, and 16S C. RuBisCO sequences varied but could also be grouped into three phylotypes: RB 1, RB 2 and RB 3.

For each gene, evolutionary distances between and within phylotypes were calculated (Table 1). Although there was clear sequence variation within phylotypes (Table S3), no measurable evolutionary distances were observed within phylotypes of either gene, likely because most discrepancies were within wobble positions. Evolutionary distances between phylotypes were greater for the RuBisCO gene than for 16S rRNA. 16S A and 16S B are more similar to each other (evolutionary distance = 0.002) than to 16S C (evolutionary distances = 0.008 and 0.007, respectively). A similar pattern is seen among RuBisCO sequences, with RB 3 being the most distant, as we calculated distance values of 0.229 and 0.267 in comparison to RB 1 and RB 2, respectively. A clear divergence in the evolutionary history of 16S rRNA and RuBisCO genes within the Thyasira cf. gouldi symbionts was evident upon examination of phylogenetic trees that included sequences from the same free-living bacteria and chemosynthetic symbionts (Figs. 2 and 3).

Evidence for multiple symbiont phylotypes within a host

Sequences with degenerate bases were commonly observed (Table S3). Of the 54 16S rRNA sequences, 11 had a single instance of a high-quality call of multiple bases, while 11 others possessed between three and 12 degenerate bases in variable gene regions. Based on alignments, sequences with degenerate bases suggested a combination of slightly varying sequences within a single phylotype; in no cases did there appear to be more than one 16S rRNA phylotype within a single host bivalve.

Among the 50 RuBisCO sequences, 29 showed no ambiguities, while 18 showed high quality calls of multiple bases at sites varying within a given symbiont phylotype (i.e., the host appeared to contain multiple strains within a particular RB phylotype). Three individuals showed multiple calls at sites that suggested a mixture of symbiont phylotypes RB 2 and RB 3 (Fig. S3).

Phylotypes and associated host or symbiont characteristics

We obtained corresponding, unambiguously categorized 16S rRNA and RuBisCO sequences from 33 host specimens and noted patterns in gene pairings within host individuals (Table 2). Phylotypes 16S C and RB 1 appeared associated with each other: two of the four host specimens with 16S C had RB 1 (no RuBisCO phylotype data could be obtained from the remaining two host specimens due to insufficient volumes of extracted DNA). Similarly, all 15 hosts containing the 16S B phylotype possessed RB 2. In contrast, six of the 16 Thyasira cf. gouldi individuals with symbiont phylotype 16S A had RB 2, while 10 had RB 3. The chi-square test showed that 16S and RB phylotypes were not statistically independent of each other within host specimens (p-value = 0.00049).

Some patterns in the spatial distribution of symbiont 16S rRNA and RuBisCO phylotypes within Bonne Bay were observed, and chi-square tests revealed that site and phylotype were not statistically independent of each other (p-value <0.003 for both genes). Phylotypes 16S A and B were identified at all sampling sites. 16S A was more common than B at Deer Arm and South East Arm, and was rare at Neddy’s Harbour (Fig. 1). Phylotypes 16S C (N = 4) and RB 1 (N = 2) were found in the same bivalve specimens, all collected in Neddy’s Harbour, suggesting site restriction. RB 3 was found at Deer Arm (N = 4) and South East Arm (N = 8), but not at Neddy’s Harbour.

Specimens were collected over a span of four years, and hosts with symbiont phylotypes 16S A, 16S B, RB 2 and RB 3 were identified in all months and years (Table S1). Symbiont phylotypes 16S C and RB 1 were only identified in 2010 and 2011. Chi-square tests revealed that sampling date and phylotypes were statistically independent of each other (p-values > 0.7 for both genes). We found no apparent correspondence between host OTU and symbiont phylotype; the more common OTU 1 associated with all symbiont phylotypes, and the three individuals of host OTU 2 contained symbionts with 16S A/RB 2, 16S C/RB 1, and 16S C/unknown RB respectively. Finally, host size (Table S1) showed no obvious relationship with symbiont phylotype (chi-square p-value > 0.1 for both genes).

Symbiont diversity and spatial-temporal patterns

The greater number of samples analyzed herein has led to a revision of the spatial distribution of 16S rRNA phylotypes since Batstone & Dufour (2016), notably through a greater 16S rRNA phylotype diversity at Deer Arm and South East Arm than previously reported. It is now apparent that 16S A and 16S B are widespread in Bonne Bay, while 16S C has only been identified in Neddy’s Harbour. Some RuBisCO phylotypes occurred at multiple sampling sites, but RB1, found only in hosts having 16S C, appeared restricted to the shallower Neddy’s Harbour site. The restriction of 16S C/RB 1 to a single site supports the symbiont environmental acquisition mode proposed for thyasirids (Duperron et al., 2012; Dufour et al., 2014).

The ecotype hypothesis proposes that bacterial strains assemble in relation to environmental factors such as sedimentary organic matter, grain size, and sulfur content. Therefore, in environmentally acquired symbionts, the distribution of symbiont phylotypes will reflect these habitat characteristics (Brissac et al., 2016). This hypothesis may explain the apparent absence of RB 3, and the possible restriction of the 16S C/ RB 1 phylotype in Neddy’s Harbour, which is the shallowest sampling site, with the lowest organic matter content and coarsest sediments (Batstone & Dufour, 2016). The Deer Arm and South East Arm sites have similar water depths and sediment characteristics (Batstone & Dufour, 2016) that may be conducive to all phylotypes but 16S C/ RB 1. The more widespread RB 2 phylotype may be able to function within a broader range of environmental conditions than RB 1 and RB 3. Environmental patchiness on mm- to cm-scales can explain why thyasirids housing different symbiont phylotypes were found within the same Peterson grab sample (Table S1).

Symbiont evolution and the relationship between 16S and RuBisCO phylotypes

The Thyasira cf. gouldi symbiont 16S rRNA and RuBisCO phylogenies are not entirely congruent, although some conservation of gene pairs was identified (Table 2). Notably, while the three symbiont 16S phylotypes are closely related and form a single cluster, RuBisCO genes within the same symbionts are more evolutionarily diverged: RB 1 and RB 2 cluster together and with other symbionts, while RB 3 clusters with free-living bacteria (Fig. 3). A similar pattern in associated 16S rRNA and RuBisCO gene sequences has been noted in autotrophic proteobacteria and cyanobacteria, and attributed to horizontal transfer of the RuBisCO gene amongst phylogenetic lineages (Delwiche & Palmer, 1996; Elsaied & Naganuma, 2001; Kleiner, Petersen & Dubilier, 2012). During their free-living existence, symbionts are exposed to other bacteria in the sediment could undergo horizontal gene transfer events (Dahlberg, Bergström & Hermansson, 1998; Davison, 1999).

While the slight variation in 16S rRNA genes amongst Thyasira cf. gouldi symbionts may not be reflective of physiological differences, RuBisCO gene variants may be biologically significant, providing fitness benefits to symbionts under particular environmental conditions. The horizontal transfer of genes in bacterial symbionts may increase their metabolic efficiency, and therefore thyasirid symbioses may be particularly flexible by acquiring symbionts that are locally adapted to their microenvironment. Bathymodiolus septemdierum, another bivalve with environmentally acquired sulfur oxidizing symbionts, shows symbiont genomic variation linked with differences in metabolic capabilities, thought to be the result of horizontal gene transfer (Ikuta et al., 2016).

Host-symbiont interaction

As observed previously, the Thyasira cf. gouldi OTUs present in Bonne Bay do not show co-speciation with their symbionts; rather, both host OTUs 1 and 2 can form symbioses with a restricted, but diverse group of bacteria present in the environment (Batstone & Dufour, 2016).

Further, the identity and position of double-peaks in our sequence alignments suggest that some T. cf. gouldi individuals may host heterogeneous symbiont populations (i.e., they may show multi-infection, or have mixed symbiont populations, either consisting of multiple strains within a phylotype, or multiple phylotypes; as in (Vrijenhoek, Duhaime & Jones, 2007). Therefore, the symbiont population of thyasirids may resemble that observed in some mytilid species with extracellular, and some with intracellular symbionts (Won et al., 2003; Ikuta et al., 2016), as opposed to the clonal symbiont populations found in lucinids (Brissac, Merçot & Gros, 2011; Brissac et al., 2016).

In bivalves that acquire symbionts from their environment, selectivity may take place as gills undergo development. In lucinids, aposymbiotic juveniles may pick up multiple symbiont strains, which are maintained in undifferentiated cells dispersed throughout the lateral zone of gill filaments (Brissac, Merçot & Gros, 2011). These gill cells later differentiate into mature bacteriocytes (Gros, Frenkiel & Mouëza, 1997). Theoretically, from this mixed infection, the bacteriocytes with the best energetic yield are kept and imprinted with that bacterial strain throughout their lifetime, coupling the adult host with a specific strain of symbiont (Gros et al., 2012; Brissac et al., 2016). In contrast, some Thyasira cf. gouldi specimens appear to be co-infected by multiple symbiont strains, even as adults (shell sizes of hosts with degenerate bases in the sequence are not smaller than those of hosts having no sequence ambiguity; Table S1). Thyasirids may be capable of acquiring new and genetically mixed symbionts over their lifetime, potentially increasing metabolic fitness of the holobiont in changing environments. The burrowing behaviour of thyasirids, combined with symbiont uptake on the mucociliary surface of the foot with subsequent transfer to gills, is a possible mechanism for the acquisition of new symbionts over the course of the host’s life (Dufour et al., 2014; Zanzerl & Dufour, 2017).

This work highlights the importance of looking past 16S rRNA diversity when investigating symbiont populations, and suggests a greater degree of extracellular symbiont diversity in thyasirids than previously recognized. The ability of Thyasira cf. gouldi to associate with different symbiont strains may lead to improved fitness within environmentally variable habitats and contribute to its phylogenetic diversity. Further work should examine differences between T. cf. gouldi symbiont phylotypes and symbiont population heterogeneity in greater detail, through use of clone libraries, genomic, or proteomic investigations (e.g., cloning, next-generation sequencing, in situ hybridization). Similar work in other hosts with extracellular symbionts would also be warranted for a more comprehensive understanding of specificity and selectivity in chemosynthetic symbioses.