Skip to main content
Download PDF

Targeted single-cell genomics reveals novel host adaptation strategies of the symbiotic bacteria Endozoicomonas in Acropora tenuis coral

Microbiome volume 10, Article number: 220 (2022) Cite this article

Abstract

Background

Endozoicomonas bacteria symbiosis with various marine organisms is hypothesized as a potential indicator of health in corals. Although many amplicon analyses using 16S rRNA gene have suggested the diversity of Endozoicomonas species, genome analysis has been limited due to contamination of host-derived sequences and difficulties in culture and metagenomic analysis. Therefore, the evolutionary and functional potential of individual Endozoicomonas species symbiotic with the same coral species remains unresolved.

Results

In this study, we applied a novel single-cell genomics technique using droplet microfluidics to obtain single-cell amplified genomes (SAGs) for uncultured coral-associated Endozoicomonas spp. We obtained seven novel Endozoicomonas genomes and quantitative bacterial composition from Acropora tenuis corals at four sites in Japan. Our quantitative 16S rRNA gene and comparative genomic analysis revealed that these Endozoicomonas spp. belong to different lineages (Clade A and Clade B), with widely varying abundance among individual corals. Furthermore, each Endozoicomonas species possessed various eukaryotic-like genes in clade-specific genes. It was suggested that these eukaryotic-like genes might have a potential ability of different functions in each clade, such as infection of the host coral or suppression of host immune pathways. These Endozoicomonas species may have adopted different host adaptation strategies despite living symbiotically on the same coral.

Conclusions

This study suggests that coral-associated Endozoicomonas spp. on the same species of coral have different evolutional strategies and functional potentials in each species and emphasizes the need to analyze the genome of each uncultured strain in future coral-Endozoicomonas relationships studies.

Video Abstract

Background

In recent years, coral bleaching has been known to have several impacts on marine ecosystems. Corals have symbionts, including intracellular symbiotic dinoflagellates, coral-associated bacteria, and viruses, which have been reported to play an important role in coral health and survival [1, 2]. The importance of consideration for the function of the coral itself as well as the interactions among the holobionts in the coral ecology has been globally recognized. In the functional analysis of coral symbionts, symbiotic dinoflagellates have been the main target, but recently, the functions of coral-associated bacteria have also been gradually clarified, including the degradation of dead symbiotic dinoflagellates [3], provision of nitrogen sources [4], and protection from pathogenic bacteria via antibiotic activities [5].

The genus Endozoicomonas has been found to be abundant in corals worldwide and is associated with an extensive range of marine organisms, including shellfish [6, 7], sea slugs [8], sponges [9, 10], and sea anemones [11, 12]. Endozoicomonas has been the most widely studied coral-associated bacterium, and it has been shown that phylogenetically different Endozoicomonas species colocalize within a single individual [3, 13]. In addition, Endozoicomonas has been reported to play several functional roles in coral survival, such as suppressing mitochondrial degradation, supplying sugars, and contributing to the coral sulfur cycle through dimethylsulphoniopropionate (DMSP) [14, 15]. However, little is known about the mechanisms by which Endozoicomonas establishes symbiotic relationships with its host corals. Recently, Ding et al. investigated the host cell entry mechanism of Endozoicomonas montiporae in Montipora coral. They hypothesized that this bacterium used eukaryotic-like genes, such as ephrin ligands, to enter the host cell [14]. This adaptation strategy is consistent with the fact that Endozoicomonas is not an obligate symbiont [14, 16]. Endozoicomonas diversity has been studied at the 16S rRNA gene level, but not yet at the genome level. Therefore, there are method-dependent gaps in our understanding of Endozoicomonas diversity. Furthermore, comparative genome analysis of Endozoicomonas spp. from the same host coral species can reveal the evolutionary and adaptive strategies of each Endozoicomonas. However, these understandings require genomic information on many Endozoicomonas from the same host coral species.

The lack of genomic information is due to the difficulty in culturing host-associated bacteria [16]. Thus, culture-independent approaches are needed for studying genome data to understand the relationship between Endozoicomonas and its host. Currently, shotgun metagenomics with binning is the mainstream method for microbial genome sequence analysis, but this method is greatly influenced by microbiome composition and diversity [17].In most cases, especially in genome analysis of closely related species with highly similar 16S rRNA gene sequences, it is impossible to identify individual strain-level draft genome sequences by metagenomic binning [18]. Furthermore, detailed genome sequence analysis of host-associated bacteria is complicated because contamination by host also significantly affects binning efficiency [16, 19,20,21]. As an alternative to metagenomics, single-cell genomics is very useful, as it can reveal the genome diversity of individual cells [22,23,24] from a mixed cell population. On the other hand, in the current single-cell genomics, single-cells are randomly sorted for subsequent analysis. Therefore, when host-derived sequences are contaminated as in this case, the efficiency of acquiring target microbe is greatly reduced. Deep sequencing was required to obtain multiple genome sequences to reveal genome diversity, which required wasteful cost and effort.

Here, we suggest a novel approach for targeted single-cell genomics by detecting the target gene using droplet microfluidics. In this study, we identified seven draft genomes of Endozoicomonas spp. from A. tenuis sampled at four different sites in Okinawa Prefecture, Japan. We then investigated the genomic differences between different clades of Endozoicomonas symbiotic with A. tenuis coral and examined the role of eukaryotic-like genes. Finally, we proposed a symbiotic strategy model for Endozoicomonas bacteria in A. tenuis corals, which revealed a diverse adaptation process.

The results of single-cell genome analysis indicate the existence of at least two strategies for microbe-host interactions, indicating the possibility of microbial regulation of host metabolic systems. These findings may provide insights into the robustness of coral reefs.

Materials and methods

Coral branch sampling

A. tenuis samples were collected at four sites, namely Sesoko Minami (26° 53.05′ N, 127° 85.77′ E), Ishikawabaru (26° 67.51′ N, 127° 87.00′ E), Onna Village (26° 50.36′ N, 127°85.39′ E ), and the coral aquaculture center of Yomitan Village (26° 40.90′ N, 127° 71.55′ E) (Fig. 1A). Official permission to collect coral samples was obtained from Okinawa Prefecture in accordance with Okinawa Prefecture Fishing Regulation. The level of coral bleaching was evaluated on a 3-point scale (0: no bleaching, 1: slight bleaching, 2: complete bleaching). Small pieces of branches (~3 cm) were collected from A. tenuis corals and kept in a sterile zipper bag (140×100×0.04 mm; SEISANNIPPONSHA Ltd., Japan) filled with seawater from each sampling site. For 16S rRNA gene amplicon sequencing, coral branches were immediately transferred into a 5-mL tube (Eppendorf Ltd., Germany) containing 4 mL of RNAlater (Thermo Fisher Scientific, Waltham, MA, USA) and then frozen with liquid nitrogen. The frozen samples were kept in a −80°C deep freezer until analysis. For single-cell genome sequencing, the branches were kept on ice until just before the subsequent process.

Fig. 1
figure 1

Distribution and diversity of the Endozoicomonas genus in A. tenuis from four sites in Japan. A Coral samples were collected at four sites in Okinawa, Japan. Red: Ishikawabaru, green: Sesoko Minami, blue: Onna, purple: Yomitan. B 16S rRNA gene composition. Each color indicates Endozoicomonas and is shown for amplicon sequence variants that are above 8% in any sample. C Heatmap plot based on 16S rRNA gene copy number denoting Endozoicomonas abundance. The copy number of 16S rRNA gene is normalized by the surface area of the coral. The uncultured Gammaproteobacteria clade (black) is displayed as the outer group. The star shows corals that did not bleach in August 2016. Bleaching was recorded in Ishikawabaru and Sesoko Minami, but not in Yomitan and Onna

Full size image

DNA extraction and library preparation

The frozen coral branch was thawed on ice. The coral branch was then picked with tweezers and washed with artificial seawater. Then, in a sterilized zipper bag, the coral tissue was blown off using a Waterpik (EW-DJ61-W, Panasonic Corp., Japan), and the suspensions were collected into a 50-mL tube (Nippon Genetics Co., Ltd., Tokyo, Japan). After centrifugation at 10,000×g, 4°C for 30 min, the supernatant was removed, and the pellet was resuspended in 500 μL of artificial seawater. The suspension was transferred into a 1.5-mL tube (SSIBio, Lodi, CA, USA) and centrifuged at 15,000×g, 4°C for 15 min. After the supernatant was removed, the pellet was stored at −80°C until DNA extraction. DNeasy Plant Mini Kit (69104; QIAGEN, Hilden, Germany) was used for DNA extraction according to the manufacturer’s instruction, with a small modification. The thawed pellet was mixed with 400 μL of Buffer AP1 and transferred into 2-mL screw cap tubes with stabilized 0.1 mm zirconia beads (Yasui Kikai Corp., Osaka, Japan). The tubes were homogenized three times at 2500 rpm for 60 s at 60-s intervals. After centrifugation, 4 μL of RNase A was added, and the sample was incubated at 65°C for 10 min before DNA extraction was performed following the kit protocol.

