Bacterial composition in the toheroa (Paphies ventricosa), a threatened surf clam from Aotearoa (New Zealand)

The toheroa (Paphies ventricosa) is an Aotearoa (New Zealand) endemic surf clam that remains threatened following population collapse due to overfishing in the twentieth century. Despite protective measures being in place for more than 4 decades, toheroa populations have inexplicably failed to recover. As part of an investigation into the possible role of disease in preventing their recovery, an exploration of the bacterial composition in toheroa was conducted over their entire geographic range. The bacterial composition in toheroa tissues was dominated by Spirochaetaceae, Mycoplasmataceae, and Endozoicomonadaceae, and varied at both large (between geographically separated sites) and small spatial scales (beds < 10 km apart). At small scales, it was habitat, in this case the presence or absence of freshwater outflows, which appeared to be a major influence on bacterial composition. Given that the decline of toheroa has also coincided with changes in land use that have reduced the amount of freshwater reaching the toheroa beaches, it is possible that habitat-related shifts in the abundance of certain bacterial symbionts are affecting the health and impeding recovery of this iconic and culturally significant species.


Introduction
The toheroa (Paphies ventricosa) is a large intertidal surf clam, typically found on exposed high-energy surf beaches. On many beaches where toheroa were once incredibly abundant, they are now absent (Muriwai and Foxton), or populations are severely diminished (Ripiro and Te-Oneroa-a-Tōhē, see Fig. 1) (reviewed in Ross et al. 2018a). Unsustainable harvesting practises (recreational and commercial) of the early twentieth century contributed to, or caused, the collapse of toheroa populations, but do not explain their failure to recover despite more than four decades of management and protection. Habitat loss and disease have been identified as potential factors preventing population recovery (Ross et al. 2018b). While toheroa are often found at higher densities in close proximity to freshwater streams (seeps or overland outflows) (Beentjes 2010;Williams et al. 2013), the physiological reasons for this habitat 'selection' are not fully understood. Recent studies have suggested thermal protection from desiccation (Cope 2018) or a basin-like beach topography at stream sites that concentrates toheroa spat and phytoplankton (Williams et al. 2013;Cope 2018;Ross et al. 2018a) as explanations for their abundance in these areas. Consequently, the loss of freshwater streams and associated topography as a consequence of exotic pine forestry (Williams et al. 2013) and water extraction for agricultural purposes (Ross et al. 2018a) may be affecting remaining toheroa populations in this region. Considering the influence of freshwater outflows on toheroa health is unknown, the implications of this habitat modification for toheroa health and remaining toheroa populations is also a mystery.
In hosts, bacteria regulate resistance to, and contribute to, disease (Bass et al. 2019), and many biogeochemical cycles within host cells are driven by symbiotic bacteria (Neave et al. 2016;Tandon et al. 2020). In this way, bacteria contribute to host survival and homeostasis. Bacteria are drivers and modulators of numerous metabolic processes within their hosts. In corals, bacteria perform functional roles in antimicrobial activity (Morrow et al. 2015), bleaching defence, and various biogeochemical processes (reviewed in Bourne et al. 2016). In other marine invertebrates, metabolic processes are underpinned by gut microbiota and their respective functional roles. For this reason, health studies of commercially exploited species often consider bacterial composition in the gut as a proxy for host health. For example, gut microbiota studies have been conducted on penaeid shrimp (Holt et al. 2020), finfish (Mouchet et al. 2012;Ghanbari et al. 2015), and molluscs (King et al. 2012;Dubé et al. 2019). Similar studies have provided an insight into threatened species too (West et al. 2019), enhancing conservation efforts. The environment can have a significant modifying effect on the composition of host microbiomes (Lokmer and Mathias Wegner 2015;Pennisi 2019;Scanes et al. 2021a, b). This can prove detrimental for host health, particularly when human influence (e.g., habitat degradation or pollution) negatively contributes to dysbiosis or a pathobiome state (reviewed in Bass et al. 2019).
Considering the limited information of toheroa health, in this study the bacterial composition in toheroa is explored in specimens across their entire geographic range for the first time (Fig. 1). The toheroa bacterial composition is described (key taxonomic constituents), and an investigation is carried out into how locally different environmental conditions (whether freshwater outflow is present or not) might modulate the bacterial composition in toheroa tissues.

