Extensive Thioautotrophic Gill Endosymbiont Diversity within a Single Ctena orbiculata (Bivalvia: Lucinidae) Population and Implications for Defining Host-Symbiont Specificity and Species Recognition

Symbiont diversity and host/symbiont functions have been comprehensively profiled for only a few lucinid species. In this work, unprecedented thioautotrophic gill endosymbiont taxonomic diversity was characterized within a Ctena orbiculata population associated with both seagrass- and alga-covered sediments. Endosymbiont metabolisms included known chemosynthetic functions and an additional conserved, previously uncharacterized C1 oxidation pathway. Lucinid-symbiont associations were not species specific because this C. orbiculata population hosted multiple endosymbiont strains and species, and other sympatric lucinid species shared overlapping symbiont 16S rRNA gene diversity profiles with C. orbiculata. Our results suggest that lucinid-symbiont association patterns within some host species could be more taxonomically diverse than previously thought. As such, this study highlights the importance of holistic analyses, at the population, community, and even ecosystem levels, in understanding host-microbe association patterns.

highlights the importance of holistic analyses, at the population, community, and even ecosystem levels, in understanding host-microbe association patterns. KEYWORDS host-microbe interactions, lucinid, metagenomics, metatranscriptomics, symbiosis C hemolithoautotrophic symbiosis, where symbiotic bacteria use inorganic chemical energy for the synthesis of organic compounds that benefit their hosts, is prevalent in marine bivalves, including the Lucinidae (1,2). All extant lucinid species examined to date host chemosynthetic bacterial endosymbionts belonging to the class Gammaproteobacteria within specialized epithelial gill cells known as bacteriocytes (3). Lucinid gill endosymbionts possess a diverse and varied suite of functions, including thioautotrophy (4), aerobic respiration (5), assimilatory and dissimilatory nitrate reduction (6), mixotrophy (7,8), hydrogenotrophy (7,8), and diazotrophy (7,9), which enable lucinids to colonize otherwise uninhabitable conditions (10). Specifically, lucinids burrow into sediments where oxygen is available from the oxic water column and dissolved sulfide is present due to the decomposition of organic matter via sulfate reduction (2,3,11). In modern seagrass habitats, a three-way symbiotic association among bacteria, lucinids, and seagrass has been proposed (11,12). Lucinid hosts acquire oxygen for respiration from oxygenated water surrounding seagrass roots, and the thioautotrophic endosymbionts oxidize sulfide while fixing inorganic carbon for the host (12). Seagrass growth is promoted because the endosymbionts remove sulfide and have the potential to fix nitrogen (7,11,12). This symbiotic system likely developed in the late Cretaceous period, based on the correspondence of adaptive radiation of shallow marine lucinid species in the fossil record to the evolutionary first appearances of marine angiosperms, including seagrasses (13).
Lucinid gill endosymbionts are related to a larger group of diverse marine thioautotrophic symbionts (1). However, unlike chemosymbiotic Solemyidae and Vesicomyidae bivalves, where vertical or mixed symbiont transmission has been observed (14)(15)(16), lucinids studied to date acquire their endosymbionts environmentally (17)(18)(19). Based on their 16S rRNA gene sequences, lucinid endosymbionts occur in three distinct clades, with the largest (clade A) originating from hosts that primarily inhabit sediments from seagrass beds and two (clades B and C) originating from lucinids that are more likely to inhabit sediments from mangrove-rich areas (20). It is possible that each clade represents a separate species (8), although an insufficient number of lucinid species has been examined at a sufficiently high genetic resolution to test this hypothesis. At present, low to no variability among clade A lucinid endosymbiont 16S rRNA gene sequences suggests that this group belongs to a single species (20)(21)(22)(23). Earlier studies of gill thioautotrophic endosymbionts from lesser Antillean lucinid hosts, including Ctena imbricatula (referred to as Codakia orbiculata in reference 21), Codakia orbicularis, Parvilucina pectinella (referred to as Codakia pectinella in reference 21), Anodontia alba, Divalinga quadrisulcata (referred to as Divaricella quadrisulcata in reference 21), and Lucina roquesana (referred to as Linga pensylvanica in reference 21), show identical 16S rRNA gene sequences (22,23) and the capability of infecting aposymbiotic Codakia orbicularis juveniles (24). Recent reanalysis of these thioautotrophic endosymbionts using five other marker genes (instead of the slow-evolving 16S rRNA gene) revealed that intraspecific symbiont strain diversity in the population was shaped by host geographic location (25). Strain-specific symbiont acquisition was also observed in the same study (25), where starved Ctena imbricatula individuals (referred to as Codakia orbiculata in reference 21) could reacquire only the exact symbiont strain that they initially hosted before starvation. This suggests that, unlike undifferentiated naive bacteriocytes, mature bacteriocytes can possibly recognize and select specific chemolithoautotrophic symbiont strains (24,25).
In spite of these findings, species-and strain-level diversity within lucinid endosymbiont clades remains unverified for most lucinid species. To date, endosymbiont functions for two clade A endosymbionts, "Candidatus Thiodiazotropha endolucinida" from Codakia orbicularis (9) and "Ca. Thiodiazotropha endoloripes" from Loripes orbiculatus (syn ϭ Loripes lucinalis) (7), as well as for the clade C endosymbiont from mangrovedwelling Phacoides pectinatus (8), have been investigated using omics approaches. However, limited attention was given to species-and strain-level diversity among the endosymbiont populations. Information about endosymbiont species-and strain-level diversity will be useful for interhost and interpopulation comparisons and could provide new insight into host-symbiont association patterns and possibly biogeographical or environmental drivers of diversity among free-living bacteria in the environment, prior to being acquired by the host or becoming part of the gill microbiome.
In this study, we focused on characterizing the taxonomic, genetic, and functional composition of endosymbiont communities within Ctena orbiculata (Montagu, 1808) living sympatrically with at least five other lucinid species in carbonate tidal flat sediment at Sammy Creek Landing, Sugarloaf Key, FL (USA). We selected C. orbiculata as the main study organism because it dominated the community. The primary goal was to investigate possible influences of environmental and spatial factors, such as occupation of seagrass-or alga-covered sediment, on endosymbiont diversity, host-symbiont functions, and strainlevel diversity by comparing symbionts from multiple individual hosts of the same species. We used C. orbiculata gill microbiomes and metagenomes to generate bacterial taxonomic profiles and metagenome-assembled genomes (MAGs) that could test the hypothesis that strain-level C. orbiculata endosymbiont diversity exists within a population. Prior research had not been at a sufficiently high genetic resolution to uncover this level of endosymbiont diversity. We also used metatranscriptomic analyses to infer host-symbiont gene expression and to identify differentially expressed genes across C. orbiculata endosymbiont taxa and the environment.

