Microbial community analysis in the gills of abalones suggested possible dominance of epsilonproteobacterium in Haliotis gigantea

Sample preparation and DNA extraction

Five cultured giant abalones, Haliotis gigantea, were collected from the Minami-ise Farming Fishery Center (Minami-ise, Mie, Japan) in October 2016. Wild abalone specimens (H. gigantea, H. discus and H. diversicolor) were obtained from a fish-market in Mie, Japan, from April 2017 to May 2019. The purchased abalones were collected offshore along the coast of Mie prefecture by scuba diver, and were transported in a water tank to the laboratory within 6 h of sampling.

Gill tissues from the cultured H. gigantea (n = 5) specimens collected from the farm and wild H. discus (n = 3) and H. diversicolor (n = 5) specimens collected from the fish-market were pooled into three individual tubes, one for each species (sample code: Hgig1, Hdis1 and Hdiv1). Gill tissues from other H. gigantea and H. discus specimens were prepared individually (sample code: Hgig2-6, Hdis2-4). Gut and foot tissues were also collected from all H. gigantea specimens, excluding Hgig1 from which no foot tissues were obtained. For all H. discus and H. diversicolor specimens, only gill tissues were used. The collected gill, gut or foot tissues were then homogenized in sterile artificial seawater using a bead beater homogenizer (4,200 rpm, 30 s; Tietech Co., Nagoya, Japan) followed by a previously described method (Tanaka et al., 2004). Host tissues were removed from the homogenate by a quick centrifugation step (1 s, 8,000 g), and the supernatant was transferred to new tubes and centrifuged for 20 min at 15,000 g to recover bacterial cells.

Seawater and stone samples were collected from the shore in Minami-ise (34.333958 N 136.692204 E) within 30 m from the Minami-ise Farming Fishery Center. Fifty milliliters of seawaters were filtered by passing through 0.22 µm filter paper and resuspended in sterile phosphate-buffered saline (PBS: 130 mM NaCl, 10 mM Na₂HPO₄/NaH₂PO₄; pH 7.4) (sample code: SW). Stones were shaken in sterile PBS (sample code: ST) and bacterial cells were recovered from the SW and ST PBS samples by centrifuging for 20 min at 15,000 g. Bacterial genomic DNA from each bacterial pellet was extracted using a Promega DNA purification system (Promega Corp., Madison, WI, USA) according to the manufacturer’s instructions.

PCR amplification, illumina MiSeq sequencing and sequence processing

All gill, SW and ST samples were used for the analysis of microbial communities. The V1–V2 region of 16S rRNA genes was amplified using Ex Taq (TaKaRa Biotechnology Corp., Kyoto, Japan). The first PCR step was performed using primers 27F-mod (5′-ACACTCTTTCCCTACACGACGCTCTTCCGATCTAGRGTTTGATYMTGGCTCAG-3′) and 338R (5′-GTGACTGGAGTTCAGACGTGTGCTCTTCCGATCTTGCTGCCTCCCGTAGGAGT-3′). The amplification conditions were as follows: initial denaturation of 2 min at 94 °C followed by 24 cycles of denaturation for 30 s at 94 °C, primer annealing for 30 s at 55 °C, and primer extension at 72 °C for 30 s. The amplified PCR products were purified using a Wizard SV Gel and PCR Clean-Up System (Promega Corp., Madison, WI, USA) and used for the second PCR step, which was performed using primers with a tag sequence. After the second PCR, the products were once again purified using a Wizard SV Gel and PCR Clean-Up System before sequencing on a MiSeq platform (Illumina Inc., San Diago, CA, USA).