16S rRNA gene sequencing and analysis

The bacterial community was identified by targeting the variable region V1-V2 of the 16S rRNA gene using the primer set 27F and 338R (27F: 5′-GAG TTT GAT CCT GGC TCA G-3′, 338R: 5′-GCT GCC TCC CGT AGG AGT-3′). PCR amplification of the 16S rRNA gene amplicon (300 bp) was performed in a 25-μL mixture, and the amplicons were sequenced using Ion Torrent PGM with 318 Chip v2 (Thermo Fisher Scientific).

16S rRNA gene analysis was performed using QIIME2 2020.2 [25]. For quality control, amplicon sequence variants (ASVs) were constructed using DADA2 (via q2-dada2) [26]. All ASVs were aligned using MAFFT (via q2-alignment) [27]. Phylogenetic tree was constructed based on masked aligned ASVs using FastTree2 (via q2-mask and q2-phylogeny) [28]. Taxonomic assignment was performed using the q2-feature-classifier [29] classify-sklearn naïve Bayes taxonomy classifier against the Silva 138 99% OTUs full-length sequences [30]. The downstream analysis was performed using the R bioconductor phyloseq [31] and ggtree [32] in R version 3.5.2.

Droplet digital PCR for quantification of 16S rRNA gene copy number

The PCR mixture contained 10 μL of 2× ddPCR Supermix for Probes (Bio-Rad, Hercules, CA, USA), 1.0 μL of 10 μM forward and reverse primers, 1.0 μL of 10 μM Taqman probe, 1.0 μL of 1 mM Dextran fluorescein (FD10S; Merck KGaA, Darmstadt, Germany), 1.0 μL of extracted DNA, and 5.0 μL of nuclease-free water (Sigma-Aldrich, St. Louis, MO, USA). A total of 20 μL of the PCR mixture was introduced into microfluidic devices (Supplemental Figure 1), and 40 μm of droplets was generated with droplet generation oil for probes (Bio-Rad). The primers and Taqman probes were synthesized by Integrated DNA Technologies (IDT, Inc., Coralville, IA). The sequence of the forward primer, reverse primer, and Taqman probe was 5′-AGA GTT TGA TCM TGG CTC AG-3′, 5′-GCT GCC TCC CGT AGG AGT-3′, and [6-FAM]-5’-CAG GCC TAA-[ZEN]-CAC ATG CAA GTC-[IBFQ]-3′, respectively. The droplets were collected into 0.2-mL tube and then thermal cycled using a T100™ Thermal Cycler (Bio-Rad) with the following protocol: 95°C for 10 min, followed by 40 cycles of 94°C for 30 s and 60°C for 1 min, 95°C for 10 min, and maintenance at 4°C. Subsequently, the droplets were reintroduced into microfluidic devices for digital counting. We have developed a system for calculating the rate of fluorescence-positive droplets using a 488-nm laser (Lambda mini, Tokyo Instruments, Tokyo, Japan) (Supplementary Fig. S1). Fluorescent signals from the droplets were detected by a photosensor module (H11902-01; Hamamatsu Photonics KK, Hamamatsu, Japan) and recorded by an oscilloscope (PICOSCOPE 2204A; Pico Technology Ltd., St. Neots, Cambridgeshire, UK). Combined with the information of droplet diameter measured with ImageJ, 16S rRNA gene copy number per 1 ng of extracted DNA was calculated [33].

Measurement of surface area of coral skeletons by CT scan

The coral branches with surface tissue removed by water picking were used for calculating the surface area. 3D images of the coral branches were obtained using an X-Ray CT scanner (TDM1300-IS; Yamato Scientific, Tokyo, Japan). The surface area was calculated by the Mimics software (Materialise, Leuven, Belgium). The 16S rRNA gene copy number per 1 cm2 of surface area was calculated.

Preparation of bacterial suspension for single-cell genome sequencing

A small piece of coral branch was collected and kept in a 25-mL tube (Iwaki, Tokyo, Japan) filled with seawater from each sampling site. The tube was kept on ice immediately. Each branch was transferred into 5 mL of 0.22-μL-filtered UV-treated seawater of each sampling site and crushed thoroughly by a scalpel (Feather Safety Razor Co. Ltd., Japan). After 5 min of standing on ice, the supernatant was collected into three 1.5-mL tubes. The tube was subsequently centrifuged at 300×g for 5 min, and the supernatant was transferred into a new 1.5-mL tube to remove symbiotic dinoflagellates and other larger particles. The collected supernatant was centrifuged at 8000×g for 5 min. The supernatant was removed, and 500 μL of 0.22-μL-filtered UV-treated seawater was added to resuspend the pellet. The centrifugation at 300×g was repeated until most of the symbiotic dinoflagellates were removed. The centrifugation at 8000×g and the washing steps were repeated three times. Next, the pellet was resuspended in 500 μL of 1x SYBR Green (Thermo Fisher Scientific, Hampton, NH) for 5 min for DNA staining. After centrifugation at 8000×g, the pellet was resuspended in 50 μL of 0.22-μL-filtered UV-treated seawater, and the cell concentration was calculated with a bacterial counter (SLGC, Tokyo, Japan) under a fluorescent microscope (CKX53; Olympus Optical Co Ltd., Tokyo, Japan).

Whole genome amplification and detection of Endozoicomonas in gel beads

Cells were encapsulated into 40-μm microfluidic droplets at a concentration of 0.3 cell/droplet using a microfabricated microfluidic device reported previously [22]. The theoretical number of encapsulated cells in a droplet was derived by the Poisson distribution (Formula 1). The percentage of doublets would be less than 3.7% of the total droplets when the cell concentration was 0.3 cell/droplet.

$$p\left(k,\lambda \right)=\frac{\lambda^k{e}^{-\lambda }}{k!}$$
(1)

Next, gel-bead-based whole genome amplification (WGA) was conducted according to the previously reported WGA method (single-cell amplified genomes (SAGs) in the gel are referred to as SAG-gel) [21, 34]. After WGA, the gel beads were stained with DAPI, and amplification of the DNA was confirmed with microscopic observation (CKX53; Olympus). Subsequently, the gel beads were washed and suspended in UV-treated Dulbecco’s phosphate-buffered saline (-) (Thermo Fisher Scientific). To detect the target gene sequence, 20 μL of PCR mixture consisting of 10 μL of PrimeTime Gene Expression Master Mix (Integrated DNA Technologies, Coralville, IA, USA), 1.0 μL of each 10 μM primer (forward and reverse, final concentration 500 nM), 0.5 μL of 10 μM Taqman probe (final concentration 250 nM), and 7.5 μL of gel beads solution was prepared. The sequence of each primer and probe is listed in Supplementary Table 1. On the basis of the 16S rRNA gene amplicon sequencing results, we designed two sets of primers and probes (sets A and B) using the Geneious software (Biomatters Ltd., Auckland, New Zealand). The primer and probe sequences were evaluated with Primer BLAST [35] to confirm that they did not exhibit any specificity to the mitochondrial genome of A. tenuis. To detect the Endozoicomonas sequence amplified in gel beads, either set A or set B was used. The PCR mixture was reintroduced into microfluidic devices, and 50 μm of gel beads-containing microfluidic droplets were generated with carrier oil (HFE-7500 oil with 2% surfactant; RAN Biotechnologies Inc., Beverly, MA, USA). The concentration of the beads was adjusted at 0.3 bead/droplet. The droplets were collected into a 0.2-mL tube, and PCR amplification (95°C for 3 min, 25 cycles of 95°C for 5 s and 50°C for 40 s, and maintenance at 10°C) was conducted. After thermal cycling, the fluorescence of the droplets was observed with a microscope (CKX53; Olympus). The droplets exhibiting both DAPI and FAM fluorescence were manually picked with a micro dispenser (Drummond Science Company Broomall, PA, USA) and transferred into a 0.2-mL tube. Next, second-round WGA was performed using a REPLI-g Single Cell Kit (QIAGEN). Buffer D2 (0.6 μL) was added to each well, followed by incubation at 65°C for 10 min. After that, 8.6 μL of WGA mixture (0.6 μL of stop solution, 1.8 μL of H2O, 5.8 μL of reaction buffer, and 0.4 μL of DNA polymerase) was added, and the mixture was incubated at 30°C for 120 min. The WGA reaction was terminated by heating at 65°C for 3 min. The amplicon yields were quantified by a Qubit dsDNA HS assay kit (Thermo Fisher Scientific). PCR amplification over the 16S rRNA gene was also performed against second-round WGA products. Primer pair sequences for the V3-V4 region were used according to Illumina’s MiSeq system protocols (Forward: 5′-TCG TCG GCA GCG TCA GAT GTG TAT AAG AGA CAG CCT ACG GGN GGC WGC AG-3′, Reverse: 5′-GTC TCG TGG GCT CGG AGA TGT GTA TAA GAG ACA GGA CTA CHV GGG TAT CTA ATC C-3′). PCR amplification was confirmed with agarose electrophoresis (100 V, 15 min), and the amplicon sequences were obtained via Sanger sequencing (Fasmac, Kanagawa, Japan) to determine whether the target Endozoicomonas sequence was acquired.