RESULTS
Site characterization. A volume of nearly 3 m 3 of carbonate sediment was hand dug and hand sieved from 13 quadrats covered with either seagrass or algae ( Fig. 1; see also Table S1 in the supplemental material) to reveal a diverse lucinid community that averaged 45.5 live lucinids per m 3 , predominately from the top 30 cm of sediment. We observed an average of 46.4 live C. orbiculata specimens from the site quadrats (n ϭ 13), 22.3 live Codakia orbicularis specimens, 6.2 live Lucinisca nassula specimens, and 4.1 live Anodontia alba specimens (Table S1). A few live individuals of Parvilucina pectinella (1.5) and Radiolucina amianta (0.8) were also recovered at the site. Lucinids were encountered from sediment that had a polymodal grain size distribution of very poorly sorted to poorly sorted coarse sand to gravelly sand with an average of 3.3% organic carbon content, regardless of sediment depth. Seagrass species encountered in the quadrats included Halodule wrightii (shoal grass), Syringodium filiforme (manatee grass), and Thalassia testudinum (turtle grass), as well as rhizophytic calcareous algae, such as Penicillus spp., Halimeda spp., other green algae, and the red algae Goniolithon spp. (Table S1). Porewater dissolved sulfide concentrations ranged from Ͻ30 mol/liter to 2.9 mmol/liter, and dissolved oxygen concentrations were lower (Յ31 mol/liter) than ocean water samples for all sampled quadrats (Table S1). Dissolved methane concentrations were generally higher near the shoreline, with the highest concentration being 2.18 mol/liter from the same quadrat as the highest sulfide measurement.
Gill microbiome diversity. From the 13 quadrats, four were chosen for detailed examination of lucinid and symbiont diversity based on seagrass or alga coverage (Fig. 1); two quadrats were dominated by Halodule wrightii and two others by Halimeda spp. To our knowledge, a potential association between lucinids and calcareous green algae has not been established previously. The V4 region of the 16S rRNA gene was sequenced from a total of 36 C. orbiculata gill samples retrieved from the four quadrats, as well as from three L. nassula gill samples, two A. alba gill samples, and one Codakia orbicularis gill sample. Of these hosts, 10 Ctena orbiculata gill samples were dissected from the same bivalves but processed and sequenced at two different institutions (see Materials and Methods), thus serving as technical replicates (Fig. 2). Sequences were clustered into operational taxonomic units (OTUs) at 99% sequence identity (26) and showed, on average, 98% Ϯ 2% Good's sequencing coverage (27). Taxonomic classification of resulting OTUs against the Silva v132 (28) database at 90% bootstrap confidence revealed six coexisting "Candidatus Thiodiazotropha"-like OTUs (OTU1 to OTU6) present at Ͼ60% relative abundances in at least one gill sample from the mixed lucinid community (Fig. 2). OTU1 dominated the gills of 18 of 26 C. orbiculata individuals (average relative abundance of 93% Ϯ 8%) and a single Codakia orbicularis individual from the site (99% relative abundance; Fig. 2). The representative sequence for OTU1 was also 100% identical to sequences from two Codakia orbicularis individuals previously found in the Bahamas (Fig. S1) (7). OTU2 was predominant in the gills of three Ctena orbiculata individuals at 85% Ϯ 7% average relative abundance and two A. alba individuals at 96% Ϯ 0.3% average relative abundance, while OTU3 dominated the gills of three C. orbiculata individuals at 95% Ϯ 4% average relative abundance (Fig. 2). The OTU3 representative sequence matched the previously sequenced V4 region of the 16S rRNA gene from L. nassula off Guadeloupe (22) (Fig. S1). OTU4 dominated the gills of three of four L. nassula individuals (80% Ϯ 7% average relative abundance) and one C. orbiculata individual (97% relative abundance) but was also detected at 19% relative abundance in another C. orbiculata individual and 0.6% Ϯ 0.4% relative abundance in six additional C. orbiculata individuals (Fig. 2). The representative sequence for OTU4 was identical to sequences retrieved from P. pectinella off Guadeloupe (GenBank accession no. MK346268) and L. nassula from Lemon Bay, FL, USA (29) (Fig. S1). OTU5 occurred at 69% relative abundance in the gills of one C. orbiculata individual and at 1% Ϯ 4% average relative abundances in 17 other individuals (Fig. 2). The OTU distribution did not follow any clear spatial trend. For instance, high relative abundances of OTU1 to OTU3 were identified in gill specimens from both seagrass-and alga-covered quadrats, and OTU4 was detected at high relative abundances (Ͼ50%) in gill specimens from a seagrass-covered quadrat (T20 at 40 m) but low relative abundances in the other three quadrats (Table S1). Incidentally, as previously observed in studies on lucinid symbiont diversity (8,(30)(31)(32), Endozoicomonas-like (order Oceanospirillales) and Spirochaeta-like bacterial OTUs were also present in lucinid gill tissues from this site (discussed in Text S1).
Metagenomic read coverage profiles of each representative OTU-specific MAG, bacterial replication rates estimated by iRep (34) from MAG data, and percentages of metatranscriptomic reads mapped to protein-encoding genes of each representative OTU-specific MAG were generally consistent with relative abundance patterns of their corresponding 16S rRNA gene OTU (Fig. 4). The only exception was OTU2-dominated gill metatranscriptome 22B from quadrat T22 (3.5 m), which showed higher percentages of reads mapped to OTU1 than to OTU2 (Fig. 4C). Symbiont-specific gill transcriptomes of OTU1-/OTU3-dominated specimens clustered together with an 0.9 Ϯ 0.07 average pairwise Pearson correlation coefficient (PCC), while the OTU2-dominated symbiont transcriptome 4D from quadrat T20 (40 m) and OTU4-dominated symbiont transcriptome 4F from quadrat T20 (40 m) appeared to be outliers sharing Ͻ0.4 PCC with the other specimens (Fig. 5). Symbiont OTU-specific patterns were not observed across the entire gill metatranscriptomes (average 0.7 Ϯ 0.4 PCC between samples) or another subset of Mollusca-related transcriptomes (average 0.8 Ϯ 0.06 between samples). Host 18S rRNA, 28S rRNA, and mitochondrial cytochrome b (cytb) gene sequences extracted from unbinned C. orbiculata gill metagenomes clustered unambiguously with reference sequences from C. orbiculata (Fig. S2), thereby confirming the host taxonomy of these specimens.
Core symbiont functions. Pangenomes of "Ca. Thiodiazotropha"-like MAGs were predicted by Rapid Annotation using Subsystem Technology (RAST) (35) to share ϳ62% gene and ϳ83% subsystem content. As noted by Lim et al. (8), gene and subsystem annotations were based on incompletely sequenced and annotated MAGs. As such, the numbers of shared genes and subsystems presented here are imprecise estimates that cannot account for missing, unbinned, or unclassifiable genes or for incomplete pathways and potential strain/cross-species contamination of the MAGs. Moreover, limitations of host-symbiont metatranscriptomic analyses also apply (8).
Ctena orbiculata symbionts also encoded and expressed genes for the valine-glycine repeat protein G (VgrG; average 0.7 Ϯ 1 TPM) but not hemolysin-coregulated protein (Hcp) and VasL, which is exclusive to the type VI secretion system 2 gene cluster (39). In addition to thiotrophy and hydrogenotrophy-related genes, a C 1 oxidation gene cluster encoding proteins involved in pyrroloquinoline quinone (PQQ) synthesis (average 0.2 Ϯ 0.6 TPM), PQQ-dependent methanol oxidation (Mdh; average 3 Ϯ 3 TPM), and tetrahydromethanopterin (H 4 MPT)-dependent formaldehyde oxidation (average 0.2 Ϯ 0.6 TPM) was conserved in all sequenced C. orbiculata symbionts and "Ca. Thiodiazotropha endolucinida" (8) (Fig. 7). Downstream of this gene cluster, another formate oxidation gene cluster encoding NADH-quinone oxidoreductase subunit F and formate dehydrogenase alpha subunit (FdhA) was predicted in all C. orbiculata symbionts, "Ca. Thiodiazotropha endolucinida" (9), and "Ca. Thiodiazotropha endoloripes" (7). Phylogenetic analyses of Mdh and FdhA protein sequences showed OTU-specific clustering patterns across C. orbiculata symbionts (Fig. S4 and S5 and Text S1) consistent with the 16S rRNA gene tree (Fig. S1) and phylogenomic tree (Fig. 3A). These sequences were most closely related to each other and more distantly associated with the giant Teredinidae bivalve Kuphus polythalamia (40) and other free-living bacterial species (Text S1). Two sets of qPCR primers targeting Mdh from OTU1 and OTU2 were designed to amplify DNA and cDNA from gill specimens (Table S2), and the qPCR cDNA copy numbers of Mdh were generally consistent with TPM values observed in five out of seven amplified gill specimens (Fig. 7D). Other genes potentially related to C 1 oxidation were also annotated in MAGs and/or transcriptomes of C. orbiculata symbionts (Text S1). C 1 assimilation genes in the ribulose monophosphate (RuMP) pathway, and four out of 10 key genes (serine hydroxymethyltransferase, malyl coenzyme A lyase [malyl-CoA-lyase], malate dehydrogenase, and enolase) in the serine-glyoxylate cycle (41), were identified in the symbiont MAGs. These four genes, as well as 11 additional accessory genes assigned by RAST (35) to the serine-glyoxylate cycle subsystem, were also predicted in other carbon-related pathways, including biosynthesis, glycolysis, photorespiration, tricarboxylic acid cycle, glyoxylate cycle, polyhydroxybutyrate metabolism, acetyl-CoA fermentation to butyrate, and phenylalkanoic acid degradation.
Symbiont species/strain differences. In addition to interspecies differences in the forms of RuBisCO identified, C. orbiculata symbiont MAGs and transcriptomes showed very little interstrain and interspecies variation of complete or nearly complete metabolic pathways. For aerobic respiration, C. orbiculata symbionts and other clade A lucinid symbionts (7,9) potentially utilize cbb3 (average 4 Ϯ 5 TPM in the former) and aa3 terminal oxidases (average 5 Ϯ 9 TPM). Furthermore, OTU2-related MAGs and transcriptomes contained genes for cytochrome bd ubiquinol oxidase (average 0.5 Ϯ 0.4 TPM) that were also detected in the unbinned metagenomic assembly and metatranscriptome of OTU4-dominated gill specimen 4F from the seagrass-covered quadrat T20 (40 m) at 0.02 TPM, as well as the metatranscriptome of OTU3-dominated gill specimen 4C from quadrat T20 (40 m) at 0.3 TPM. Ctena orbiculata symbionts also likely use distinct types of clustered regularly interspaced short palindromic repeats (CRISPR)-associated genes (Text S1).