The raw paired-end FASTQ reads were demultiplexed using the Fastq barcode splitter (http://hannonlab.cshl.edu/fastx_toolkit/index.html) and imported into the Quantitative Insights Into Microbial Ecology 2 program (QIIME2, ver. 2019.7, https://qiime2.org/). Demultiplexed reads were quality filtered, denoised, chimera checked and dereplicated using a DADA2 denoise-paired plugin (Callahan et al., 2016). To be equal to sampling-depth, sequences were rarefied at 23,000 reads using qiime feature-table rarefy (Weiss et al., 2017). For preparing Upset diagram, sequences from H. gigantea and H. discus were also combined into one for each host species using feature-table merge plugin on Qiime2 in addition to prepare individual samples, and then the sequences were resampled at 23,000 reads using qiime feature-table rarefy plugin on Qiime2 (Weiss et al., 2017). Next, the align-to-tree-mafft-fasttree pipeline from the q2-phylogeny plugin was used to perform multiple sequence alignment, remove phylogenetically uninformative or ambiguously aligned sequences, and to generate unrooted and rooted phylogenetic trees (Lane, 1991; Price, Dehal & Arkin, 2010; Katoh & Standley, 2013). Diversity metrics were calculated using the core-metrics-phylogenetic pipeline from the diversity plugin on QIIME2. For alpha diversity, Shannon index (Shannon, 1948), observed OTU number (DeSantis et al., 2006), Chao 1 index (Chao, 1984) and Good’s coverage (Good, 1953) were calculated using the diversity alpha command on QIIME2. For PCoA plots, unweighted UniFrac data and Shannon diversity values were calculated by the core-metrics-phylogenetic pipeline on Qiime2, and they were combined and visualized using the R packages “qiime2R”, “tidyverse”, “phyloseq” (McMurdie & Holmes, 2013) and “ggplot2” (Wickham, 2009) in R (v3.5.0). Taxonomic assignments were performed using the qiime feature-classifier classify-sklearn on Greengenes v_13.8 (McDonald et al., 2012). Taxa bar plots were constructed using the plugin qiime taxa bar plot. For analysis of core or unique microbiome (Reich et al., 2018; Neu, Allen & Roy, 2019), Upset diagram was generated using R package UpSetR (version 1.3.3) with grouped next-generation sequencing (NGS) sequence data by each host species or environment (Lex et al., 2014). All of the data were deposited at the Sequence Read Archive (SRA) under the accession number PRJDB8953.

Obtaining the 16S rRNA gene sequences of the uncultured epsilonproteobacterium by cloning

PCR amplification of the bacterial 16S rRNA gene was performed using primer 27F (5′-AGAGTTTGATCCTGGCTCAG-3′, Lane, 1991) and primer 1492R (5′-GGTTACCTTGTTACGACTT-3′, Lane, 1991). The bacterial 16S rRNA gene clone library was constructed using PCR amplicons obtained from Hgig1. PCR reaction mixtures contained 1 × PCR reaction buffer, 200 µM dNTP, 5 pmol of each primer, 2.5 units Ex Taq polymerase (TaKaRa Biotechnology Corp., Kyoto, Japan), and 10–100 ng of community DNA in a total volume of 50 µl. PCR reactions were performed on a thermal cycler (iCycler; Bio-Rad Laboratories, Hercules, CA, USA) using the following amplification conditions: initial denaturation of 4 min at 95 °C followed by 25 cycles of denaturation for 30 s at 95 °C, primer annealing for 30 s at 55 °C, and primer extension at 72 °C for 1.5 min. This was followed by a final extension reaction at 72 °C for 7 min. PCR product was ligated into the TOPO TA cloning vector (Invitrogen Corp., Carlsbad, CA, USA) according to the manufacturers’ instructions. Ligation products were transformed into Escherichia coli One Shot TOP10 cells (Invitrogen Corp., Carlsbad, CA, USA) and clones were amplified by PCR using vector-specific primers. Plasmid DNA with insertions was sequenced using primer 27F. Partial sequencing was performed using the BigDye terminator cycle sequencing method, and an ABI 3730 sequencer (Applied Biosystems, Foster City, CA, USA). The resulting chromatograms of DNA sequences were examined using Chromas 2.33. Homology searches were performed using sequences and close relatives were identified using a BLAST search of the GenBank database on the National Center for Biotechnology Information website (http://www.ncbi.nlm.nih.gov/). The sequences identified as Epsilonproteobacteria by BLAST search were selected from random clones generated during cloning. The phylogenetic characterization of the selected sequence was estimated using the phylogenetic tree with representative sequence of Epsilonproteobacteria by the Maximum Likelihood method in MEGA 7.0 software (Kumar, Stecher & Tamura, 2016).

Nested PCR analysis and sequencing of the uncultured epsilonproteobacterium

The nested PCR assay targeted 356 bp of the uncultured epsilonproteobacterial 16S rRNA gene. At the first PCR, bacterial universal primers 27F and 1492R were used, and then the uncultured epsilonproteobacterium-specific primers, Eps222F (5′-CGCTAAGAGATTGGACTATAT-3′) and Eps578R (5′-GACTTAATAGGACACCTACATACC-3′) (designed in this study) were used at the second PCR. DNA samples extracted from gill, gut or foot tissue of each H. gigantea specimen were used as DNA templates. The first-round master mixture contained the following: primers (27F and 1492R, Lane, 1991) at 0.2 µM (each), 25 µl volume of EmeraldAmp^® PCR Master Mix (TaKaRa Biotechnology Corp., Kyoto, Japan), approximately 50 ng of DNA, and up to 50 µl with distilled H₂O. Distilled H₂O was used as the template for negative controls. Cycle parameters were 98 °C for 10 s; 10 cycles of 98 °C for 10 s, 55 °C for 30 s and 72 °C for 90 s; and 72 °C for 5 min. PCR reactions were performed using an iCycler (Bio-Rad Laboratories, Hercules, CA, USA). A 2.5-µl aliquot of the amplified PCR product was transferred to a new master mixture containing primers (Eps222F and Eps578R) at 0.2 µM (each) in a 25 µl volume of EmeraldAmp^® PCR Master Mix (TaKaRa Biotechnology Corp., Kyoto, Japan) made up to a volume of 50 µl with distilled H₂O. Cycle parameters were 98 °C for 10 s; 25 cycles of 98 °C for 10 s, 57 °C for 30 s and 72 °C for 60 s; and 72 °C for 5 min. A 5-μl aliquot of the amplified PCR product was analyzed by agarose (1.5% (wt/vol) prepared in TAE buffer) gel electrophoresis and Midori Green Direct (Nippon Genetic Europe GmbH, Düren, Germany) staining, and the gel was photographed under UV light. The presence of a band at around 350 bp was interpreted as a positive result. Sequencing and phylogenetic analysis were performed as previously described in the subsection “Obtaining the 16S rRNA gene sequences of the uncultured epsilonproteobacterium using cloning method” above.

Fluorescence in situ hybridization localization of the uncultured epsilonproteobacterium

Fluorescence in situ hybridization (FISH) was performed on the Hgig4 sample using a previously described method (Tanaka et al., 2016). Gill tissue was obtained from Hgig4 and placed in 75% ethanol at −30 °C before being fixed in 4% (v/v) paraformaldehyde/PBS at 4 °C overnight. Following fixation, specimens were rinsed and dehydrated in a 50%, 70%, 80%, 85%, 90%, 95% and 100% ethanol series, followed by 100% xylene. The fixed gill specimens were then embedded in paraffin and sliced into 10-μm transverse sections using a microtome (RV240, Yamato, Japan), before being placed on APS-coated microscope slides (Matsunami, Japan), and stored in slide boxes at room temperature until deparaffinization. Wax was removed, and tissue was rehydrated in a decreasing ethanol series. The sections were dewaxed and hybridized at 47 °C for 3 h using a solution containing 20% formamide, 0.9 M NaCl, 20 mM Tris–HCl (pH 7.4), 0.1% SDS and 0.5 μM FITC-labeled Eub338 probe (5′-GCTGCCTCCCGTAGGAGT-3′, Amann, Krumholz & Stahl, 1990) or 0.5 μM TAMRA-labeled Eps222 probe (5′-CGCTAAGAGATTGGACTATAT-3′), which was designed to specifically target the 16S rRNA gene of the uncultured epsilonproteobacterium (designed in this study). After hybridization, the sections were washed twice with washing buffer containing 20 mM Tris-HCl (pH 7.4), 180 mM NaCl and 0.01% SDS for 30 min at 48 °C, and the sections were then rinsed with ddH₂O and air-dried. An epifluorescence light microscope (Eclipse 400; Nikon Instech., Tokyo, Japan), was used to observe the stained cells.

Microbial community analysis by NGS

The microbial communities on the gills of various marine invertebrates were characterized by sequencing the V1–V2 region of the 16S rRNA gene. A total of 473,349 quality-filtered sequence reads were obtained from 13 samples (PRJDB8953). Alpha diversity index values for each sample are shown in supplementary Table 1. For all of the samples in this study, observed OTUs and Chao1 were approximately the same value, and Good’s coverage value was higher than 99.99%, indicating that sequencing depth sufficient for capturing all bacterial species and for downstream analysis. We identified 322 OTUs from 15 phyla, which we taxonomically assigned with the Greengenes database v. 13_8. At the phylum level, the microbial taxonomic composition of most samples, excluding Hgig5, had a high relative abundance of Proteobacteria (67.8% on average). Hgig5 showed almost complete dominance of Spirochaetes (68.4%).

In Proteobacteria, sequence reads were related to Alphaproteobacteria, Gammaproteobacteria or Epsilonproteobacteria (Fig. 1). Alphaproteobacteria was dominant in Hgig1, Hgig5, Hdis2, Hdis4, SW and ST, (12.0–39.5%, Fig. 1). The most abundant Alphaproteobacteria in abalone were aligned with unclassified Rickettsiales (11.2% on average, two OTUs), while an unclassified Rhodobacteraceae (9.1%, one OTU) and Nautella sp. (12.5%, one OTU) were the most abundant Alphaproteobacteria-related sequences in the ST and SW samples, respectively. The relative abundance of Gammaproteobacteria-related sequences was more than 10% in all samples, except Hgig2. The genus Vibrio (Gammaproteobacteria) was observed in all samples. In addition, the relative abundance of Vibrio spp. (six OTUs) was more than 10% in SW and all H. discus specimens, especially in Hdis4 (58.4%, three OTUs). In Hdiv1, an unclassified Endozoicimonaceae (one OTU) accounted for 25.1% of the relative abundance and were more frequently recovered than Vibrio spp. In H. gigantea specimens, Alcanivorax dieselolei (one OTU) showed 11.8% relative abundance in Hgig1, while an unclassified Francisellaceae (one OTU) and an unclassified Gammaproteobacteria (one OTU) were distributed at 9.4% and 5.9% in Hgig6, respectively. Gammaproteobacteria-related sequences in other Hgig samples however were classified into OTUs with relative abundance less than 5%. Sulfur-oxidized Gammaproteobacteria, which was reported as symbiotic bacteria, was not detected in any samples. Epsilonproteobacteria-related sequences were dominantly recovered from SW and all H. gigantea specimens, except Hgig5 (18.3–62.2%, Fig. 1). All Epsilonproteobacteria-related sequences in SW were affiliated with Arcobacter sp. (28.7%, one OTU). Arcobacter sp. was also detected in the others samples, although the relative abundance was less than 5%. Other Epsilonproteobacteria-related sequences were most closely related to unclassified Epsilonproteobacteria (two OTUs). One of these was detected in all H. gigantea specimens and dominated in Hgig1 to Hgig4, accounting for 50.5 ± 0.05% of the relative abundance. The unclassified Epsilonproteobacteria was not found in all Hdis, SW and ST, although the related sequences represented 0.2% of the relative abundance in Hdiv1. Another unclassified Epsilonproteobacteria was detected in Hgig4, Hgig5, Hgig6, Hdis1 and Hdiv1, accounting for up to 7.1% of the relative abundance in Hgig4. Spirochaetes-related sequences were aligned with an unclassified Spirochaetaceae (one OTU). These sequences were observed in all abalone specimens except in Hdis4. The unclassified Spirochaetaceae was abundant in Hgig5, Hdis1 and Hdiv1, accounting for 68.4%, 28.1%, 19.9%, respectively. No Spirochaetes-related sequences were detected in SW and ST. All members of the phylum Tenericutes were affiliated with one OTU in the genus Mycoplasma. Mycoplasma sp. was present in every abalone, and most prevalent in Hdis1 (21.9%), followed by Hdis3 (21.3%). In contrast, the sequences related to Mycoplasma were never observed in SW and ST. Chloroplast-related sequences were only detected in ST and SW, accounting for 22.6% (six OTUs) and 4.5% (six OTUs), respectively. Most Chloroplast-related sequenceare associated with an unclassified Stramenopiles (one OTU) and accounted for up to 12.6% of the relative abundance in ST.

The similarity in the microbiota was supported by principal coordinates analysis (PCoA, Fig. 2). The microbiota in each invertebrate differed from those in the environment, such as from the seawater and stones. The microbial community structure of H. gigantea was similar to one of H. discus, while the microbe of H. diversicolor formed a cluster distinct from other abalone species.

The idea of core and unique microbe was used to clarify the characteristics among samples (Reich et al., 2018; Neu, Allen & Roy, 2019). The largest numbers of unique OTUs was obtained from Hgig with 64 OTUs (Fig. 3). In descending order, the unique OTUs were 52, 49 and 28 OTUs in ST, Hdis and SW, respectively. The number of unique OTUs in Hdiv was the smallest where eight OTUs were obtained. The relative abundance of the unique OTUs showed low value of 3.5%, 2.4%, 3.9% and 7.6% in Hgig, Hdis, Hdiv and SW, respectively (Table 1). On the other hand, the unique OTUs in ST showed high proportion (26.2%), although the highest relative abundance showed only 4.5%. 10 OTUs were obtained from all samples as core microbe, and the relative abundance of the core microbe showed 12.6%, 36.3%, 12.2%, 44.5% and 20.7% in Hgig, Hdis, Hdiv, SW and ST, respectively (Fig. 3; Table 1). Arcobacter sp. and a member of genus Vibrio (one OTU) were shared between all specimen as the core microbe. A total of 14 OTUs were detected as the core microbe in abalones, which was commonly obtained from Hgig, Hdis and Hdiv, and their proportions were 31.3%, 32.6% and 44.0%, in Hgig, Hdis and Hdiv, respectively (Table 1). The core microbe in abalones included the unclassified Spirochaetaceae and Mycoplasma sp. The core microbe in environments, SW and ST, was eight OTUs, including the unclassified Stramenopiles, and their proportions were 2.7% and 25.3% in SW and ST, respectively (Table 1).

Phylogenetic analysis of 16S rRNA gene of the uncultured epsilonproteobacterium

Cloning and sequencing of the long 16S rRNA gene sequence of the uncultured epsilonproteobacterium found in H. gigantea were performed to clarify the phylogenetic position of this bacterium. Four out of 16 cloned sequences completely matched to the sequence of the uncultured epsilonproteobacterium obtained from the 16S rRNA gene amplicon sequencing undertaken in this study. The longest sequence length was 1,436 bp (LC511979).

BLAST analysis of the uncultured epsilonproteobacterial sequence revealed that the sequence could be assigned to various described Epsilonproteobacteria species, such as Arcobacter canalis strain F138-33 (87.43% identity), Arcobacter marinus strain CL-S1 (87.19% identity), Sulfurovum lithotrophicum strain 42BKT (87.01% identity) and Helicobacter pullorum strain ATCC51801 (86.78% identity), but these identity scores were low. Among the environmental clones in public databases, the uncultured epsilonproteobacterium specimen was most closely related to an uncultured clone sequence found in the gill tissue of a wood-feeding gastropod in genus Pectinodonta (96.12–96.78% identity, Zbinden et al., 2010), followed by an uncultured clone sequence from gill tissue of thyasirid bivalves from a cold seep in the Eastern Mediterranean (94.86% identity, Brissac et al., 2011). In the phylogenetic analysis with Maximum Likelihood method, the uncultured epsilonproteobacterium grouped within the order Campylobacterales (Fig. 4). However, the uncultured epsilonproteobacterium was not able to be classified below the order, since the phylogenetic position of the sequence was at the root of the family Helicobacteraceae.

Detection of the uncultured epsilonproteobacterium from each part of H. gigantea

Results of the 16S rRNA gene amplicon sequencing analysis revealed that the uncultured epsilonproteobacterium was dominant only in the gill tissues of H. gigantea. Therefore, in order to determine whether the uncultured epsilonproteobacterium was restricted to the gill tissues, a PCR assay of gill, gut and foot tissues from H. gigantea was performed using uncultured epsilonproteobacterium-specific primers. PCR products were obtained from gill tissue from Hgig2, Hgig3 and Hgig6 (Table 2). Amplified products were also obtained from the gill tissue from Hgig1, Hgig4, gut tissues from Hgig3 and foot tissues from Hgig4, but band visibility was low. All amplification product sequences matched that of the uncultured epsilonproteobacterial sequence obtained by cloning.

Localization of the uncultured epsilonproteobacterium on gill tissue of H. gigantea

In semi-thin section FISH analysis, Eps222-positive bacteria were observed on the gill surface of Hgig4 (Fig. 5) where they formed microcolonies on gill tissue (Figs. 5B and 5C). Eps222-positive bacteria showed a coccus-like morphology and measured ~1.0 µm in diameter (Fig. 5C).