Endozoicomonas single-cell genome sequencing and genome assembly

Second-round WGA products from single-cell samples were used for next-generation sequencing library preparation with Nextera XT DNA sample prep kit (Illumina, San Diego, CA, USA) according to the manufacturer’s instructions. Each SAG library was sequenced using an Illumina MiSeq 2 × 75bp or 150bp (Supplementary Table 2).

All raw reads were trimmed to remove adapter sequences by fastp v0.20.0 with default parameters [36]. Trimmed reads were assembled using SPAdes v3.12.0 with “-k auto --sc --careful” parameters [37]. All SAGs were concatenated using ccSAG [38]. SAGs with an average similarity of > 99.9% in single-copy marker genes were co-assembled and used for subsequent analyses. Contigs less than 1000 bp in the SAG were removed. We removed host- and human-derived contaminations, which were annotated as eukaryote by BLAST against the NCBI nt database, from the SAG [39]. The quality of the SAGs was checked using CheckM v1.0.13 [40]. For the subsequent analyses, only high-quality SAG data were used (85% > completeness, < 5% contamination).

Phylogenetic and comparative genome analysis

Gene prediction was performed to estimate protein coding genes (Prodigal v. 2.6.3) [41], tRNAs (ARAGORN v. 1.2.38) [42], and rRNAs (Barrnap v. 0.9) [43] by Prokka v.1.14.5 for the SAGs [44]. Predicted proteins were annotated by diamond v0.9.14 [45] against the NCBI nr database [46]. We also ran InterProScan v.5.40-77.0 [47] and eggNOG-mapper v.1.0.3 [48] to annotate domain information and GO terms and KEGG orthology. Enrichment analysis was performed using the clusterProfiler package [49].

The orthogroup was estimated using OrthoFinder v2.3.0 with the default parameters [50]. Single-copy marker genes were extracted based on the OrthoFinder results. Each gene was aligned by MAFFT v7.407 and trimmed by trimAl v1.4 [27, 51]. A phylogenetic tree was constructed using the single-copy marker genes by IQ-TREE v1.6.7 [52] with ultrafast bootstrap [53] and ModelFinder [54]. Average Nucleotide Identity (ANI) and Average Amino Acid Identity (AAI) were calculated by the all-against-all ANI/AAI matrix calculator [55]. Synteny plots were created by gggenomes [56].

Eukaryotic-like protein analysis

All the Pfam IDs with a Z-score of 10,000 in EffectiveELD version 5.2 were extracted from the InterProScan results [57]. For extracting eukaryotic-like genes based on sequence homology, genes were selected by taxonomic annotation using eggNOG-mapper v.1.0.3 [48]. In addition, we performed a taxonomy annotation of eukaryotic-like genes using MMseqs2 with “--tax-lineage 1 --vote-mode 0” options, using the eukaryotic-only NCBI nr database [58]. We also checked that these genes are not annotated with other bacteria to confirm that Endozoicomonas acquired the genes from eukaryotes during host adaptation using the NCBI nr database without the Endozoicomonas genus. Major hits of annotation result were segmented annelid worms Capitella teleta. These genes were considered to be derived from contamination and removed from the results.

Results

Diversity analysis of the Endozoicomonas genus in Acropora tenuis

A. tenuis were collected at four locations (Ishikawabaru, Sesoko Minami, Yomitan, and Onna) in Okinawa in May and August 2016 (Fig. 1A). These four locations are placed within a 15-km radius. Through 16S rRNA gene sequencing, 1418 ASVs were detected as coral-associated bacteria, and 209 ASVs were classified as Endozoicomonas spp. In Ishikawabaru, Yomitan, and Onna, Endozoicomonas spp. was dominant, representing 70.7–99.7% of the detected bacteria (Fig. 1B). However, in Sesoko Minami, the proportion of Endozoicomonas spp. was only 1.8–50.3%, and order Rickettsiales was dominant instead (Supplementary Figure 2A), including Rickettsiales bacteria strain SESOKO1, which was previously reported by our group [59]. The results of principal coordinate analysis (PCoA) based on the Bray-Curtis distance of Endozoicomonas composition showed spatial specificity for A. tenuis (Supplementary Figure 2B and C).

Phylogenetic analysis of the 16S rRNA gene revealed that the Endozoicomonas species can be divided into three major clades: Clade A1 (Red) and Clade A2 (Green) constituted the Onna and Yomitan samples, whereas Clade B (Blue) comprised the Ishikawabaru and Sesoko Minami samples (Fig. 1C, phylogenetic tree and heatmap). Absolute abundance of 16S rRNA gene was calculated by the 16S rRNA gene copy number measured by ddPCR and the surface area value obtained by CT scan. The average (±sd) number of Endozoicomonas was 37,952 ± 75,836 copies/mm2 (Ishikawabaru), 2306.8 ± 4241.6 copies/mm2 (Sesoko Minami), 30,180 ± 30,165 copies/mm2 (Onna), and 37,287 ± 29,083 copies/mm2 (Yomitan), respectively (Fig. 1C, bar plot). These results suggest that the cell numbers of Endozoicomonas spp. vary drastically from site to site. Especially at Sesoko Minami, where the copy number of 16S rRNA gene was exceedingly low, the absolute number of Endozoicomonas was significantly lower than that at other sites. In August 2016, large-scale coral bleaching was observed around Okinawa Prefecture; however, the bleaching level of coral colonies at Sesoko Minami was lower than that at Ishikawabaru [60]. In particular, the single coral colony, in which Endozoicomonas Clade A1 was dominant, did not show any bleaching in August 2016 (Fig. 1C; star). In our previous study, the density per surface area of zooxanthellae in corals was greatly reduced during coral bleaching in August 2016 [60]. However, in the same coral, the absolute number of Endozoicomonas did not decrease significantly from May to August in 2016 (Fig. 1C).

Single-cell genome sequencing of Endozoicomonas bacteria

To selectively obtain genomic information of target bacteria from crushed coral tissue suspensions, we have developed a method for targeted single-cell genomics using droplet microfluidics. Bacterial suspensions prepared from coral tissues were encapsulated into agarose gel droplets and processed to gel-bead-based WGA to obtain SAGs [24] (Fig. 2A). Microscopic observation revealed that 24.4–26.3% of the gel beads exhibited DAPI fluorescence, which corresponded to the Poisson distribution value. Subsequently, PCR-based detection of target genes with FAM probes was conducted, which enabled the identification of gel beads containing target SAGs (Fig. 2B). After PCR amplification, there were droplets displaying both DAPI and FAM fluorescence (Fig. 2C). The percentage of DAPI fluorescence-positive droplets is 0.8–1.2%, which represents the rate of droplets encapsulating the target bacterium.

Fig. 2
figure 2

Targeted single-cell genomics with droplet microfluidics. A Microbial cells were encapsulated into 40-μm microfluidic droplets at the single-cell level (0.3 cell/droplet) with ultra-low melting temperature agarose solutions. Collected droplets were solidified, treated with lysis solution, and subjected to whole genome amplification (WGA) using a previously described single-cell amplified genome (SAG)-gel platform. After washing, whole genome amplified gel beads were resuspended in the PCR mixture. B Gel beads suspended in PCR mixture were reintroduced into microfluidic devices, and water in oil (W/O) droplets encapsulating gel beads were generated. After thermal cycling, droplets with target sequence exhibited fluorescence derived from Taqman probes. The fluorescence-positive droplets were isolated into multi-well plates and then subjected to second-round WGA and library preparation. C Microscopic images of W/O droplets after PCR. Droplets with genome amplification by WGA exhibit blue fluorescence derived from DAPI. Droplets with target gene amplification by PCR exhibit green fluorescence derived from FAM. The droplet indicated by the arrowhead contains genome amplicons of the target bacterium

Full size image

The results of 16S rRNA gene sequencing followed by second-round WGA showed that 64% (28 out of 44) of the samples corresponded to the Endozoicomonas sequence. Among the remaining 16 samples, 13 exhibited insufficient DNA yield (< 2 ng/μL) after the second-round WGA. The remaining three were confirmed to be derived from host mitochondria. When droplets that were positive only for DAPI were isolated, all the sequence corresponded to the sequence of host mitochondria (n = 30), suggesting that our targeted single-cell genomics has a low false-negative rate. In order to obtain 28 target bacterial SAGs by random sampling, 2333–3500 droplets need to be isolated according to the calculation, which means that our method showed 53.0–79.5 times higher screening efficiency than the random sampling.