DE analyses of host-symbiont gene expression across quadrats.
Besides taxonspecific clustering patterns that were evident from 16S rRNA gene profiles and sequence content of the MAGs, OTU1-dominated symbiont-specific transcriptomes formed two separate subclusters based on whether specimens were collected from alga-covered and seagrass-covered quadrats ( Fig. 2 and 5), although these transcriptomes shared 0.9 Ϯ 0.01 average pairwise PCC with each other. Exploratory DE analyses on these transcriptomes showed only five symbiont-related DE genes predicted by one (DESeq2) (45) of four DE analysis algorithms used. Of these, the dissimilatory sulfite reductase transferase protein DsrC was upregulated in the alga-covered quadrat, but three transcript clusters that were upregulated in the seagrass-covered quadrat encoded hypothetical proteins and contained one transcript cluster homologous to a cell wall-associated hydrolase from an alphaproteobacterium. Host-related transcripts did not exhibit the same clustering patterns between alga-and seagrass-covered quadrats but instead showed 73 putative DE genes predicted by both DESeq2 (45) and edgeR (46) across OTU1-dominated metatranscriptomes in these quadrats. According to Uni-Prot (47) annotations, these host-related DE genes were involved in a variety of muscle-related, cytoskeletal, cochaperone, transcriptional, translational, and other functions. Notably, host-related transcripts encoding cytochrome b-c1 complex, mitochondrial succinate dehydrogenase, and mitochondrial sulfide:quinone oxidoreductase (Sqr) were predicted to be upregulated if originating from seagrass-covered quadrats.