Specimen collection
Sampling toheroa populations is challenging given their protected status and extreme scarcity in some regions (Foxton and Te-Oneroa-a-Tōhē). This is reflected in the small sample sizes used in the present study. Toheroa were sampled from Te-Oneroa-a-Tōhē (Ninety Mile Beach) in the far north (November 2019; n = 6), Ripiro Beach (Kaipara; September 2019; n = 15), Foxton Beach in the lower North Island (March 2019; n = 1), and Oreti Beach in the far south (February 2019; n = 6) ( Fig. 1 and Table 1). At Ripiro Beach, five toheroa were collected from each of three sites (Island, Mahuta Gap, and Kopawai). Dissections were carried out on return to the laboratory (Sulphur Point, Tauranga). Sampled toheroa were sectioned from anterior to posterior (encompassing all major organs). Tissues were placed in 70% EtOH for DNA extraction. Toheroa are protected, and their collection is restricted under Fisheries Regulations 2013 (SR 2013/482 r25). Specimens were subsequently gathered under special permit issued by Fisheries New Zealand (Special Permit no.706-2).

DNA extraction and PCR amplification
Tissue from toheroa gills (20 mg) and digestive glands (5 mg) was excised for DNA extraction. Gill and digestive gland tissues (25 mg combined) were macerated using a scalpel and transferred into 1.5 mL microcentrifuge tubes. Tissue digest buffer (180 µL) and digest enzyme (20 µl proteinase K) were added to each sample and lysed overnight at 56 °C on a heated-shaking platform. DNA extraction was carried out using the QIAamp ® DNA Mini Kit (Qiagen, CA, USA), following manufacturer's instructions. The presence of amplifiable DNA within the extract was confirmed using an 18S rRNA internal control (Life Technologies, Ribosomal 18S rRNA Endogenous Control). Extracted DNA was quantified using the Qubit fluorometer high sensitivity double stranded DNA protocol (Life Technologies, Auckland, NZ). DNA was stored at − 20 °C until further analyses.
Polymerase chain reaction (PCR) was used to amplify the V4-V5 regions of the 16S rRNA gene from extracted DNA using the primer set F515 (5′CCA TCT CAT CCC TGC GTG TCT CCG ACT CAG -unique IonCode barcode-GAT , and 2 µL of genomic DNA template (5 ng µL −1 ). A lower DNA template concentration (1 ng µL −1 ) was used for three samples to account for PCR inhibition. Reactions were performed in triplicate for each sample with the following thermocycler conditions: 3 min initial denaturation at 94 °C, followed by 30 cycles of 45 s at 94 °C, 1 min at 50 °C, and 1.5 min at 72 °C, with a final 10 min elongation at 72 °C. A positive and a negative control was run in every PCR. Triplicate reactions were pooled, and the PCR products were verified on a 1% SYBR Safe stained agarose gel. Wells were loaded with 5 µL PCR product and 2 µL 3 × loading dye (Invitrogen). Gels were run for 25 min at 70 V with a 1 KB + ladder and visualised using an Alpha Imager. Amplicons were stored at − 20 °C until sequencing.
PCR amplicon concentration for each sample was normalised and purified using a SequelPrep Normalization Kit, following the manufacturer's instructions (Life Technologies, Auckland). DNA sequencing was undertaken at the University of Waikato DNA Sequencing Facility using an Ion-Torrent Personal Genome Machine (PGM) DNA sequencer with an Ion 318v2 chip (Life Technologies). Raw sequences were filtered in mothur to remove long and short reads and reads with excessive homopolymers (Schloss et al. 2009). A total of 1,519,764 sequences (250 bp) were quality filtered using USEARCH (Edgar 2010). High-quality reads (274,069) were then clustered into operational taxonomic units (OTUs) at a similarity threshold of 97%. OTUs (FASTA) were classified using SILVA Incremental Aligner (SINA) (v. 1.2.11) (Pruesse et al. 2007;Quast et al. 2013). The minimum query identity was 80% and taxonomic placements at < 70% identity match were assigned 'unknown'. Taxonomic classifications were obtained at kingdom, phylum, class, order, family, and genus level. Sequences were filtered for chloroplast, mitochondria, Archaea, Eukaryota, singletons and 'unknown', leaving a total of 553 OTUs, with an average of 9990 reads per sample (n = 24). Four samples were removed from any further analyses due to considerably lower reads (< 2000). A 'phyloseq' object was created using taxonomic data, OTU count table, and metadata (McMurdie and Holmes 2013). Several OTUs were compared to published sequences using the NCBI nucleotide BLAST tool (Johnson et al. 2008). Statistical analyses of bacterial composition data were carried out in R (R Core Team 2013).

Data analysis and visualisation
Prior to data analysis, OTU counts were normalised using scaling with ranked subsampling (SRS) (Beule and Karlovsky 2020). Following taxonomic assignment to 16S rRNA sequences, taxonomy data, OTU counts, and metadata were explored using the 'phyloseq' and 'microbiome' packages (McMurdie and Holmes 2013; Lahti et al. 2017). Using the same packages in R, stacked barplots were created based on the relative abundance of the top 10 taxa at the genus level; all other taxa were merged as 'Other'. The Bray-Curtis dissimilarity index was used to assess dissimilarity between the bacterial communities in toheroa tissues sampled from different populations (locally and geographically). Normalised (SRS) OTU counts were log 10 (x + 1) transformed prior to dissimilarity index use. Hierarchical clustering was used to visualise dissimilarity between samples based on the Bray-Curtis index. To establish whether dissimilarities of the bacterial composition in toheroa was statistically Several other statistical and graphical packages were used to produce figures and manipulate data in the R statistical environment, including 'ggplot2' (Wickham 2009) and 'ggdendro' (de Vries and Ripley 2020).

Taxonomic composition
This study provides the first detailed description of the bacterial composition in toheroa tissues, which we have demonstrated to vary both between beaches and between beds at a single beach (Figs. 2 and 3). At the family level, Spirochaetaceae (OTUs 1, 2, and 3) was the most highly represented taxa in the toheroa bacterial composition regardless of study site (mean 75%, range 24-97%). Mycoplasmataceae (OTUs 4,9,11,and 23) was the second most abundant taxa (mean 6%, range 1-38%) in Ripiro, Oreti, and Foxton specimens, while Endozoicomonadaceae (OTU5) (Bartz et al. 2018), was second most abundant in Te-Oneroa-a-Tōhē toheroa. Pseudoalteromonadaceae was comparatively moderately abundant in Oreti specimens. Fusobacteriaceae was comparatively highly abundant in the single Foxton specimen (n = 1). Vibrionaceae, while relatively highly abundant in Oreti specimens, was recorded at lower levels at Foxton and Te-Oneroa-a-Tōhē but was not recorded in Ripiro toheroa. At the genus level, unclassified Spirochaetaceae and Spirochaeta-2 were most abundant taxa at sites other than Oreti where Spirochaeta-2, Vibrio, and Mycoplasma were the most abundant taxa. Mycoplasma was present in all specimens and the most abundant taxa in MG.3 from Ripiro (Fig. 2). This same sample also hosted the highest abundance of OTU6, Psychrilyobacter (Fig. 2).
Both Spirochaetaceae and Mycoplasmataceae have previously been reported as dominant constituents in the bacterial compositions of bivalve molluscs (Pierce and Ward 2018;Dubé et al. 2019;Nguyen et al. 2020;Offret et al. 2020). For instance, King et al. (2020) found that Spirochaetaceae bacteria were conserved across tissue types (gill, digestive gland, mantle, and adductor muscle) and sampling environments in Pacific oysters, Crassostrea gigas. This includes P. australis, a congener of the toheroa, where Spirochaetaceae and Mycoplasmataceae were found to be the most abundant taxa present in the siphon and digestive gland, respectively (Biessy et al. 2020). As the present study was based on the composite of gill and digestive gland it is not possible to discern tissue-specific bacterial composition, where different tissues can host very different bacterial compositions (Lokmer et al. 2016b;King et al. 2020). Synechococcus CC9902 was also found to be in the top taxa (in terms of relative abundance) in several specimens from Ripiro and Te-Oneroa-a-Tōhē (Fig. 2). Biessy et al. (2020) identified several possible bacteria that could be associated with tetrodotoxin (TTX) in pipi, including Synechococcus Many studies cited herein examined the V3-V4 variable regions of the 16S rRNA gene (Biessy et al. 2020;Offret et al. 2020;Scanes et al. 2021a), whereas the V4-V5 regions were examined in the present study. It should be noted that primer choice can effect results (Willis et al. 2019) obtained, and, therefore, the comparability of results between studies. Nevertheless, the bacterial composition reported herein does align well with results obtained by Biessy et al. (2020) for a congener of toheroa (P. australis) despite differing methodologies.

Bray-Curtis dissimilarity and indicator OTUs analysis
While between and within-beach differences in the toheroa microbiome were expected due to differences in the environmental conditions between geographically distant sites (e.g., Lokmer et al. 2016a, b;King et al. 2020;Nguyen et al. 2020;Offret et al. 2020), some of the patterns of groupings that emerged from our analyses were unanticipated. Hierarchical clustering placed specimens geographically separated (c. 1240 km) closer together (more similar) than specimens from the same beach (c. 10 km apart) based on their bacterial composition (Fig. 3). For example, the majority of Ripiro samples clustered together, but were more similar to samples from Oreti, the most distant site, than to those collected at either Te-Oneroa-a-Tōhē or Foxton (Fig. 1). The single Foxton sample grouped with Te-Oneroa-a-Tōhē, but also with three of the five samples from Island at Ripiro (Fig. 3).
An examination of the relative abundances of the key indicator OTUs (as identified by 'indicspecies' (De Cáceres et al. 2020; Fig. 4), revealed the identity of the OTUs are driving dissimilarity among hierarchical groupings (Fig. 3). Foremost, Endozoicomonas spp. (OTU5) were largely absent from both the Kopawai and Mahuta Gap (stream associated sites) and relatively abundant elsewhere at sites not associated with consistently flowing freshwater outflow (Table 1). It is, therefore, suggested that these streams are having a modulating effect on the abundance of certain taxa within the bacterial composition in toheroa (see OTUs 5, 37, and 38, Fig. 4). This is probably attributed to differing environmental conditions at these sites driven by the presence or absence of a freshwater stream. As shown by Cope (2018), the presence of a freshwater stream at toheroa habitats on Ripiro Beach (Mahuta Gap and Island) can reduce pore water salinity, lower sub-surface sediment temperature, and modify beach topography.
A PERMANOVA model identified site as a considerable contributor to dissimilarity (F = 3.3, df = 3, p = 0.001). Indicator species (or in this case OTU) analysis was subsequently used to identify key taxa contributing to site-specific dissimilarity. The most indicator OTUs associated to a group or combination of groups was for Kopawai and Mahuta Gap combined (35 taxa). OTU74 (family: Neisseriaceae), OTU70 (family: Desulfocapsaceae), and OTU37 (family: Flavobacteriaceae) were the top three taxa for this combination of groups with an IndVal score > 0.9 (p < 0.001). When Island was added to the group (i.e., all Ripiro Beach sites), 21 taxa were identified as indicator OTUs compared to Te-Oneroa-a-Tōhē. The top two were OTU38 (family: Desulfocapsaceae) (IndVal = 0.961, p = 0.006) and OTU83 (family: Desulfosarcinaceae) (IndVal = 0.952, p = 0.001). Separately, grouping Island and Kopawai, the top indicator OTU identified was OTU99 (family: Saprospiraceae) (IndVal = 0.862, p = 0.028). OTU118 was identified as the top indicator taxa for Te-Oneroa-a-Tōhē (IndVal = 0.707, p = 0.038), taxonomically assigned to Cobetia (family: Halomonadaceae). Finally, for a combination of groups: Island, Kopawai, and Te-Oneroa-a-Tōhē (i.e., excluding Mahuta Gap) OTU5 (Endozoicomonas spp.) was identified as the sole indicator OTU (IndVal = 0.957, p = 0.002). The relative abundances of some key indicator taxa have been used to construct a heatmap (Fig. 4). Despite not being included in indicator OTU analysis, Oreti and Foxton samples were included for comparison of indicator OTUs. Interestingly, when Kopawai, Mahuta Gap, and Te-Oneroa-a-Tōhē were grouped, one indicator OTU, OTU20 (Pseudoalteromonas spp.) was identified. When compared with other samples (Fig. 4), there is clear difference in abundance between sites, but interestingly, the Oreti Beach samples both host the highest relative abundance despite not being included in the indicator OTU analysis (Fig. 4).

Sulfur cycling symbionts
For several taxa identified as indicator OTUs (i.e.,OTUs 99,74,70,37,38,and 83) there is a clear link to Kopawai and Mahuta Gap compared to all other sites (Fig. 4). Of these, three are particularly interesting. OTUs 70, 38, and 83  (De Cáceres et al. 2020). Oreti and Foxton Beach specimens were not included in indicator OTU analysis but are included here for comparison. Relative abundance > 2% is represented by a grey tile for respective samples were taxonomically assigned to the phylum Desulfobacterota (sulfate-reducing bacteria). It is possible that this is attributed to algal blooms that are typical in autumn (austral) on Ripiro Beach (Williams et al. 2013). Dense mats of decaying phytoplankton cover the mid to low intertidal area of Ripiro Beach during this time (pers. obs., at time of sampling in Sept-19), which likely increases bacterial decomposition following an influx of decaying organic matter. Further, Cope (2018) suggests that the modification of toheroa habitat (topography) by freshwater streams creates an embayment that aggregates phytoplankton to toheroa beds in the mid to low intertidal (Ross et al. 2018a). This could explain the higher relative abundance of these taxa (OTUs 70, 38, and 83) in toheroa from stream-bearing sites (Mahuta Gap and Kopawai) compared to toheroa from sites without a freshwater stream. Increased bacterial respiration probably explains the presence of sulfate-reducing bacteria (Desulfobacterota) within the bacterial composition in toheroa as oxygen is depleted within toheroa beds and sulfur (SO 4 2− ) reduction rates increase, e.g., van Erk et al. (2020). It is, therefore, suggested that the detection of Desulfobacterota in the bacterial composition of toheroa tissues (composite: gill and digestive gland) is due to environmental exposure (i.e., Desulfobacterota bacteria being highly abundant in the sediment in toheroa beds) rather than Desulfobacterota being fundamental constituents of the toheroa microbiome. This is supported by the considerably low relative abundance of Desulfobacterota OTUs within the toheroa composite bacterial composition (Fig. 4). Similarly, elsewhere the bacterial composition in the digestive gland of the Manila clam (Ruditapes philippinarum) from Brittany, France, was dominated by taxa of the orders Mycoplasmatales and Rickettsiales, not reflective of the surrounding sediment that was dominated by Desulfobacterales (36%) (Offret et al. 2020). Thus, the influence of environment on the bacterial composition of shellfish hosts can be significant. For example, elevated pCO 2 (Scanes et al. 2021a(Scanes et al. , 2021b and Vibrio spp. infection and elevated temperature (Lokmer and Mathias Wegner 2015;Green et al. 2019) have been shown to considerably modify the bacterial composition of oyster spp. hosts.
Interestingly, links to the sulfur cycle were made elsewhere too. When the sequence of OTU3 (average relative abundance: 20.2%) was compared to published sequences in the NCBI database (BLAST) a moderate identity was found to a spirochaete endosymbiont (88.8%, AJ620512.1) from a gut-less worm that relies on sulfur cycling symbionts for nutrition (Blazejak et al. 2005;Cano et al. 2020). Additionally, some Endozoicomonas spp. (OTU5) are associated with sulfur cycling roles (endosymbiosis) in coral host tissues (Tandon et al. 2020). Both Spirochaetes and Endozoicomonas taxa together have been shown to be dominant symbionts of other marine organisms i.e., red coral Corallium rubrum (van de Water et al. 2016). Given the multiple links between the toheroa-associated bacterial composition and well-known symbionts (sulfur cycling) future work in this area should consider their potential importance for toheroa health and potentially nutrition.

Closing remarks
Compositional variation in the bacterial taxa present in toheroa appears to be driven by both geography and habitat, with the presence of freshwater outflows explaining differences between beds across a single beach. This could be due to the modification of environmental conditions (e.g., salinity, pH, temperature, and nutrient input), linked to freshwater outflow from terrestrial sources. The decline of toheroa populations across the North Island (Te Ika-a-Māui) has previously been linked to decreases in the flow of freshwater to the toheroa beaches resulting from forestry and water extraction. Williams et al. (2013) reported 40 and 64% reductions in the number of streams reaching Ripiro and Te-Oneroa-a-Tōhē, respectively, since 1950. Considering the apparent influence of these freshwater outflows on the bacterial composition in toheroa, a reduction in freshwater quantity or quality could have significant implications for toheroa health and population recovery.
Compositional differences found herein could be evidence of dysbiosis in toheroa at sites where toheroa habitat requirements are not being fully met. On the other hand, these differences might represent adaptation to different habitat types (sites with and without streams), with considerably different environmental conditions (Cope 2018). A better understanding of the mechanisms driving these apparent habitat-related differences of the bacterial composition in toheroa tissues, and knock-on impacts to physiology (homeostasis and resistance) and performance (survival and growth), could provide information crucial for future restoration efforts.