As a result, seven draft SAGs ranging from 2.1 to 6.4 Mbp were acquired from four sampling sites, three of which were classified as high-quality and three of which were classified as medium-quality according to the criteria by the Genomic Standards Consortium (GSC) (Table 1).

Table 1 Overview of Endozoicomonas
Full size table

Comparative genome analysis of Endozoicomonas genus

We conducted a comparative genome analysis using our five draft SAGs, which include four from this study and one from our previous report [60], collected from A. tenuis, and 10 public Endozoicomonas isolates collected from various kinds of marine organisms. Phylogenetic trees constructed with single-copy marker genes showed that Endozoicomonas spp. represented single clades that corresponded to symbiotic host organisms. However, Endozoicomonas spp. collected from Acropora coral (Endozoicomonas sp. ISHI1, SESOKO1, ONNA1, ONNA2, and YOMI1) were divided into two clades (clade A1 and B) (Fig. 3A), which was also confirmed by the ANI and average AAI (Supplementary Figure 3).

Fig. 3
figure 3

Comparative genome analysis of Endozoicomonas. A Plot of lineage-specific genes based on ortholog analysis. Bar plots indicate the number of orthogroup types. Clade A1 and clade B correspond to the clades in Fig. 1C based on the similarity of the 16S rRNA genes. The color of the circles and species names indicate the type of host organism. The colors of the bar plot indicate the genes used in the KEGG enrichment analysis shown in B (red: clade A1 specific) and C (blue: clade B specific). B Enrichment analysis results of Endozoicomonas clade A1-specific genes for KEGG pathways. The top 10 pathways are shown. The size of the circle indicates the number of genes. The thickness of the lines between pathways indicates the number of shared genes. C Enrichment analysis results of Endozoicomonas Clade B-specific genes for KEGG pathways

Full size image

A similar trend was observed in the ortholog analysis. Out of a total of 9414 orthogroups, the number of orthogroups commonly detected in all Endozoicomonas spp. was 1148. We confirmed that the bacterial secretion systems are evolutionarily conserved in all Endozoicomonas symbiotic relationships with marine organisms (Supplementary Table 3), which are important for Endozoicomonas to secrete proteins related to eukaryotic pathways in the host. The prediction of effector proteins showed that 12.5–20.4% (T3SS) and 2.97–4.22% (T4SS) of proteins were released by bacterial secretion systems (Supplementary Table 4). On the other hand, the remaining 8266 orthogroups and 8061 unassigned genes were clade- or strain-specific, suggesting that Endozoicomonas spp. possess a diverse orthogroup profile. It appears that the profile of the orthogroups depends on the host of each Endozoicomonas sp. However, clade A1 and clade B, which were symbiotic with the same host (A. tenuis), shared no unique orthogroups.

Enrichment analysis of accessory genes to the KEGG pathway revealed the functional characteristics and differences of each clade (clade A1 and clade B). Of the top ten enriched pathways (nine pathway in Clade B) of each clade, eight pathways (except bacterial chemotaxis and flagellar assembly) in clade A1 as well as seven pathways (except sphingolipid metabolism and terpenoid backbone biosynthesis) in clade B were assigned as eukaryotic pathways (Fig. 3B, C). Primarily, in Endozoicomonas clade A1, immune-related pathways were enriched, such as the TRAF signaling pathways and NF-κB pathways (Fig. 3B), which are generally not conserved in bacteria. Conversely, in clade B, these immune-related pathways were not enriched, but proteasomes and proteolysis pathways including eukaryotic-like E3 ubiquitin ligase gene were enriched (Fig. 3C).

Different symbiotic strategies utilizing eukaryotic-like genes

We further conducted a comparative genome analysis focusing on eukaryotic-like domains and found that all Endozoicomonas spp. used in this study had at least one eukaryotic-like domain (Fig. 4A). The types of eukaryotic-like domains were different between clade A1 (YOMI1, ONNA1, and ONNA2) and clade B (ISHI1 and SESEKO1) Endozoicomonas; namely, the eukaryotic-like domain found in common between clade A1 and clade B was only RING-type zinc finger (RINC-type zinc-finger, Zinc-finger double domain, Zinc finger, and C3HC4 type were grouped), and the remaining six domains (SOCS box, Inhibitor of Apoptosis domain, ephrin ligand, Ets domain, MATH domain, and FANCLE C-terminal domain) were found only in either species (Fig. 4A, Supplementary Table 5). We have previously reported that Endozoicomonas sp. ISHI1 has coral-like ephrin ligand genes [60]. In addition, there were other eukaryotic-like genes (such as the SNARE domain, Vps4 C terminal oligomerization domain-containing genes). In contrast, Endozoicomonas sp. YOMI1, ONNA1, and ONNA2, which belong to clade A1, did not possess ephrin ligand, but had eukaryotic-like genes, including MATH/TRAF, SOCS box, an Ets-domain, and an inhibitor of apoptosis domain (Fig. 4A).

Fig. 4
figure 4

Eukaryotic-like genes distribution and horizontal gene transfer in Endozoicomonas. A Properties of eukaryotic-like domains predicted using EffectiveDB in Endozoicomonas (Z-score = 10,000). The phylogenetic tree based on the single-copy marker gene is the same as that in Fig. 3A. The hierarchical clustering at the top is based on the profiles of eukaryotic-like genes. B Insertion of the inhibitor of apoptosis domain-containing genes of Endozoicomonas sp. YOMI1 and Endozoicomonas acroporae Acr-14 into the genome. The inserted eukaryotic-like gene is marked by a yellow square. The region where the GC content of Endozoicomonas ascidiella DSM22380 is flat is the gap region (N)

Full size image

Furthermore, we found 532 genes that did not contain eukaryotic-like domains but were highly homologous against eukaryote in A. tenuis-symbiotic Endozoicomonas (YOMI1: 125 genes, ONNA1: 134 genes, ONNA2: 47 genes, SESOKO1: 90 genes, ISHI1: 136 genes) (Supplementary Table 6). As genes closely related to host A. tenuis coral genes, TRAF/MATH domain-containing genes (match: 26.0–62.2%), low-density lipoprotein receptor gene (match: 30.9–67.6%), and transcription factor 15-like gene (match: 27.5–30.4%) were detected in Endozoicomonas clade A1. In contrast, the genes for ephrin ligand (match: 31.3–46.0%) and RING finger and CHY zinc finger domain-containing gene (match: 46.1–64.3%) were detected in clade B. Eukaryotic-like E3 ubiquitin-protein ligase was detected as a common gene in both clades (Supplementary Table 5). Some of these eukaryotic-like genes were inserted into the genome (e.g., survivin gene, ephrin gene) (Fig. 4B). Furthermore, these genes did not contain introns, despite their similarity to the intron-containing host genes (Supplementary Figure 4).

Discussion

The absolute abundance of Endozoicomonas varies greatly among coral habitats

Previous studies have focused solely on bacterial composition and failed to establish a relationship between the amount of coral symbiotic bacteria and their host conditions. To overcome this limitation, this study quantified the absolute 16S rRNA gene copy number of bacteria and revealed the abundance of Endozoicomonas spp., which varied by more than 72-fold among sampling sites. Several researchers have considered that the abundance of Endozoicomonas is related to health and robustness of corals [65, 66]. However, in August 2016, corals from Ishikawabaru, which had a higher abundance of Endozoicomonas than Sesoko Minami, showed far extensive bleaching, which was not consistent with their findings [67]. Our previous studies also showed that inflammation-related genes are upregulated in A. tenuis, which is occupied by clade B Endozoicomonas spp. In this study, we measured quantitative Endozoicomonas abundance, which may be involved in the host side response at the gene expression level. The individual adaptive capacity of Endozoicomonas may also have an effect on the host [68]. For a deeper understanding of Endozoicomonas function in Acropora coral, the difference of the clade or absolute amount of Endozoicomonas should be considered.

Targeted single-cell genomics enables selective accumulation of low-abundance bacterial genomes