DISCUSSION
To date, the taxonomic and functional diversities of lucinid symbiont species remain largely unexplored because, out of Ͼ100 identified lucinid species listed in the NCBI database (48), only three gill endosymbionts, each from a different host species, have been comprehensively sequenced (7)(8)(9). In spite of previous marker gene-based diversity studies (20,22,23,25), and the recently characterized gill microbiome, metagenome, and metatranscriptome from "Ca. Sedimenticola endophacoides" from P. pectinatus (8), in-depth investigations into symbiont variation within a single lucinid host population are currently lacking. Our evidence for endosymbiont species-and strainlevel taxonomic, genetic, and functional diversity within a single lucinid host expands previously reported findings of strain-level diversity within a single thioautotrophic endosymbiont species (25). Similar functionally homogenous but taxonomically heterogenous marine symbiont communities have been reported, for example, in the light organs of the squids Sepiola affinis and Sepiola robusta, where closely related Vibrio fischeri and Vibrio logei were detected at different abundances (49). The C. orbiculata gill metatranscriptomic analyses further reveal differentially expressed host and/or endosymbiont genes of different bacterial taxa. Although gill endosymbiont OTUs from C. orbiculata shared a high number of core genes and functions, which included a well-conserved and previously undiscovered C 1 oxidation pathway, the discovery of other genetic and functional differences among coexisting "Ca. Thiodiazotropha"-like endosymbiont taxa supports our hypothesis that strain-and even species-level endosymbiont diversity can exist within a single host population living under variable environmental conditions. Moreover, because we also identified "Ca. Thiodiazotropha"like OTUs among other sympatric lucinid species, genetic and functional trait similarities with C. orbiculata endosymbionts likely exist within other lucinid populations, which emphasizes the importance of future investigations of these chemolithoautotrophic symbiotic systems at broader population, community, and ecosystem levels.
The thioautotrophic C. orbiculata endosymbionts shared a high number of common genes and functions, including previously characterized lithoautotrophy, diazotrophy, and potential heterotrophy (7)(8)(9), as well as the oxidation of C 1 compounds, including methanol, formaldehyde, and formate. These potential C 1 compound oxidation functions were also discovered in MAGs of "Ca. Thiodiazotropha endolucinida" (9) that had been unrecognized and undescribed until our analyses. Ctena orbiculata endosymbionts may use C 1 compounds as an energy source because they encoded and expressed formate dehydrogenase O that enables the use of formate as an electron donor during respiration (50). Because there was no substantial genetic evidence for RuMP and serine-glyoxylate C 1 assimilation pathways, we hypothesize that the CO 2 end product of C 1 oxidation for these endosymbionts could be fixed via the autotrophic Calvin-Benson-Bassham cycle, as demonstrated in the diazotrophic alphaproteobacterial Xanthobacter strain 25a (51). Additionally, many transcripts involved in temperature, oxidative, envelope, and nutritional stress responses were highly expressed in C. orbiculata endosymbionts and likely reflect stresses caused by the external environment or stresses due to intracellular host selection mechanisms (8), the latter being similar to the Euprymna-Vibrio symbiotic system (52)(53)(54). Defense-related transcripts involved in type I, II, and VI secretion systems, such as colicin-related transcripts, general secretion (Sec)-related transcripts, and twin-arginine translocation (Tat)-related transcripts, may reduce symbiont-symbiont competition by killing closely related strains (8,55), as well as contributing to host infection (56,57). Also, differentially expressed genes for secretion, stress response, transport, and biosynthesis proteins across the C. orbiculata endosymbiont taxa may also be due to variable habitat conditions. For aerobic respiration, C. orbiculata endosymbionts and other seagrass-associated lucinid symbionts (7,9) can potentially use the high-affinity cbb3-type and low-affinity aa3-type terminal oxidases under low and high oxygen concentrations (58,59), respectively. Additionally, cytochrome bd ubiquinol oxidase, adapted to microaerobic environments (60), was encoded and expressed mainly by the "Ca. Thiodiazotropha"-like OTU2. This enzyme may facilitate symbiotic nitrogen fixation via oxygen scavenging and respiratory protection, as demonstrated in other free-living and plant-associated nitrogen-fixing bacteria (61)(62)(63), along with glutathione peroxidase that confers similar protective functions by scavenging H 2 O 2 , as shown in the legume-Rhizobium symbiosis (64). Conversely, RuBisCO form Iaq, detected in "Ca. Thiodiazotropha"-like OTU1 to OT4, is more efficient than form II, detected in OTU2 and OTU4, at distinguishing between the competing substrates oxygen and CO 2 (65). Differences in aerobic oxidase and RuBisCO forms, as well as highly expressed antioxidant glutathione peroxidase transcripts, indicate that the symbionts experience higher or variable levels of oxygen in their environments, which could be during their intracellular free-living or extracellular host-associated growth stages. Analyses of host-related transcripts revealed 100-foldhigher expression of sulfide-reactive hemoglobin 1 (66) in C. orbiculata compared to oxygen-reactive hemoglobins 2 and 3 (66). This is in contrast with previous observations in P. pectinatus, where consistently similar high expression levels of hemoglobins 1, 2, and 3 were detected (8). Compared to alga-covered quadrats, C. orbiculata host-related transcripts in seagrass-covered quadrats showed candidate upregulated genes that encoded cytochrome b-c1 complex involved in aerobic respiration and oxidative stress-triggered apoptosis (67), mitochondrial succinate dehydrogenase that connects the tricarboxylic acid cycle to the electron transport chain (68), and mitochondrial Sqr that oxidizes sulfide to thiosulfate (69). Collectively, these findings indicate that complex host-symbiont cooperative regulation of intracellular oxygen and sulfide levels facilitates thioautodiazotrophic processes, which are distinct from the P. pectinatus chemolithoautotrophic symbiotic system (8).
Unlike previous studies on Caribbean lucinid symbiont species (25) and Sepiola squid symbioses (70), where species distribution is influenced by host geography and temperature, respectively, we did not uncover extensive spatial controls or apparent nonrandom patterns of strain-level C. orbiculata endosymbiont species distribution. The observed endosymbiont diversity from C. orbiculata and other sympatric lucinid species could be due to constant symbiont turnover, because adult C. imbricatula clams (referred to as Codakia orbiculata in reference 21) have been reported to lose their chemosynthetic endosymbionts when starved for sulfide and reacquire endosymbionts from their natural habitat. This behavior points to continuous horizontal symbiont acquisition throughout the life span of this host species (71,72). However, a recent cross-infection experiment demonstrates that starved C. imbricatula (referred to as Codakia orbiculata in reference 21) adults with mature bacteriocytes could reacquire only symbiont strains that they initially hosted (25). In contrast, Codakia orbicularis juveniles, presumably with naive bacteriocyte precursors, acquire symbionts from purified gill-symbiont sections of their own species, as well as other lucinid species hosting thioautotrophic symbionts with identical 16S rRNA gene sequences (24,25). On the other hand, starvation-related symbiont loss has not been observed in Lucina roquesana (referred to as Lucina pensylvanica in references 21 and 73). The contrasting evidence for symbiont acquisition, loss, and reacquisition suggests that such mechanisms could be host species specific in lucinid bivalves (73). Based on our results, it does not appear that a single symbiont species occupies a single host. Instead, a single host may harbor multiple symbiont species and strains. This is likely due to either higher taxonomic symbiont diversity encountered by the host in the environment (20) or less stringent host regulation of symbiont acquisition in the mixed lucinid community (25).
Interhost symbiont transmission could also be a possible reason for the observed symbiont diversity in the mixed lucinid community. However, previous studies on C. imbricatula (referred to as Codakia orbiculata in reference 21) and Codakia orbicularis ruled out the possibility of host-to-host symbiont transfer by showing that adults do not release their thioautotrophic symbionts into the environment (74). In addition, opportunistic host-microbe associations have been proposed for lucinid bivalves (20) and observed in the green macroalga Ulva australis (75). Therefore, symbiont acquisition by this mixed lucinid community may be based on random encounters during the juvenile stage, as well as the diversity and distribution of "Ca. Thiodiazotropha" species in this habitat, which are under investigation. Nevertheless, our results suggest that a taxonomically heterogenous symbiont community, rather than a monospecific, homogenous symbiont community, exists in some lucinid species.
Overall, our study uncovered remarkable taxonomic, genetic, and functional thioautotrophic gill endosymbiont diversity in C. orbiculata and furthers the current understanding of lucinid-symbiont association patterns, physiology, and interactions. Despite the small sample size and replication limitations that potentially affect the statistical significance of DE analysis, as well as the lack of pure symbiotic monocultures, our findings highlight the intriguing, poorly understood complexity of lucinidbacterium chemolithoautotrophic symbioses and generate a range of new testable hypotheses to evaluate the establishment, persistence, stability, and distribution of symbiont communities within hosts and their environment. Future studies to understand symbiont taxonomy, intracellular spatial distribution, functional diversity, and host-symbiont specificity across environmental gradients should include cross-infection experiments, imaging experiments involving fluorescence in situ hybridization, controlled aquarium experiments, and large-scale field studies coupled with -omics analyses. This work also proposes future in-depth investigations on C 1 metabolism in thioautotrophic lucinid symbionts and three-way interactions between the environment, lucinid hosts, and their symbionts. Additionally, targeted experimental approaches will be useful in examining similarities and differences in phenotypic effects that various symbiont taxa and lucinid hosts exert on each other. In each of the previous -omics studies of lucinid chemolithoautotrophic symbiotic systems, a unique species name for the lucinid endosymbiont has been proposed, such as the two clade A endosymbionts, "Ca. Thiodiazotropha endolucinida" from Codakia orbicularis (9) and "Ca. Thiodiazotropha endoloripes" from Loripes orbiculatus (7). For the clade C endosymbiont from the mangrove-dwelling P. pectinatus (8), which does not fix nitrogen and instead belongs to the genus Sedimenticola, the candidate name for the endosymbiont was justifiably proposed as "Ca. Sedimenticola endophacoides" (8). However, based on our analyses here, we think it is too premature to assign a species name for the C. orbiculata symbionts. Even if "Ca. Thiodiazotropha endolucinida" is likely appropriate for the OTU2 and OTU4 strains, the OTU1 and OTU3 strains likely belong to a different species, but this species is not specific to C. orbiculata or even to lucinid hosts having clade A endosymbionts. Therefore, we contend that future lucinid symbiont taxonomic assignments and classification need to be based on detailed investigations of more lucinid symbionts from all three endosymbiont clades and from different habitats.