In addition to the 16S rRNA gene-based phylogenetic analysis, we conducted single-cell whole genome sequencing of target bacteria and collected seven SAGs belonging to clade A1 (Endozoicomonas sp. YOMI1, ONNA1, and ONNA2) and clade B (Endozoicomonas sp. SESOKO1 and ISHI1). Although SAGs consisted of fragmented contigs, and the total length of some SAGs were shorter compared to the genome size of reference Endozoicomonas, the completeness was comparable to cultured strains reported previously [15], providing evidence that this method can expand genomic information of uncultured bacteria. It is noteworthy that our method enabled selective acquisition of draft SAGs even when the abundance ratio of target bacteria was approximately 1%, thereby overcoming the disadvantage of metagenomics. Moreover, even in the conventional single-cell genomics, the composition of sequenced bacterial taxa is highly dependent on the proportion of bacteria in the target environment because single cells are randomly sorted [69]. In contrast, in our proposed method, WGA was followed by detection of target gene sequence (Endozoicomonas 16S rRNA gene) in gel beads to achieve selective isolation and subsequent genome analysis. To the best of our knowledge, this is the first report of targeted single-cell genomics on environmental microbes. Our method showed > 50 times higher screening efficiency of target bacteria than the random sampling, thus drastically reducing the amount of reagents and labor wasted on non-target samples. Our targeted single-cell genomics can also be used to detect specific gene sequences obtained from metagenomic shotgun analysis or 16S rRNA gene sequencing by designing specific primers and probes. As a further application, specific single-cell genomics with unique characteristics will be possible by targeting the sequences of biosynthetic gene clusters and drug-resistant genes. As the potential of rare bacteria to exert important functions in microbial populations is gaining increasing attention [70], our method can be a technological innovation and breakthrough for the analysis of single-cell genomics. In the single-cell genomics of environmental bacteria, it is often important to process the samples immediately after collection to obtain high-quality draft SAGs. Because the microfluidic system used in this study is portable, on-site single-cell isolation and genome amplification can be performed with standard biological laboratory equipment, which is critical for obtaining high-quality SAGs.

Diversity of host adaptation strategies in Endozoicomonas

Several 16S rRNA gene amplicon and shotgun metagenomic analyses have suggested that coral-associated Endozoicomonas spp. may improve host coral health [65, 66]. However, owing to the lack of reference genomes, detailed analysis of clade-specific functions in Endozoicomonas in the same coral species has not been performed. Therefore, we used our SAG data and compared the gene functions detected in each clade of Endozoicomonas. Comparative genome analysis revealed that the orthogroups of Endozoicomonas were conserved in each clade. Notably, even when the different clades of Endozoicomonas coexisted with A. tenuis, the shared orthogroups were only the core gene between clade A1 (Endozoicomonas sp. YOMI1, ONNA1, and ONNA2) and clade B (Endozoicomonas sp. ISHI1 and SESOKO1). These results suggest that Endozoicomonas independently adapts to its host during the process of evolution.

A comprehensive search for eukaryotic-like genes in Endozoicomonas revealed that various eukaryotic-like genes were enriched in a lineage-specific manner. Symbiotic bacteria are known to harbor eukaryotic-like domains, such as ankyrin repeat and WD40 repeat, which are involved in protein-protein interactions and are thought to be necessary for interaction with the host [71, 72]. For example, in plants and sponges, many symbiotic bacteria and pathogens intervene in host signaling pathways via eukaryotic-like genes that mimic host genes, thereby influencing host conditions [73, 74]. There have been few reports on how symbiotic bacteria acquire these eukaryotic-like genes, but in our study, some of these genes (e.g., ephrin ligand and survivin) were found to be inserted with spliced host coral genes (Supplementary Figure 4). This result suggests the bacterial-eukaryotic horizontal gene transfer is observed between host coral and coral-associated Endozoicomonas spp. One of many possible hypotheses is that the insertion of these genes may be mediated by host mRNA. Horizontal gene transfer between eukaryotes and prokaryotes has been studied in recent years, and more detailed studies are needed in the future.

A recent study reported that the copy number of ankyrin repeats is higher in E. acroporae genomes than in other marine organism-associated Endozoicomonas genomes [10], implying that some eukaryotic-like genes play essential roles in the symbiosis of Endozoicomonas with Acropora coral. On the basis of this idea, we hypothesized the functions and host adaptation strategies of each Endozoicomonas clade based on the functions of enriched eukaryotic-like domains. Recently, it was hypothesized that ephrin ligands are required for E. montiporae to transition the symbiotic mechanisms [14]. In our previous report, we have also confirmed a gene expansion of ephrin ligand in Endozoicomonas Clade B (Endozoicomonas sp. ISHI1 and SESOKO1) with a higher degree than that in E. montiporae [60]. Interestingly, this is consistent with the fact that the number of ephrin-like receptor genes is expanded in host Acropora coral compared with that in Montipora coral [75]. These results reinforce the hypothesis that Endozoicomonas uses ephrin ligand gene as a part of its adaptation strategy to its host coral. The consistency of ephrin ligand and host ephrin receptor gene expansion suggests the co-evolution between Endozoicomonas and host coral. Furthermore, Endozoicomonas clade B encodes genes associated with intracellular invasion, such as SNARE domains, which contribute to the inhibition of degradation by phagosomes in Chlamydia and Legionella [76]. However, Endozoicomonas clade B does not encode immune-related genes. This result is consistent with our recent finding that Endozoicomonas sp. ISHI1 stimulates immune-related pathways in corals through infection (Fig. 5A) [60]. In addition, genes containing the MATH/TRAF, SOCS box, and inhibitor of apoptosis domain were present in clade A1, but not in clade B. These genes suppress apoptosis induced by the innate immune pathway, such as the NF-κB and JAK-STAT pathways [77]. Recently, several studies reported that the innate immune pathways are activated in host corals as stress responses to bleaching [78]. In particular, the NF-κB pathway is conserved in cnidarians and has attracted much attention in the analysis of bleaching [79,80,81]. Therefore, we speculate that a group of eukaryotic-like genes that are significantly enriched in Endozoicomonas clade A1 may affect host innate immune pathways and increase tolerance to environmental stress.

Fig. 5
figure 5

Proposed symbiotic strategy of Endozoicomonas acroporae in A. tenuis, as estimated by comparative genome analysis. A Relationship between clade B Endozoicomonas and A. tenuis corals. Endozoicomonas possesses many genes related to intracellular invasion, such as SNARE domain-containing genes, in addition to the eukaryotic-like gene ephrin ligand, but it does not possess any genes related to apoptosis. B Relationships between Endozoicomonas clade A1 and host A. tenuis corals. Endozoicomonas clade A1 possesses several eukaryotic-like genes involved in the suppression of apoptosis but does not have the coral-like ephrin ligand

Full size image

Homology-based analysis also revealed several eukaryotic-like genes with unknown functions in Endozoicomonas clade A1 and B. These genes are phylogenetically conserved and may have some function in host adaptation and symbiosis. One example is the low-density lipoprotein receptor (LDLR) gene conserved in clade A1. Recently, the function of LDLR during oogenesis in A. tenuis has been reported; it has been shown that the transcript of AtLDLR is present in the mesentery, mesenteric membrane, and mesenteric filaments surrounding an oocyte, but not in the membrane of a developing oocyte. AtLDLR can bind to major yolk lipoproteins and may be transported into the cytoplasm of an oocyte as it develops [82]. This indicates that the LDLR-like gene of Endozoicomonas clade A1 has the binding function to host egg lipoprotein. Therefore, we speculate that the LDLR-like gene, which is involved in egg development, may contribute to vertical transfer of Endozoicomonas, as other infection-related genes such as ephrin ligand gene and SNARE gene have not been identified in clade A1. The speculation can be verified by future experiments.

Conclusions

We succeeded in the targeted single-cell genomics of Endozoicomonas from A. tenuis corals with droplet microfluidics. Thus, we accomplished obtaining high-quality Endozoicomonas SAGs (clade A1 and B) from the same coral species, which were proved to be phylogenetically and functionally diverse. Comparative genome analysis revealed the difference in the eukaryotic-like gene profiles, suggesting that the symbiotic strategy of Endozoicomonas in A. tenuis is divided into two types. Our findings suggest that even in the same bacterial genus, each bacterial lineage has different gene profiles and functions. Therefore, it is necessary to conduct strain-level genome sequencing to analyze the functions and strategies for host adaptation. Although in vitro studies in the field environment are required to prove this hypothesis, our targeted single-cell genomics could be a powerful tool for acquiring SAGs of uncultured bacterium to infer its function. Analysis of various coral-associated Endozoicomonas SAGs will clarify their host adaptation strategies in more detail and shed light on some aspects of the symbiotic mechanisms.

Availability of data and materials

The assembled genomes of seven Endozoicomonas were deposited in DDBJ/ENA/GenBank under BioProject number PRJNA721663. The raw read data of 16S rRNA gene are available under BioProject number PRJNA722798.

References

  1. Muscatine L, Porter JW. Reef corals: mutualistic symbioses adapted to nutrient-poor environments. Bioscience. 1977;27:454–60.

    Article  Google Scholar 

  2. Liu H, Stephens TG, González-Pech RA, Beltran VH, Lapeyre B, Bongaerts P, et al. Symbiodinium genomes reveal adaptive evolution of functions related to coral-dinoflagellate symbiosis. Commun Biol. 2018;1:95.

    Article  Google Scholar 

  3. Neave MJ, Rachmawati R, Xun L, Michell CT, Bourne DG, Apprill A, et al. Differential specificity between closely related corals and abundant Endozoicomonas endosymbionts across global scales. Isme J. 2016;11:186–200.

    Article  Google Scholar 

  4. Lema KA, Willis BL, Bourne DG. Corals form characteristic associations with symbiotic nitrogen-fixing bacteria. Appl Environ Microb. 2012;78:3136–44.

    Article  CAS  Google Scholar 

  5. Raina J-B, Tapiolas D, Motti CA, Foret S, Seemann T, Tebben J, et al. Isolation of an antimicrobial compound produced by bacteria associated with reef-building corals. Peerj. 2016;4:e2275.

    Article  Google Scholar 

  6. Cano I, van Aerle R, Ross S, Verner-Jeffreys DW, Paley RK, Rimmer GSE, et al. Molecular characterization of an Endozoicomonas-like organism causing infection in the king scallop (Pecten maximus L.). Appl Environ Microb. 2018;84:e00952–17.

    Article  Google Scholar 

  7. Hyun DW, Shin NR, Kim MS, Oh SJ, Kim PS, Whon TW, et al. Endozoicomonas atrinae sp. nov., isolated from the intestine of a comb pen shell Atrina pectinata. Int J Syst Evol Micr. 2014;64:2312–8.

    Article  CAS  Google Scholar 

  8. Kurahashi M, Yokota A. Endozoicomonas elysicola gen. nov., sp. nov., a γ-proteobacterium isolated from the sea slug Elysia ornata. Syst Appl Microbiol. 2007;30:202–6.

    Article  CAS  Google Scholar 

  9. Nishijima M, Adachi K, Katsuta A, Shizuri Y, Yamasato K. Endozoicomonas numazuensis sp. nov., a gammaproteobacterium isolated from marine sponges, and emended description of the genus Endozoicomonas Kurahashi and Yokota 2007. Int J Syst Evol Micr. 2013;63:709–14.

    Article  CAS  Google Scholar 

  10. Alex A, Antunes A. Comparative genomics reveals metabolic specificity of Endozoicomonas isolated from a marine sponge and the genomic repertoire for host-bacteria symbioses. Microorg. 2019;7:635.

    Article  CAS  Google Scholar 

  11. Du Z, Zhang W, Xia H, Lü G, Chen G. Isolation and diversity analysis of heterotrophic bacteria associated with sea anemones. Acta Oceanol Sin. 2010;29:62–9.

    Article  CAS  Google Scholar 

  12. Schuett C, Doepke H, Grathoff A, Gedde M. Bacterial aggregates in the tentacles of the sea anemone Metridium senile. Helgoland Mar Res. 2007;61:211–6.

    Article  Google Scholar 

  13. Glasl B, Smith CE, Bourne DG, Webster NS. Disentangling the effect of host-genotype and environment on the microbiome of the coral Acropora tenuis. Peerj. 2019;7:e6377.

    Article  Google Scholar 

  14. Ding J-Y, Shiu J-H, Chen W-M, Chiang Y-R, Tang S-L. Genomic insight into the host–endosymbiont relationship of Endozoicomonas montiporae CL-33T with its coral host. Front Microbiol. 2016;7(e1000376):e1000376-15.

    Google Scholar 

  15. Tandon K, Lu C-Y, Chiang P-W, Wada N, Yang S-H, Chan Y-F, et al. Comparative genomics: Dominant coral-bacterium Endozoicomonas acroporae metabolizes dimethylsulfoniopropionate (DMSP). Isme J. 2020;14:1290–303.

    Article  CAS  Google Scholar 

  16. Neave MJ, Michell CT, Apprill A, Voolstra CR. Endozoicomonas genomes reveal functional adaptation and plasticity in bacterial strains symbiotically associated with diverse marine hosts. Sci Rep-uk. 2017;7(1):1–12.

    Google Scholar 

  17. Yue Y, Huang H, Qi Z, Dou H-M, Liu X-Y, Han T-F, et al. Evaluating metagenomics tools for genome binning with real metagenomic datasets and CAMI datasets. BMC Bioinformatics. 2020;21:334.

    Article  CAS  Google Scholar 

  18. Quince C, Delmont TO, Raguideau S, Alneberg J, Darling AE, Collins G, et al. DESMAN: a new tool for de novo extraction of strains from metagenomes. Genome Biol. 2017;18:181.

    Article  Google Scholar 

  19. van de Water JAJM, Allemand D, Ferrier-Pagès C. Host-microbe interactions in octocoral holobionts - recent advances and perspectives. Microbiome. 2018;6:64.

    Article  Google Scholar 

  20. Feehery GR, Yigit E, Oyola SO, Langhorst BW, Schmidt VT, Stewart FJ, et al. A method for selectively enriching microbial DNA from contaminating vertebrate host DNA. Plos One. 2013;8:e76096.

    Article  CAS  Google Scholar 

  21. Schmieder R, Edwards R. Fast identification and removal of sequence contamination from genomic and metagenomic datasets. Plos One. 2011;6:e17288.

    Article  CAS  Google Scholar 

  22. Nishikawa Y, Hosokawa M, Maruyama T, Yamagishi K, Mori T, Takeyama H. Monodisperse picoliter droplets for low-bias and contamination-free reactions in single-cell whole genome amplification. Plos One. 2015;10:e0138733.

    Article  Google Scholar 

  23. Hosokawa M, Nishikawa Y, Kogawa M, Takeyama H. Massively parallel whole genome amplification for single-cell sequencing using droplet microfluidics. Sci Rep-uk. 2017;7:5199.

    Article  Google Scholar 

  24. Chijiiwa R, Hosokawa M, Kogawa M, Nishikawa Y, Ide K, Sakanashi C, et al. Single-cell genomics of uncultured bacteria reveals dietary fiber responders in the mouse gut microbiota. Microbiome. 2020;8:5.

    Article  CAS  Google Scholar 

  25. Bolyen E, Rideout JR, Dillon MR, Bokulich NA, Abnet CC, Al-Ghalith GA, et al. Reproducible, interactive, scalable and extensible microbiome data science using QIIME 2. Nat Biotechnol. 2019;37:852–7.

    Article  CAS  Google Scholar 

  26. Callahan BJ, McMurdie PJ, Rosen MJ, Han AW, Johnson AJA, Holmes SP. DADA2: High-resolution sample inference from Illumina amplicon data. Nat Methods. 2016;13:581–3.

    Article  CAS  Google Scholar 

  27. Katoh K, Misawa K, Kuma K, Miyata T. MAFFT: a novel method for rapid multiple sequence alignment based on fast Fourier transform. Nucleic Acids Res. 1911;30:3059–66.

    Article  Google Scholar 

  28. Price MN, Dehal PS, Arkin AP. FastTree 2 – approximately maximum-likelihood trees for large alignments. Plos One. 2010;5:e9490.

    Article  Google Scholar 

  29. Bokulich NA, Kaehler BD, Rideout JR, Dillon M, Bolyen E, Knight R, et al. Optimizing taxonomic classification of marker-gene amplicon sequences with QIIME 2’s q2-feature-classifier plugin. Microbiome. 2018;6:90.

    Article  Google Scholar 

  30. Quast C, Pruesse E, Yilmaz P, Gerken J, Schweer T, Yarza P, et al. The SILVA ribosomal RNA gene database project: improved data processing and web-based tools. Nucleic Acids Res. 2013;41:D590–6.

    Article  CAS  Google Scholar 

  31. McMurdie PJ, Holmes S. phyloseq: an R package for reproducible interactive analysis and graphics of microbiome census data. Plos One. 2013;8:e61217.

    Article  CAS  Google Scholar 

  32. Yu G, Smith DK, Zhu H, Guan Y, Lam TT. ggtree: an r package for visualization and annotation of phylogenetic trees with their covariates and other associated data. Methods Ecol Evol. 2017;8:28–36.

    Article  Google Scholar 

  33. Schneider CA, Rasband WS, Eliceiri KW. NIH Image to ImageJ: 25 years of image analysis. Nat Methods. 2012;9:671–5.

    Article  CAS  Google Scholar 

  34. Nishikawa Y, Kogawa M, Hosokawa M, Wagatsuma R, Mineta K, Takahashi K, Ide K, Yura K, Behzad H, Gojobori T, Takeyama H. Validation of the application of gel beads-based single-cell genome sequencing platform to soil and seawater. ISME commun. 2022;2:92.

  35. Ye J, Coulouris G, Zaretskaya I, Cutcutache I, Rozen S, Madden TL. Primer-BLAST: a tool to design target-specific primers for polymerase chain reaction. BMC Bioinformatics. 2012;13:134.

    Article  CAS  Google Scholar 

  36. Chen S, Zhou Y, Chen Y, Gu J. fastp: an ultra-fast all-in-one FASTQ preprocessor. Bioinformatics. 2018;34:i884–90.

    Article  Google Scholar 

  37. Bankevich A, Nurk S, Antipov D, Gurevich AA, Dvorkin M, Kulikov AS, et al. SPAdes: a new genome assembly algorithm and its applications to single-cell sequencing. J Comput Biology J Comput Mol Cell Biol. 2012;19:455–77.

    Article  CAS  Google Scholar 

  38. Kogawa M, Hosokawa M, Nishikawa Y, Mori K, Takeyama H. Obtaining high-quality draft genomes from uncultured microbes by cleaning and co-assembly of single-cell amplified genomes. Sci Rep-uk. 2018;8:2059.

    Article  Google Scholar 

  39. McGinnis S, Madden TL. BLAST: at the core of a powerful and diverse set of sequence analysis tools. Nucleic Acids Res. 2004;32:W20–5.

    Article  CAS  Google Scholar 

  40. Parks DH, Imelfort M, Skennerton CT, Hugenholtz P, Tyson GW. CheckM: assessing the quality of microbial genomes recovered from isolates, single cells, and metagenomes. Genome Res. 2015;25:1043–55.

    Article  CAS  Google Scholar 

  41. Hyatt D, Chen G-L, LoCascio PF, Land ML, Larimer FW, Hauser LJ. Prodigal: prokaryotic gene recognition and translation initiation site identification. BMC Bioinformatics. 2010;11:119–11.

    Article  Google Scholar 

  42. Laslett D. ARAGORN, a program to detect tRNA genes and tmRNA genes in nucleotide sequences. Nucleic Acids Res. 2004;32:11–6.

    Article  CAS  Google Scholar 

  43. Seemann T. barrnap 0.9 : rapid ribosomal RNA prediction. https://github.com/tseemann/barrnap.

  44. Seemann T. Prokka: rapid prokaryotic genome annotation. Bioinformatics. 2014;30:2068–9.

    Article  CAS  Google Scholar 

  45. Buchfink B, Xie C, Huson DH. Fast and sensitive protein alignment using DIAMOND. Nat Methods. 2015;12:59–60.

    Article  CAS  Google Scholar 

  46. Coordinators NR, Agarwala R, Barrett T, Beck J, Benson DA, Bollin C, et al. Database resources of the National Center for Biotechnology Information. Nucleic Acids Res. 2017;46:D8–D13.

    Article  Google Scholar 

  47. Jones P, Binns D, Chang HY, Fraser M, Li W, McAnulla C, et al. InterProScan 5: genome-scale protein function classification. Bioinformatics. 2014;30:1236–40.

    Article  CAS  Google Scholar 

  48. Huerta-Cepas J, Forslund K, Coelho LP, Szklarczyk D, Jensen LJ, von Mering C, et al. Fast genome-wide functional annotation through orthology assignment by eggNOG-Mapper. Mol Biol Evol. 2017;34:2115–22.

    Article  CAS  Google Scholar 

  49. Yu G, Wang L-G, Han Y, He Q-Y. clusterProfiler: an R package for comparing biological themes among gene clusters. Omics J Integr Biology. 2012;16:284–7.

    Article  CAS  Google Scholar 

  50. Emms DM, Kelly S. OrthoFinder: phylogenetic orthology inference for comparative genomics. Genome Biol. 2019;20:238.

    Article  Google Scholar 

  51. Capella-Gutierrez S, Silla-Martinez JM, Gabaldon T. trimAl: a tool for automated alignment trimming in large-scale phylogenetic analyses. Bioinformatics. 2009;25:1972–3.

    Article  CAS  Google Scholar 

  52. Minh BQ, Schmidt HA, Chernomor O, Schrempf D, Woodhams MD, von Haeseler A, et al. IQ-TREE 2: New models and efficient methods for phylogenetic inference in the genomic era. Mol Biol Evol. 2020;37:1530–4.

    Article  CAS  Google Scholar 

  53. Hoang DT, Chernomor O, von Haeseler A, Minh BQ, Vinh LS. UFBoot2: improving the ultrafast bootstrap approximation. Mol Biol Evol. 2017;35:518–22.

    Article  Google Scholar 

  54. Kalyaanamoorthy S, Minh BQ, Wong TKF, von Haeseler A, Jermiin LS. ModelFinder: fast model selection for accurate phylogenetic estimates. Nat Methods. 2017;14:587–9.

    Article  CAS  Google Scholar 

  55. Rodriguez-R LM, Konstantinidis KT. The enveomics collection: a toolbox for specialized analyses of microbial genomes and metagenomes. PeerJ Preprints. 2016:4:e1900v1.

  56. Thomas Hackl and Markus J. Ankenbrand (2020). gggenomes: a grammar of graphics for comparative genomics. R package version 0.0.0.9000. https://github.com/thackl/gggenomes

    Google Scholar 

  57. Eichinger V, Nussbaumer T, Platzer A, Jehl M-A, Arnold R, Rattei T. EffectiveDB—updates and novel features for a better annotation of bacterial secreted proteins and Type III, IV, VI secretion systems. Nucleic Acids Res. 2016;44:D669–74.

    Article  CAS  Google Scholar 

  58. Steinegger M, Söding J. MMseqs2 enables sensitive protein sequence searching for the analysis of massive data sets. Nat Biotechnol. 2017;35:1026–8.

    Article  CAS  Google Scholar 

  59. Ide K, Nishikawa Y, Kogawa M, Iwamoto E, Samuel AZ, Nakano Y, et al. High-quality draft genome sequence of a Rickettsiales bacterium found in Acropora tenuis coral from Okinawa, Japan. Microbiol Resour Announc. 2020;9:e00848–20.

  60. Maruyama T, Ito M, Wakaoji S, Okubo Y, Ide K, Fujimura H, et al. Multi-omics analysis highlights cross-organism interactions in coral holobiont. bioRvix 2021:10.25.465660. https://doi.org/10.1101/2021.10.25.465660.

  61. Tandon K, Chiang P-W, Chen W-M, Cang S-L. Draft Genome Sequence of Endozoicomonas acroporae Strain Acr-14T, Isolated from Acropora Coral. Genome Announc. 2018;6:e01576–17.

  62. Appolinario LR, Tschoeke DA, Rua CPJ, Venas T, Campeão ME, Amaral GRS, et al. Description of Endozoicomonas arenosclerae sp. nov. using a genomic taxonomy approach. Antonie Van Leeuwenhoek. 2016;109:431–8.

  63. Schreiber L, Kjeldsen KU, Obst M, Funch P, Schramm A. Description of Endozoicomonas ascidiicola sp. nov., isolated from Scandinavian ascidians. Syst Appl Microbiol. 2016;39:313–8.

  64. Neave MJ, Michell CT,Apprill A, Voolstra CR. Whole-Genome Sequences of Three Symbiotic Endozoicomonas Strains. Genome Announc. 2014;2:e00802–14.

  65. Bourne D, Iida Y, Uthicke S, Smith-Keune C. Changes in coral-associated microbial communities during a bleaching event. Isme J. 2008;2:350–63.

    Article  CAS  Google Scholar 

  66. Pogoreutz C, Rädecker N, Cárdenas A, et al. Dominance of Endozoicomonas bacteria throughout coral bleaching and mortality suggests structural inflexibility of the Pocillopora verrucosa microbiome. Ecol Evol. 2018;8:2240–52.

    Article  Google Scholar 

  67. Neave MJ, Apprill A, Ferrier-Pagès C, Voolstra CR. Diversity and function of prevalent symbiotic marine bacteria in the genus Endozoicomonas. Appl Microbiol Biot. 2016;100:8315–24.

    Article  CAS  Google Scholar 

  68. Tandon K, Chiou Y-J, Yu S-P, Hsieh HJ, Lu C-Y, Hsu M-T, et al. Microbiome restructuring: dominant coral bacterium Endozoicomonas species respond differentially to environmental changes. Msystems. 2022;7:e00359-22.

    Article  Google Scholar 

  69. Lan F, Demaree B, Ahmed N, Abate AR. Single-cell genome sequencing at ultra-high-throughput with microfluidic droplet barcoding. Nat Biotechnol. 2017;35:640–6.

    Article  CAS  Google Scholar 

  70. Lynch MDJ, Neufeld JD. Ecology and exploration of the rare biosphere. Nat Rev Microbiol. 2015;13:217–29.

    Article  CAS  Google Scholar 

  71. Mosavi LK, Cammett TJ, Desrosiers DC, Peng Z. The ankyrin repeat as molecular architecture for protein recognition. Protein Sci. 2004;13:1435–48.

    Article  CAS  Google Scholar 

  72. Jain BP, Pandey S. WD40 repeat proteins: signalling scaffold with diverse functions. Protein J. 2018;37:391–406.

    Article  CAS  Google Scholar 

  73. Díez-Vives C, Moitinho-Silva L, Nielsen S, Reynolds D, Thomas T. Expression of eukaryotic-like protein in the microbiome of sponges. Mol Ecol. 2017;26:1432–51.

    Article  Google Scholar 

  74. Levy A, Gonzalez IS, Mittelviefhaus M, Clingenpeel S, Paredes SH, Miao J, et al. Genomic features of bacterial adaptation to plants. Nat Genet. 2018;50:138–50.

    Article  CAS  Google Scholar 

  75. Cunning R, Bay RA, Gillette P, Baker AC, Traylor-Knowles N. Comparative analysis of the Pocillopora damicornis genome highlights role of immune system in coral evolution. Sci Rep-uk. 2018;8:16134.

    Article  CAS  Google Scholar 

  76. Paumet F, Wesolowski J, Garcia-Diaz A, Delevoye C, Aulner N, Shuman HA, et al. Intracellular bacteria encode inhibitory SNARE-like proteins. Plos One. 2009;4:e7375–6.

    Article  Google Scholar 

  77. Wong BR, Josien R, Lee SY, Vologodskaia M, Steinman RM, Choi Y. The TRAF family of signal transducers mediates NF-κB activation by the TRANCE receptor*. J Biol Chem. 1998;273:28355–9.

    Article  CAS  Google Scholar 

  78. Pinzón JH, Kamel B, Burge CA, Harvell CD, Medina M, Weil E, et al. Whole transcriptome analysis reveals changes in expression of immune-related genes during and after bleaching in a reef-building coral. Roy Soc Open Sci. 2015;2:140214.

    Article  Google Scholar 

  79. Williams LM, Fuess LE, Brennan JJ, Mansfield KM, Salas-Rodriguez E, Welsh J, et al. A conserved Toll-like receptor-to-NF-κB signaling pathway in the endangered coral Orbicella faveolata. Dev Comp Immunol. 2018;79:128–36.

    Article  CAS  Google Scholar 

  80. Mansfield KM, Carter NM, Nguyen L, Cleves PA, Alshanbayeva A, Williams LM, et al. Transcription factor NF-κB is modulated by symbiotic status in a sea anemone model of cnidarian bleaching. Sci Rep-uk. 2017;7:16025.

    Article  Google Scholar 

  81. Mansfield KM, Cleves PA, Vlack EV, Kriefall NG, Benson BE, Camacho DJ, et al. Varied effects of algal symbionts on transcription factor NF-κB in a sea anemone and a coral: possible roles in symbiosis and thermotolerance. bioRxiv. 2019:640177. https://doi.org/10.1101/640177.

  82. Tan ES, Izumi R, Takeuchi Y, Isomura N, Takemura A. Molecular approaches underlying the oogenic cycle of the scleractinian coral, Acropora tenuis. Sci Rep-uk. 2020;10:9914.