MATERIALS AND METHODS
Study area and sediment characterization. Field sampling was done on the carbonate tidal flat at Sammy Creek Landing, on Sugarloaf Key, FL (USA), from 27 June to 1 July 2016. Twelve quadrats were separated by 5 m along two 50-m-long transects (six quadrats each) that ran perpendicular to the shoreline; a 13th quadrat was established 3.5 m away from shore to accommodate equal sampling of seagrass versus benthic green and red algal cover (Fig. 1). The seagrass or algal cover was estimated for each quadrat before hand digging and hand sieving the sediment to uncover lucinids and other infauna. In general, depths of sediment excavation varied, from 0.29 to 0.5 m, and generally stopped at a very coarse, organic-material-rich sediment layer below which lucinids were not recovered. Intermixed with the sediment was abundant road (e.g., concrete and asphalt) and other anthropogenically sourced debris, such as manufactured wood, roofing material, and mixed-use plastic, likely due to past storm activity that impacted the shoreline and tidal flat. Methods for geochemical and sedimentological analyses are in Text S1 in the supplemental material.

Diversity and Specificity of Lucinid Symbionts
Sample collection and sequencing. Prior to dissection, all live-collected lucinid specimens were temporarily stored and transported in Whirl-Pak bags (Nasco, Fort Atkinson, WI, USA) at ambient temperature filled with surface water from the habitat (8). Dissection was performed within 30 min or within the same day of collection (explained below), where gill and foot tissues were separated from other body tissues and preserved in RNAlater (8). Dissected gill tissues were divided for molecular and sequence analyses between Clemson University (CU; Clemson, SC, USA) and University of Tennessee-Knoxville (UTK; Knoxville, TN, USA) such that, for each specimen, one fixed gill tissue went to each institution or both gill tissues went to UTK. Valves and the remaining preserved tissues were archived and cataloged at the Museum of Geology at the South Dakota School of Mines and Technology (Rapid City, SD, USA). This study focused on specimens from four of 13 quadrats that were predominantly covered by either seagrass or algae: T20 (20 m), T20 (40 m), T21 (0 m), and T22 (3.5 m) (Fig. 1). Ctena orbiculata specimens from these four quadrats used for 16S rRNA gene, metagenomic, and metatranscriptomic analyses at CU were preserved within 30 min of collection. Nucleic acids extraction, fluorometric DNA/RNA quantification, cDNA synthesis, and 16S rRNA gene library preparation from this collection were performed as described in the work of Lim et al. (8). At CU, the V4 region of the 16S rRNA gene from 24 C. orbiculata gill tissues, as well as metagenomic libraries prepared from DNA extracted from eight C. orbiculata gill samples using NEBNext dsDNA Fragmentase plus the NEBNext Ultra II DNA library prep kit for Illumina or the NEBNext Ultra II FS DNA library prep kit for Illumina (New England Biolabs, Ipswich, MA, USA), were sequenced on the Illumina MiSeq V2 2-by 250-bp platform (San Diego, CA, USA). RNA extracted from 11 C. orbiculata gill samples at CU was prepared for metatranscriptomic sequencing, using procedures in the work of Lim et al. (8), on the Illumina HiSeq 4000 2-by 150-bp platform at the Duke Center for Genomic and Computational Biology (Durham, NC, USA). All library concentrations were quantified with the Qubit dsDNA HS assay (Life Technologies, Austin, TX, USA), and their insert sizes were determined with the Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). Specimens of additional C. orbiculata individuals (n ϭ 3) and other sympatric lucinid species (n ϭ 7) from the same four quadrats were also used for only 16S rRNA gene analysis at UTK. These specimens were preserved within the same day of collection. DNA was extracted from the gill tissues of 12 C. orbiculata (including 9 of the gill samples fixed by CU), four L. nassula, two A. alba, and one Codakia orbicularis specimen using methods described in the work of Lim et al. (8). 16S rRNA gene V4 region libraries of these samples were prepared and sequenced by Molecular Research LP (Shallowater, TX, USA) using the Illumina MiSeq V3 2-by 300-bp platform.
Data analyses. 16S rRNA gene sequence reads from CU and UTK were combined, quality trimmed, and processed with MOTHUR v1.41.1 (76) using the pipeline documented in the work of Lim et al. (8). Briefly, all reads were quality trimmed at Q ϭ 25 using Mothur's trim_seq and remove_seq commands and processed according to MOTHUR's MiSeq SOP available at https://www.mothur.org/wiki/MiSeq_SOP. OTUs clustered at 99% sequence identity were taxonomically classified using the Silva v132 database (28) with 90% bootstrap confidence. The data set was subsampled to 6,676 sequences, which was the number of sequences in the smallest library. OTU relative abundances and Good's coverage (27) within each sample were calculated using the MOTHUR get.relabund and summary.single commands, as described in the work of Lim et al. (8). Representative sequences of the 10 most abundant OTUs were obtained using MOTHUR's get.otu rep command, which identifies a representative sequence sharing the smallest maximum or average distance (in the event of a tie) to other sequences classified within the same OTU. Reads from all eight metagenomic libraries were trimmed and assembled individually (8). Additionally, two metagenomic libraries of gill specimens dominated by OTU1 and two libraries of gill specimens dominated by OTU2 were coassembled separately. Read mapping, binning, MAG quality assessment, MAG annotation, and ANI and AAI calculations were performed as detailed in the work of Lim et al. (8). Methods for phylogenetic analyses are outlined in Text S1 in the supplemental material. Bacterial replication rates were estimated using the iRep software (34) by mapping metagenomic reads to representative OTU-specific MAGs (22G, 22B ϩ 4D, 21D, and 4F) with Ն75% completeness and Յ2% contamination, using Bowtie 2 v2.3.4.1 (77) in very sensitive local and dovetail mode and SAMtools v1.7 (78). For metatranscriptomic analyses, a pangenome for each symbiont OTU was created using the method in the work of Lim et al. (8). De novo metatranscriptomic assembly, transcript cluster (gene) abundance estimation, count normalization, transcript-to-gene mapping, and transcript annotation were performed using Trinity v2.6.6 (79), Trinotate v3.1.1 (https://trinotate.github.io/), and web and local Blast searches (48), as described in the work of Lim et al. (8). For symbiont abundance estimation, trimmed reads were mapped to representative MAGs using the Bowtie 2 (77) -no-unal -no-mixed -no-discordant -gbar 1000 -end-to-end -k 200 options. Nonchimeric fragments mapped to protein-encoding genes were calculated with featureCounts v1.5.2 (80) using the -c and -p options. DE and functional gene ontology (GO) enrichment analyses were performed on a raw read count matrix processed with Trinity's remove_batch_ef-fects_from_count_matrix.pl script. Batch-removed read counts of transcript clusters mapped to thioautotrophic symbiont MAGs and those with eukaryotic homologs were analyzed separately using DESeq2 (45), edgeR (46), Reproducibility-Optimized Test Statistic (ROTS) (81), voom (82), and GOSeq (83) software incorporated within Trinity. DESeq2 (45) and edgeR (46) rely on the negative binomial distribution for DE analyses but differ in their normalization and statistical testing methods. ROTS uses properties from gene expression data to optimize a modified t statistic that maximizes the overlap of top-ranked genes in group-preserving bootstrapped data sets (81), while voom models the mean-variance relationship of depth-normalized log counts nonparametrically, incorporates it into a precision weight for each individual observation, and then fits a linear model to the weighted data using the limma package (84) for DE analyses (82). On the other hand, GOSeq (83) identifies Gene Ontology (GO) (85) terms, extracted from Trinotate's annotations, that are overand underrepresented in DE gene data sets predicted by each software program. To identify DE genes and