Download references

Acknowledgements

We would like to thank Dr. Mitsuo Umezu for CT scan analysis, Naoko Okada for 16S rRNA gene sequencing, and Dr. Hisashi Anbutsu and Eisuke Iwamoto for useful suggestions. The super-computing resource was provided by the Human Genome Center (University of Tokyo). This study was supported by the Collaborative Research of Tropical Biosphere Research Center, University of the Ryukyus.

Funding

This work was supported by the “Construction of the environmental risk mathematical model by the meta-omics analyses of marine unculturable microbes based on single cell genome information” grant (JPMJCR12A4) from JST-CREST and by JSPS KAKENHI grant number JP17H06158.

Author information

Author notes

  1. Keigo Ide, Yohei Nishikawa and Toru Maruyama contributed equally to this work.

Authors and Affiliations

  1. Department of Life Science and Medical Bioscience, Waseda University, Tokyo, Japan

    Keigo Ide, Toru Maruyama, Yuko Tsukada, Masato Kogawa, Hiroki Takeda, Haruka Ito, Ryota Wagatsuma, Rimi Miyaoka, Masahito Hosokawa, Kei Yura & Haruko Takeyama

  2. Computational Bio Big-Data Open Innovation Laboratory (CBBD-OIL), National Institute of Advanced Industrial Science and Technology, Tokyo, Japan

    Keigo Ide, Yohei Nishikawa, Yuko Tsukada, Masato Kogawa, Ryota Wagatsuma, Masahito Hosokawa, Kei Yura & Haruko Takeyama

  3. Research Organization for Nano and Life Innovation, Waseda University, Tokyo, Japan

    Keigo Ide, Yohei Nishikawa, Masato Kogawa, Masahito Hosokawa, Kei Yura & Haruko Takeyama

  4. Tropical Biosphere Research Center, University of the Ryukyus, Okinawa, Japan

    Yoshikatsu Nakano & Michihiro Ito

  5. Marine Science Section, Research Support Division, Okinawa Institute of Science and Technology Graduate University, Okinawa, Japan

    Yoshikatsu Nakano

  6. Umino-Tane Co., Ltd., Okinawa, Japan

    Koji Kinjo

  7. Institute for Advanced Research of Biosystem Dynamics, Waseda Research Institute for Science and Engineering, Tokyo, Japan

    Masahito Hosokawa & Haruko Takeyama

  8. Graduate School of Humanities and Sciences, Ochanomizu University, Tokyo, Japan

    Kei Yura

  9. Faculty of Science, University of the Ryukyus, Okinawa, Japan

    Shoichiro Suda

Authors
  1. Keigo Ide
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Yohei Nishikawa
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Toru Maruyama
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Yuko Tsukada
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Masato Kogawa
    View author publications

    You can also search for this author in PubMed Google Scholar

  6. Hiroki Takeda
    View author publications

    You can also search for this author in PubMed Google Scholar

  7. Haruka Ito
    View author publications

    You can also search for this author in PubMed Google Scholar

  8. Ryota Wagatsuma
    View author publications

    You can also search for this author in PubMed Google Scholar

  9. Rimi Miyaoka
    View author publications

    You can also search for this author in PubMed Google Scholar

  10. Yoshikatsu Nakano
    View author publications

    You can also search for this author in PubMed Google Scholar

  11. Koji Kinjo
    View author publications

    You can also search for this author in PubMed Google Scholar

  12. Michihiro Ito
    View author publications

    You can also search for this author in PubMed Google Scholar

  13. Masahito Hosokawa
    View author publications

    You can also search for this author in PubMed Google Scholar

  14. Kei Yura
    View author publications

    You can also search for this author in PubMed Google Scholar

  15. Shoichiro Suda
    View author publications

    You can also search for this author in PubMed Google Scholar

  16. Haruko Takeyama
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

M.I., M.H., S.S., and H.T. designed the experiments. K.I., Y.Ni., Y.T., K.Y., and H.T. wrote all the manuscripts. M.I., K.I., Y.Ni., H.I., K.K., and Y.Na. collected the data. Y.Ni. and M.K. conducted the single-cell genomics experiments for obtaining E. sp. YOMI1, SESOKO1-4, and ONNA1-2 SAGs. H.T., R.M., and H.I. conducted the ddPCR analysis. H.I. conducted CT scan analysis. K.I., M.K., R.W., and T.M performed the bioinformatics analysis. All authors have read and approved the final manuscript.

Corresponding author

Correspondence to Haruko Takeyama.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

M.H. and H.T. are shareholders in bitBiome, Inc., which provides single-cell genomics services using the SAG-gel workflow as bit-MAP. K.I., Y.Ni., K.M., M.H., and H.T. are inventors on patent applications covering the technique for single-cell sequencing.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1: Supplementary Figure 1.

Droplet Digital PCR to quantify 16S rRNA gene copy number of bacteria. Supplementary Figure 2. Order level barplot and PCoA plot based on UniFrac distance of Endozoicomonas genus in four sampling points. Supplementary Figure 3. Average Nucleotide Identity (ANI) and Average Amino Acid Identity (AAI) of Endozoicomonas genomes. Supplementary Figure 4. Comparison of gene structures of Coral-like ephrin ligand genes.

Additional file 2: Table S1.

List of Primer sets used to identify Endozoicomonas. Table S2. Raw read statistics of our acquired Endozoicomonas genus genome. Table S3. List of bacterial secretion systems conserved in Endozoicomonas. Table S4. Percentage of proteins released by the bacterial secretion system. Table S5. List of eukaryotic-like domains in Endozoicomonas. Table S6. List of all eukaryotic-like genes in Endozoicomonas.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and Permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Ide, K., Nishikawa, Y., Maruyama, T. et al. Targeted single-cell genomics reveals novel host adaptation strategies of the symbiotic bacteria Endozoicomonas in Acropora tenuis coral. Microbiome 10, 220 (2022). https://doi.org/10.1186/s40168-022-01395-9

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s40168-022-01395-9

Microbiome

ISSN: 2049-2618

Contact us