Genomic Insight into the Host-Endosymbiont Relationship of Endozoicomonas montiporae CL-33(T) with its Coral Host.

The bacterial genus Endozoicomonas was commonly detected in healthy corals in many coral-associated bacteria studies in the past decade. Although, it is likely to be a core member of coral microbiota, little is known about its ecological roles. To decipher potential interactions between bacteria and their coral hosts, we sequenced and investigated the first culturable endozoicomonal bacterium from coral, the E. montiporae CL-33(T). Its genome had potential sign of ongoing genome erosion and gene exchange with its host. Testosterone degradation and type III secretion system are commonly present in Endozoicomonas and may have roles to recognize and deliver effectors to their hosts. Moreover, genes of eukaryotic ephrin ligand B2 are present in its genome; presumably, this bacterium could move into coral cells via endocytosis after binding to coral's Eph receptors. In addition, 7,8-dihydro-8-oxoguanine triphosphatase and isocitrate lyase are possible type III secretion effectors that might help coral to prevent mitochondrial dysfunction and promote gluconeogenesis, especially under stress conditions. Based on all these findings, we inferred that E. montiporae was a facultative endosymbiont that can recognize, translocate, communicate and modulate its coral host.

Introduction

The bacterial genus Endozoicomonas was first proposed by Kurahashi and Yokota (2007). All culturable species of this genus were isolated from marine invertebrates, including sea slug (Kurahashi and Yokota, 2007), encrusting coral (Yang et al., 2010), sponge (Nishijima et al., 2013), octocorals (Pike et al., 2013), and comb pen shell (Hyun et al., 2014). In culture-independent bacterial community studies, Endozoicomonas were frequently detected in various marine invertebrates and have been considered an important affiliate of the host holobiont, due to their positive association with healthy individuals, particularly in corals. For example, Vezzulli et al. (2013) compared bacterial communities on Mediterranean gorgonians (Paramuricea clavata) and reported that Endozoicomonas were a predominant group on healthy individuals, but declined greatly when their host was diseased. Similarly, Meyer et al. (2014) also noted that the Endozoicomonas was a predominant group in mucus of healthy corals, Porites astreoides, but was significantly less abundant when hosts were injured or diseased. Moreover, relative abundance of Endozoicomonas was altered following increased pCO2 in surrounding seawater or elevated seawater temperature (Morrow et al., 2015; Tout et al., 2015). Furthermore, the predominance of Endozoicomonas might also link to their host's prevalence in nature. In a survey of the distribution of the fungid coral Ctenactis echinata, Roder et al. (2015) reported that when C. echinata was most abundant, the microbiome was highly structured and dominated by Endozoicomonas. In contrast, when C. echinata was less abundant, bacterial communities became diverse and proportion of Endozoicomonas also decreased. A similar phenomenon was also reported by Roterman et al. (2015) regarding the oyster Spondylus in the Mediterranean Sea. In addition, two reports suggested that Endozoicomonas or Endozoicomonas-like bacteria might have a role in sulfur cycling (Raina et al., 2009) and production of antimicrobial compounds (Rua et al., 2014). Taken together, although numerous studies have highlighted the potential importance of Endozoicomonas for corals or reef invertebrates, the genetics, physiology, and ecology of these bacteria are mostly unknown.

Endozoicomonas montiporae CL-33T was the first culturable coral-associated endozoicomonal bacteria, which was isolated from the encrusting pore coral Montipora aequituberculata at the southern coast offshore of Taiwan (Yang et al., 2010). That this species is culturable facilitates in-depth studies, particularly on genetics and physiology. Hence, in this study, we presented a high quality, nearly complete genome sequence of E. montiporae, and deciphered possible ecological functions of this bacterium, based on comparative genomic approaches and physiological tests. Our discoveries provided evidence-based insights regarding how E. montiporae interacted with its host from outside to inside of the coral cell. In addition, based on unique genomic features and newly discovered potential functions, we inferred that E. montiporae can be a facultative endosymbiont with the capability of helping its host to overcome environmental stresses.

Materials and methods

Bacterial strains, medium, and culture conditions

The strain of E. montiporae (CL-33T) was obtained from a previous study (Yang et al., 2010). Type strains of Endozoicomonas elysicola (DSM 22380T), Endozoicomonas numazuensis (DSM 25634T), and Vibrio corallilyticus (DSM 19607T) were purchased from Leibniz-Institut DSMZ—Deutsche Sammlung von Mikroorganismen und Zellkulturen GmbH (Braunschweig, Germany). Routine cultivation was done on MMB medium (per liter: 19.45 g of NaCl, 18.79 g of MgCl2·6H2O, 3.24 g of Na2SO4, 0.55 g of KCl, 0.16 g of NaHCO3, 0.12 g of CaCl2, 2 mg of KNO3, 8 mg of NaHPO4, 5 g of peptone, 1 g of yeast extract, 1 ml of trace element (Guillard and Ryther, 1962), 5.95 g of HEPES, pH 7.2), whereas for high-density cultivation of E. montiporae, 0.1% (w/v) of glucose/maltose was added to MMB. For sole carbon usage test, the MM medium (the same ingredients as MMB except the NaHPO4, peptone, and yeast extract were replaced by 100 μg of cyanocobalamin and 1 ml of 1M K-phosphate buffer) was used and test carbon sources, including glucose and testosterone, were added at 1% (w/v) and 500 μM, respectively. All bacterial strains were aerobically cultured at 25°C with 200 rpm agitation.

Genomic DNA extraction, sequencing, and genome assemble

Genomic DNA of E. montiporae was extracted by the CTAB method, as described (Wilson, 2001). A physical map of E. montiporae genome was constructed from fresh cultures using an optical mapping technology, the ARGUS System (OpGen, Gaithersburg, ML, USA), at Yourgene Bioscience (New Taipei City, Taiwan). Purified genomic DNA was sequenced on three platforms: Illumina HiSeq 2000 (Illumina Inc., San Diego, CA, USA), Roche 454 GS Junior (Roche 454 Life Sciences, Branford, CT, USA) and Pacific Biosciences RS II (Pacific Biosciences, CA, USA). Three DNA libraries were prepared for Illumina sequencing: paired-end library (fragment size ~500 bp), PCR-free library (fragment size ~350 bp), and mate-paired library (fragment size ~4 Kbp). For pyrosequencing (454) and Pacific Biosciences (PacBio) sequencing, one flow cell and one SMRT cell were used for sequencing, respectively. The total read length from high-throughput sequencing was ~10,718 Mbp. Furthermore, CLC Workstation 6 (CLC bio, Aarhus, Denmark) was used to trim low-quality bases off Illumina and 454 raw reads. Illumina and PacBio reads were de novo assembled using ALLPATHS-LG (Butler et al., 2008) and SMRT-Analysis 2.2 (Pacific Biosciences, CA, USA), respectively. Scaffolds from Illumina and PacBio assemblies were analyzed in MapSolver™ (OpGen) to construct chromosomal scaffold by aligning the AflII cutter (CTTAAG) identified in contig sequences to the optical map. Both GMcloser (Kosugi et al., 2015) and GapFiller (Nadalin et al., 2012) were used in gap filling with PacBio and Illumina reads, respectively. When using GMcloser, 454 reads were used to correct PacBio reads with the PBcR pipeline (Berlin et al., 2015) of Celera Assembler 8.2 (Myers et al., 2000) before gap filling. The sequences of rRNA operons were confirmed with molecular cloning method (see below).

Genome annotation and bioinformatic analysis

The assembled genome was first subjected into non-coding gene predictions with tRNAscan-SE 1.3.1 (Lowe and Eddy, 1997), CRISPRFinder (Grissa et al., 2007), and using an HMM search with Rfam database 12 (Nawrocki et al., 2015). After masking predicted non-coding regions, coding genes were predicted with Prodigal-2.6 (Hyatt et al., 2010). Functional assignment for coding genes was accomplished with CDD (Marchler-Bauer et al., 2011), Pfam release 27 (Finn et al., 2014), COG 2014 update (Galperin et al., 2015), and NCBI PSI-BLAST against the Reference Sequence (RefSeq) database with the e-value threshold of 10−6. Metabolic pathways were analyzed with KEGG (Kanehisa and Goto, 2000) and MetaCyc (Karp et al., 2000). Gene annotations of the genome were manually revised and corrected using Artemis (Carver et al., 2011).

Prophage prediction was done with PHAST (Zhou et al., 2011). Repeat sequences were detected with UGENE (Okonechnikov et al., 2012), with a minimum length threshold of 500 bp and minimum identity of 98%. Regions containing rRNA genes (16S and 23S) were excluded from analyses. Repeat sequences were further identified using IS Finder database (https://www-is.biotoul.fr) with the e-value threshold of 10−6. Prokaryotic and eukaryotic subcellular location predictions were performed using Phobius (Käll et al., 2007) and BacelLo (set for prediction in animals; Pierleoni et al., 2006), respectively. The bacterial secretome was predicted using EffectiveT3 (set “gram-” for SignalIP, “type III effector prediction with animal set” for classification module and “selective” for cut-off; Arnold et al., 2009) and SSPred (set to use a Hybrid-II approach; Pundhir and Kumar, 2011). Coiled-coil prediction and homology modeling were performed using COILS (Lupas et al., 1991) and SWISS-MODEL (Biasini et al., 2014), respectively. GlycoEP (Chauhan et al., 2013) was used to predict N-glycosylation sites in eukaryotic proteins. A genome map was generated using Circos (Krzywinski et al., 2009). Visualization of protein models was done with UCSF Chimera v1.10.1 (Pettersen et al., 2004). Reference genomes/transcriptomes used in this study are listed in Table 1.

Table 1

Reference genomes/transcriptomes used for comparative analysis in this study.

Organism name	Taxon (Phylum)	Lineage (Class/Order/Family)	References/Data source
ANIMAL
Acropora digitifera	Cnidaria	Anthozoa/Scleractinia/Acroporidae	Shinzato et al., 2011
Acropora millepora	Cnidaria	Anthozoa/Scleractinia/Acroporidae	Hemmrich and Bosch, 2008
Aiptasia pallida	Cnidaria	Anthozoa/Actiniaria/Actiniidae	Hemmrich and Bosch, 2008
Anemonia viridis	Cnidaria	Anthozoa/Actiniaria/Actiniidae	Hemmrich and Bosch, 2008
Montastrea faveolata	Cnidaria	Anthozoa/Scleractinia/Montastraeidae	Hemmrich and Bosch, 2008
Nematostella vectensis	Cnidaria	Anthozoa/Actiniaria/Edwardsiidae	Hemmrich and Bosch, 2008
Porites astreoides	Cnidaria	Anthozoa/Scleractinia/Poritidae	Hemmrich and Bosch, 2008
Elysia chlorotica	Mollusca	Gastropoda/-/Plakobranchidae	Bhattacharya et al., 2013
Amphimedon quennslandica	Porifera	Demospongiae/Haplosclerida/Niphatidae	Hemmrich and Bosch, 2008
Leucosolenia complicata	Porifera	Calcarea/Leucosolenida/Leucosoleniidae	Hemmrich and Bosch, 2008
Oscarella carmela	Porifera	Homoscleromorpha/Homosclerophorida/Plakinidae	Hemmrich and Bosch, 2008
Sycon ciliatum	Porifera	Calcarea/Leucosolenida/Sycettidae	Hemmrich and Bosch, 2008
BACTERIA
Alcanivorax borkumensis SK2	Proteobacteria	Gammaproteobacteria/Oceanospirillales/Alcanivoracaceae	Schneiker et al., 2006
Alcanivorax dieselolei B5	Proteobacteria	Gammaproteobacteria/Oceanospirillales/Alcanivoracaceae	Lai et al., 2012
Alcanivorax hongdengensis A-11-3	Proteobacteria	Gammaproteobacteria/Oceanospirillales/Alcanivoracaceae	Lai and Shao, 2012
Bermanella marisrubri RED65	Proteobacteria	Gammaproteobacteria/Oceanospirillales/Oceanospirillaceae	Integrated microbial genomes
Chromohalobacter israelensis ATCC 43985	Proteobacteria	Gammaproteobacteria/Oceanospirillales/Halomonadaceae	Integrated microbial genomes
Chromohalobacter salexigens DSM 3043	Proteobacteria	Gammaproteobacteria/Oceanospirillales/Halomonadaceae	Integrated microbial genomes
Endozoicomonas elysicola DSM 22380	Proteobacteria	Gammaproteobacteria/Oceanospirillales/Hahellaceae	Neave et al., 2014
Endozoicomonas numazuensis DSM 25634	Proteobacteria	Gammaproteobacteria/Oceanospirillales/Hahellaceae	Neave et al., 2014
Hahella chejuensis KCTC 2396	Proteobacteria	Gammaproteobacteria/Oceanospirillales/Hahellaceae	Jeong et al., 2005
Hahella ganghwensis DSM 17046	Proteobacteria	Gammaproteobacteria/Oceanospirillales/Hahellaceae	Integrated microbial genomes
Halomonas daqiaogenesis CGMCC 1.9150	Proteobacteria	Gammaproteobacteria/Oceanospirillales/Halomonadaceae	Integrated microbial genomes
Halomonas elongata DSM 2581	Proteobacteria	Gammaproteobacteria/Oceanospirillales/Halomonadaceae	Schwibbert et al., 2011
Halomonas halocynthiae DSM 14573	Proteobacteria	Gammaproteobacteria/Oceanospirillales/Halomonadaceae	Integrated microbial genomes
Kangiella aquimarina DSM 16071	Proteobacteria	Gammaproteobacteria/Oceanospirillales/Alcanivoracaceae	Integrated microbial genomes
Kangiella koreensis DSM 16069	Proteobacteria	Gammaproteobacteria/Oceanospirillales/Alcanivoracaceae	Han et al., 2009
Marinomonas mediterranea MMB 1	Proteobacteria	Gammaproteobacteria/Oceanospirillales/Oceanospirillaceae	Integrated microbial genomes
Marinomonas posidonica IVIA-Po-181	Proteobacteria	Gammaproteobacteria/Oceanospirillales/Oceanospirillaceae	Integrated microbial genomes
Thalassolituus oleivorans MIL-1	Proteobacteria	Gammaproteobacteria/Oceanospirillales/Oceanospirillaceae	Integrated microbial genomes
Vibrio coralliilyticus ATCC BAA-450	Proteobacteria	Gammaproteobacteria/Vibrionales/Vibrionaceae	Kimes et al., 2012
Vibrio harveyi ATCC BAA-1116	Proteobacteria	Gammaproteobacteria/Vibrionales/Vibrionaceae	Integrated microbial genomes
Vibrio shilonii AK1	Proteobacteria	Gammaproteobacteria/Vibrionales/Vibrionaceae	Integrated microbial genomes

Amplification, cloning, and sequencing rrn operons

Five operon-specific forward primers (rrnA_f, rrnB_f, rrnC_f, rrnD_f, and rrnE_f) and one universal reverse primer (rrnUniv_r; Supplementary Table S1) were designed from flanking regions of rrnA-rrnE. Target DNA fragments were amplified by using Phusion® High-Fidelity DNA Polymerase (New England BioLabs, Ipswich, MA, USA) and PCR reactants were prepared according to manufacturer's instruction. The amplification was done with 35 cycles of: 98°C for 10 s, 55°C for 30 s, 72°C for 190 s. Amplified products were separated with agarose gel electrophoresis and DNA bands with expected size (~5.5 Kbp) were eluted by QIAquick Gel Extraction Kit (Qiagen, Hilden, Germany). Purified DNAs were subjected into an A-tailing reaction with Ex Taq DNA Polymerase (TaKaRa, Shiga, Japan) and TA-cloning with pGEM-T vector system (Promega, Fitchburg, WI, USA) according to manufacturer's instruction. Cloned DNAs were transformed into One Shot® TOP10 Chemically Competent E. coli and transformants were screened with β-galactosidase α-complementation. Purification of plasmid DNAs and Sanger sequencing the insert DNAs were carried out at Genomics Inc. (New Taipei City, Taiwan). For complete sequencing the cloned rRNA operons, 16 universal walking primers were designed (Supplementary Table S1).

Steroids degradation and detection

Test strains were aerobically cultured in MMB supplied with 0.1% (w/v) of maltose and 500 μM of testosterone. On each day, 1 ml of culture sample was taken, vigorously mixed with 0.8 ml of ethyl acetate and centrifuged (10,000 × g, 4°C, 10 min). After centrifugation, 0.6 ml of the upper-most layer was withdrawn and vacuum-dried. For detection of testosterone-derived compounds, dried samples were dissolved in 40 μl of ethyl acetate and 5 μl applied on TLC Silica gel 60 F254 plates (MERCK, Whitehouse Station, NJ, USA) and developed in dichloromethane:ethyl acetate:ethanol (7:2:0.5, v/v; Chiang et al., 2010). The separated metabolites were visualized by spraying on plates with 40% sulfuric acid solution, followed by heating at 150°C.

Nucleotide sequence accession number

The genome sequence of E. montiporae and its annotation were deposited at GenBank under BioProject Number PRJNA66389 (Registration date: 28-Apr-2011).

Results and discussion

General genomic features

The extracted genomic DNA of E. montiporae was sequenced though five libraries in three next-generation platforms and reads were assembled into one circular chromosome ~5.43 Mbp and ~1900-fold in coverage depth. Two gaps (~9.8 and 1.0 Kbp) remained (genome completeness was approximately 99.80%). At the beginning, lengths of assembled contigs from Illumina/454 reads were highly fragmented and in silico AflII patterns of contigs were not properly matched to the physical map generated by optical mapping technology. We also noticed that a previously released E. montiporae genome (Neave et al., 2014) had the same limitation, which was attributed to sequence assemblies being greatly affected by numerous repeat sequences (which are longer than sequencing reads). Therefore, PacBio long-read sequencing technology was used. After incorporating long-reads, the final assembly had good agreement in AflII cutting pattern (indicates the genome was accurately assembled; Figure 1). Detailed genome statistics are listed (Table 2).

Figure 1

Physical and genetic maps of the genome of . The track numbers represent: (1), IS linkages; (2), pseudogenes; (3), rRNA operons; (4), tRNAs; (5), eukaryotic domain proteins; (6), testosterone degrading gene cluster; (7), prophage regions; (8), unique genes in E. montiporae; (9), conserved genes in three Endozoicomonas; (10), genes shared with E. numazuensis; (11), genes shared with E. elysicola; (12), GC profile (red, > 50%; blue, < 50%); (13), ALLPATH-LG assembled scaffolds (green, orientated by MpSolver; light blue, orientated in gap filling process); (14), SMRTAnalysis assembled scaffolds (yellow, orientated by MpSolver; light blue, orientated in gap filling process); (15), AflII pattern from assembled sequences; (16) physical map (AflII cuts) generated with optical mapping technology. The black lines that connect track 15 and 16 indicated the alignment of AflII cuts from assemble genome sequence and physical map. The two assembly gaps are indicated by black near position 1.6 and 2.8 Mbp.

Table 2

Generic features of the .

Parameter	E. montiporae
Genome size (estimated from optical mapping)	5,432,010 bp
Total size of genome assembly	5,430,252 bp
GC content	48.37%
Number of coding sequences	4572
Number of rRNA operons	7
5S rRNA gene	8
16S rRNA gene	7
23S rRNA gene	7
Number of tRNAs	109
Number of non-coding genes	55
Number of pseudogenes/frame-shifted genes	198
Number of CPISPR arrays	5

The genome contained seven rRNA operons. The 16S-23S rRNA gene copy number in E. montiporae differed from that in other reported endozoicomonal genomes (16S/23S: E. elysicola, 6/6; E. numazuensis, 5/2; E. montiporae, 4/4; Neave et al., 2014). In contrast, the copy number of 5S rRNA gene was 8 (E. elysicola and E. montiporae, 8; E. numazuensis, 3; Neave et al., 2014). The rRNA operon rrnD was the only operon that contained two copies of 5S rRNA genes (Supplementary Information). In the orthology analysis, E. montiporae shared ~59.49 and ~63.09% coding genes with E. elysiola and E. numazuensis genomes, respectively. In the metabolic pathway analysis, the three species had the potential to synthesize all proteinogenic amino acids and many cofactors, prosthetic group and electron carriers required for growth, except vitamin B12.

In the comparative analysis, there were 1179 coding genes that could only be detected in E. montiporae (Figure 1). Based on the COG functional profile, E. montiporae had a higher proportion of genes (2.5% of all COGs) related to the mobilome (Figure 2). Many genes within the mobilome were similar to phage structural proteins (i.e., phage capsid or tail proteins), indicating frequent interactions between E. montiporae and its infectious agents, bacteriophages. Based on prophage prediction, there were eight prophages in the genome (ranging from 15.5 to 75 Kbp; Table 3). Four of those prophages were predicted as intact lysogenic phages, similar to Vibrio phage vB VpaM MAR (region 1 in the genome), Vibrio phage VPUSM 8 (regions 2 and 3), and Escherichia phage vB_EcoM-ep3 (region 4). In contrast, no prophage was detected in the draft genome of E. elysicola and E. numazuensis.

Figure 2

Functional COG profile of the . Genes shared in three Endozoicomonas, with E. numazuensis, with E. elysicola, or unique in E. montiporae are highlighted as gray, brown, green, and red with respectively. Corresponding locations are shown in Figure 1 (track 8–11).

Table 3

Prophage regions identified in the genome of .

Region	Position	Length (Kbp)	Completeness	PHAST score	ORF number	Unique ORFs	Possible phage
1	1624494–1699528	75	Intact	150	99	69	Vibrio phage vB VpaM MAR (NC_109722)
2	2025905–2046887	20.9	Intact	140	25	25	Vibrio phage VPUSM 8 (NC_022747)
3	2717260–2741388	24.1	Intact	140	32	31	Vibrio phage VPUSM 8 (NC_022747)
4	4621590–4662332	40.7	Intact	150	67	62	Escherichia phage vB_EcoM-ep3 (NC_025430)
5	1858945–1891753	32.8	Incomplete	60	31	13	Staphylococcus prophage phiPV83 (NC_002486)
6	2179474–2203283	23.8	Incomplete	40	21	10	Stenotrophomonas phage S1 (NC_011589)
7	2774561–2790090	15.5	Incomplete	50	28	24	Pseudoalteromonas phage PM2 (NC_000867)
8	4748489–4767285	18.7	Incomplete	30	27	21	Rhizobium phage RR1-B (NC_021557)

The E. montiporae genome contained at least 452 mobile elements, which include 46 IS1-like IS elements (most active), 36 of Group II introns and 370 of integrase/transposase coding genes from various IS families. The estimated repeat coverage in E. montiporae was ~9.07%, greatly higher than E. elysicola and E. numazuensis, which were only ~0.11 and ~4.05%, respectively (Figure 1, also see Supplementary Figure S1). Those mobile elements were likely still active in E. montiporae, because 56 genes were inactivated due to insertions (Supplementary Table S2).

Endozoicomonas montiporae is a testosterone degrader

The three endozoicomonal species in the present study shared high similarities in metabolic capacities. It was noteworthy that they all had genes for degrading testosterone (male sex hormone), although that was rarely reported in oceanospirial genomes. Key genes required for degrading testosterone were detected within five of seven conserved gene clusters (I, II, V, VI, and VII) in the endozoicomonal genomes (Figure 3A, Supplementary Table S3). From a reconstructed pathway, testosterone could be degraded into propionyl-CoA and pyruvate through the 9,10-seco pathway (Coulter and Talalay, 1968). Furthermore, the first two metabolites could be further utilized (via propionate metabolism and TCA cycle) in the Endozoicomonas species (Figure 3B). Clusters III and IV were mainly composed of genes related to fatty acid β-oxidation (Supplementary Table S3). Interestingly, those genes formed conserved synteny with known aerobic testosterone-degrading bacteria (Horinouchi et al., 2012), indicating those genes were probably acquired via horizontal gene transfer. In our degradation tests, all three Endozoicomonas were able to completely degrade testosterone (Figures 3C–E, with Endozoicomonas; Figure 3G, medium only control). In contrast, there were no genes related to testosterone degradation identified in Vibrio species and the coral pathogen V. corallilyticus had a negative reaction regarding testosterone degradation (Figure 3F). When minimal medium was used to determine whether testosterone could serve as sole carbon source for growth, only E. elysicola and E. numazuensis, but not E. montiporae, grew (albeit with a very low biomass; i.e., OD600 ~0.1). Low testosterone concentrations can be detected in coral reef seawater, with highest concentrations in coral tissues (42.7 ng/g tissue in conjugated form; Twan et al., 2006). However, the concentration was still ~5,300-fold lower than our tests, suggesting that testosterone might not be a nutrient but more like an “animal cue” for Endozoicomonas.

Figure 3

Testosterone degradation in . The gene order of key genes required for degradation (A) and their corresponding reactions (B) presented in Endozoicomonas species. Degrading capabilities of test bacteria were confirmed by supplying 500 μM of testosterone in culture medium and sampled daily in the previous culture. Test bacteria are: E. montiporae (C), E. elysicola (D), E. numazuensis (E), V. corallilyticus (F), and a blank control without inoculum (G). Abbreviations used in TLC analysis are: 1-DT, 1-dehydrotestosterone; AD, androst-4-ene-3,17-dione; ADD, androsta-1,4-diene-3,17-dione; S, standard mixture; SF, solvent front; T, testosterone. Numbers stand for the day of culture.

Endozoicomonas montiporae possesses a N-deglycosylation enzyme that might help it to penetrate host's mucus layer

EZMO1_4380 (endo-AEmo) was a unique gene in E. montiporae. Its translated product was similar to endo-β-N-acetylglucosaminidase A (PDB number: 2VTF) from Arthrobacter protophormiae (Supplementary Figure S2). Endo-β-N-acetylglucosaminidases (ENGases) have been reported in a wide variety of organisms, with roles in hydrolyzing glycosidic bond present in N-linked sugar chains in glycoproteins (Yin et al., 2009). In contrast to those cytoplasmic ENGases in plants and animals, bacterial and fungal ENGases are secreted enzymes (Suzuki et al., 2002). The Endo-AEmo contained a signal peptide sequence at its N-terminus (1–25 residues), a hallmark feature of a secreted enzyme. The predicted function of this protein might be related to digest glycoproteins, which were commonly present in the natural habitats of E. montiporae, e.g., coral mucus. Glycoprotein mucins are the major structural component of coral mucus. However, only the O-linkage was experimentally detected in coral mucins (Meikle et al., 1987). In genomic surveys, the starlet sea anemone Nematostella vectensis has an N-glycosylated mucin coding gene (XP_001636262; Lang et al., 2007). In addition, in the draft genome of the coral Acropora digitifera, open reading frames were detected in two regions located inside the scaffold 125 and 132 that resembled the mucin protein of N. vectensis and contained putative N-glycosylation sites. Therefore, we inferred that coral mucins could be N-glycosylated and physical characteristics might be affected due to deglycosylation by Endo-AEmo. Blocking N-glycosylation in rat Muc2 reduced gel formation and/or mucin viscoelasticity (Bell et al., 2003). Hence, relaxation of mucin by Endo-AEmo would facilitate E. montiporae to pass through coral mucus. Actually bacteria degrade mucin by two types of reactions, namely O-linkage hydrolysis and proteolysis (McGuckin et al., 2011). However, no similar homologs related to enzymes of the two reactions were detected in Endozoicomonas genomes, suggesting the bacteria might not be able to deeply decompose their host mucus. Perhaps the role of Endo-AEmo was not solely related to mucus dissociation, but it could help the bacterium to reach to specific host receptors on the cell membrane (Figure 4), similar to other bacteria (McGuckin et al., 2011).

Figure 4

Proposed infection model of the . When E. montiporae attaches to the host mucus layer, secreted Endo-AEMO could reduce mucus viscosity, enabling the bacterium to penetrate through the mucus (1). Once the bacterium reached the host plasma membrane, the EfnB2_2 could bind with host Eph receptors (2), and initiate endocytosis (3). When E. montiporae enters into a host cell, it might use extracellular enzymes to disrupt the endosome and interfere with phagolysosome maturation by chelating Rabs with EZMO1_3398 (4). After escaping the endosome, E. montiporae could utilize nutrients in the host cytoplasm, particularly glucose produced by endosymbiotic zooxanthellae (5).

Host's ephrin receptors might be involved in endocytosis of E. montiporae

Two genes, EZMO1_1051 (efnB2_1) and EZMO1_3739 (efnB2_2), encoded proteins similar to eukaryotic proteins, but lacking orthologs in other bacteria. In that regard, EfnB2_1 and EfnB2_2 corresponded to NEMVEDRAFT_v1g212995 and NEMVEDRAFT_v1g236436 from N. vectensis, with 39 and 34% amino acid identities, respectively. Both EfnB2_1 and EfnB2_2 contained ephrin ectodomain (cd02675) and secretion signals at their N-terminus. Besides, EfnB2_2 encompassed transmembrane regions and was predicted as a membrane protein that could adhere on a cell surface, whereas EfnB2_1 was not. Tertiary structures of the two proteins, predicted by homology modeling, were very similar to each other, and they also resembled type B ephrin-B2 from mouse (1IKO; Supplementary Figure S3). However, it is noteworthy that the ephrin ligand protein was similar to the cupredoxin domain, a protein domain containing copper binding sites and usually responsible for electron-transfer reactions in bacteria (De Rienzo et al., 2004). Distinct from cupredox domain proteins, ephrin acts as a signal molecule in animals, with no metal-binding capability (Toth et al., 2001). For the two proteins, EfnB2_1 and EfnB2_2, all lacked copper binding residues by aligning with bacterial cupredoxins (azurins and auracyanins), indicating that the functions of EfnB2_1 and EfnB2_2 were more similar to eukaryotic ephrin ligand.

In vertebrates, ephrins are ligands for ephrin receptors (Ephs) and binding can activate intracellular signal pathways involved in a variety of cellular functions, including differentiation, migration, segmentation, and endocytosis (Kullander and Klein, 2002; Pitulescu and Adams, 2010). In cnidarian, ephrins, and Ephs in Hydra vulgaris belonged to type B (Tischer et al., 2013). Furthermore, ephrin ligands and receptors were strongly expressed on endoderm of tentacles and buds, suggesting ephrin-Eph signaling might have roles in cell adhesion and tissue boundary formation in H. vularis (Tischer et al., 2013). Ephrin-Eph signaling in corals has never been discussed; however, ephrin/Eph receptor genes were present in the genome of A. digitifera and the ephrin ligand gene was differentially expressed in A. palmata larvae (Polato et al., 2013), indicating that ephrin-Eph signaling could have roles in corals.

Although, the reason why E. montiporae has ephrin genes remains unknown, the rare occurrence of eukaryotic cell surface ligand coding genes might have a role in targeting host receptor and initiating internalization. In contrast, pathogenic bacteria can use azurin, a cupredoxin protein similar to ephrin in structure, to invade host cells. For example, bacterial azurin can solely enter the murine macrophage J774, as well as several types of cancer cells and trigger programed cell death of mammalian cell lines by interfering with Eph signal pathways (Yamada et al., 2005; Chaudhari et al., 2007). Furthermore, azurin in the human pathogen Neisseria gonorrhoeae is essential for survival inside host cells (Wu et al., 2005). Azurin othologous genes were only detected in three pathogenic Vibrio genomes, but not in the three Endozoicomonas. Accordingly, we speculated that the ephrin ligand in E. montiporae could serve as a “cloak” and facilitate this bacterium to enter host cells, but would not impair its host (Figure 4).

A giant secreted protein of E. montiporae might modulate host's intracellular vesicle trafficking

To survive inside a host cell, intracellular pathogens and symbionts have to evade the host defense system, which relies on intracellular vesicle trafficking (Davy et al., 2012; Ashida et al., 2015). In the present study, it was determined that E. montiporae contained at least one gene that may be essential for this purpose. In that regard, EZMO1_3398 is a giant extracellular protein-coding gene, which encodes 3872 amino acids and comprises the myosin tail (pfam01576) and Mycoplasma helix-rich (TIGR04523) domains at its N-terminus (Supplementary Figure S4A). Furthermore, EZMO1_3398 was predicted as a Sec-dependent secreted protein but not a membrane protein. The myosin tail domain contains α-helices that can form a coiled-coil structure for dimerization (Krendel and Mooseker, 2005). Regarding the Mycoplasma helix-rich domain, it is still functionally unclear, although it is common in Mycoplasma species.

Based on structure prediction, EZMO1_3398 could be a chimeric protein, encompassing bacterial and eukaryotic protein features. Potential functions could involve modulating host's intracellular trafficking processes, particularly during initial stages of infection (Figure 4). Three protein substructures were detected in the EZMO1_3398, including the tail of myosin (1I84), LidA (3TNF), and stalk with microtubule-binding domain from dynein 2 (4RH7), at its first half (Supplementary Figure S4A). The LidA-like structure was a core domain that connected the myosin tail-like and the dynein 2 stalk-like structures to its helix-1 and helix-7, respectively (Figures S4C–E). However, no myosin/dynein motor domains were identified, neither in sequence nor in structure similarities, in EZMO1_3398. The LidA protein of Legionella penumophila is a Rab supereffector protein that strongly binds to host's Rab proteins and might help L. penumophila to survive inside a host cell by preventing autophagosome maturation (Joshi and Swanson, 2011; Schoebel et al., 2011). Interfering with Rab proteins inside hosts seems to be a common strategy also used by endosymbionts to establish a symbiotic relationship. Even though a detailed mechanism is unclear, study of Symbiodinium-Aiptasia pulchella symbiosis suggests the symbiotic zooxanthellae could manipulate endosomal trafficking and prevent lysosome fusion (Davy et al., 2012). Rab proteins can perform several cellular functions when binding to various effector proteins. Regarding vesicle movement, Rabs can recruit myosins or dyneins to exert movement along microtubes. Binding with effectors and GTPase activating protein can activate the GTPase activity of Rab and hydrolysis of GTP can trigger movement (Hutagalung and Novick, 2011). Interestingly, EZMO1_3398 contained substructures similar to myosin and dynein and both are connected with the Rab-binding central domain. Therefore, the EZMO1_3398 might be able to attach host's microtubes via binding certain Rab proteins for some functional purposes, e.g., to prevent lysosome fusion (Figure 4).

Type III secretion effector proteins might help E. montiporae to infect coral host and mediate host's metabolism

A common feature of the Endozoicomonas species in this study is that they all had a host-dependent secretion system, the type III secretion system (T3SS), a well-known bacterial secretion system that is functionally involved in bacteria–host interactions, either symbiosis or pathogenesis (Preston, 2007). The presence of T3SS in the Endozoicomonas suggests they are capable of interacting with their hosts. To identify potential interactions between E. montiporae and its coral host, we analyzed predicted T3SS secretomes and the possible subcellular locations of those secreted proteins (i.e., effectors). There were 242 proteins predicted to be T3SS effectors, including many hypothetical proteins (Supplementary Table S4). Among those possible T3SS effectors, five proteins might be involved in survival inside hosts, regulating host's metabolism and/or increasing host's fitness.

EZMO1_0953 is a catalase gene of E. montiporae and its product was one of the T3SS effectors. Besides, catalase was a common T3SS-screted protein in other Endozoicomonas. Secretion of catalase by a host-dependent system indicated that the H2O2-scavenging enzyme might be crucial for Endozoicomonas to live within their hosts. For animal-associated bacteria, including pathogens, parasites, or symbionts, catalase activity is important for them to survive in their hosts (Bishai et al., 1994; Rocha et al., 1996; Visick and Ruby, 1998). Interestingly, instead of catalase, we identified another type of H2O2-removing enzyme, the thiol peroxidases (Tpx), that could be secreted via T3SS in pathogenic bacteria Vibrio. Tpx is a part of T3SS assembly in the human pathogen Yersinia species and also connects to effector delivery (Wang et al., 2011). In Shiga toxin-producing Escherichia coli O157:H7, Tpx is required for biofilm formation on the human HT-29 cell line (Kim et al., 2006). The same gene was detected and predicted as a T3SS-secreted protein in our analysis. Collectively, we speculated that the T3SS-secreted antioxidant profile could be a distinguishable signature between pathogenic and non-pathogenic bacteria.

EZMO1_3421 encodes isocitrate lyase (ICLEmo), a key enzyme of the glyoxylate cycle, and was a predicted T3SS secreted protein in E. montiporae. The glyoxylate cycle is common in bacteria, fungi and plants, and can bypass carbon flow from fatty acid degradation to gluconeogenesis. Glyoxylate cycle enzymes are also present in nematodes and cnidarians (Kondrashov et al., 2006) and expression of ICL was detected in transcriptomes of A. millepora larvae and A. palmate (Meyer et al., 2009; Polato et al., 2013).

In addition to ICLEmo, there were two other T3SS effectors in E. montiporae that could concurrently work with the host enzymes for gluconeogenesis, namely EZMO1_2385 and EZMO1_2390 that encode the aconitate hydratase (AcoNEmo) and phosphoenolpyruvate synthase (PpsAEmo), respectively. In bacteria, AcoN and PpsA can redirect metabolites from the TCA cycle to gluconeogenesis. Up-regulation of fatty acid β-oxidation, glyoxylate cycle and gluconeogenesis have been proposed to enable corals to convert stored fatty acids to carbohydrates for surviving stress-induced starvation (Kenkel et al., 2013). Hence, coupling of T3SS secreted ICLEmo, AcoNEmo, and PpsAEmo from E. montiporae might increase metabolic efficiency of gluconeogenesis in the coral host. Furthermore, such metabolism integration might promote host survival under stressful conditions (Figure 5).

Figure 5

Proposed pathway modulation model of . When the coral host is under stress, loss of symbiotic zooxanthellae (indicated by gray) makes host cell starts to use fatty acids and convert them into glucose. Active energy production and metabolite exchange could lead to accumulation of reactive oxygen species (ROS) in mitochondria matrix, which can oxidize deoxyribose purine triphosphates (dRTP) into 8-oxo-dRTP and cause mutation. The T3SS secreted enzymes could hasten carbon flow from fatty acids to glucose and hydrolyze 8-oxo-dRTP inside mitochondria. The abbreviations of host enzymes are: ACO, aconitate hydratase; MS, malate synthase; ME1, NADP-dependent malic enzyme; MDH2, malate dehydrogenase (cytoplasmic).

Interestingly, we identified ICL orthologs that did not belong to T3SS effectors in E. elysicola and E. numazuensis. Notably, glyoxylate cycle enzymes were not detected in the representative genomes for their hosts (i.e., sea slug and sponge; Table 4). Besides, the predicted T3SS secretomes of Endozoicomonas were less common (Supplementary Table S4). Therefore, Endozoicomonas might have evolved unique strategies to interact with their own specific hosts.

Table 4

Glyoxylate cycle enzymes (ICL and MS) in selected marine invertebrates.

Type	Organism name	Data source	Number of BLAST hit
			ICL	MS
Coral	Acropora digitifera	Genome	1	2
	Acropora millepora	ESTs	1	2
	Montastrea faveolata	ESTs	3	1
	Porites astreoides	ESTs	2	2
Sea anemone	Aiptasia pallida	ESTs	1	2
	Anemonia viridis	ESTs	1	0
	Nematostella vectensis	Genome	1	2
Sea slug	Elysia chlorotica	SRA	0	0
Sponge	Amphimedon quennslandica	Genome	0	0
	Leucosolenia complicata	Genome	11	2
	Oscarella carmela	Genome	0	0
	Sycon ciliatum	Genome	2	0

1. ICL and MS sequences from Saccharomyces cerevisiae were used as representitive sequences (queries).

2. When searching in EST database, tblastn was used and thresholds were set: e-value < 1e-15; bitscore > 100.

3. When searching in SRA database, tblastn was used first. Reads may affiliate with query sequences were recovered and subjected into blastn search against to nr/nt database. In blastn search, reads shared high identities and coverages to bacterial subjects will be excluded.

The EZMO1_3450 (mthEmo) encodes the 7,8-dihydro-8-oxoguanine triphosphatase (MTH) and was another T3SS effector. The predicted subcellular location in host cell of MTHEmo was in the mitochondria. The MTH enzyme can hydrolyze damaged purine nucleoside triphosphates caused by attacks from reactive oxygen species and therefore can confer protection to the cell against various oxidative stresses. In animal models, MTH can prevent mitochondrial dysfunction, attenuate stress-induced cell death and enhance vitality (Ichikawa et al., 2008; De Luca et al., 2013). Even though the detailed mechanism is still unknown, bacterial endosymbiosis can prevent mitochondrial dysfunction in a mycorrhizal host and improve fitness of its fungal host during the pre-symbiotic rhizospheric phases (Salvioli et al., 2015). Exporting MTHEmo from E. montiporae to the host's mitochondria might promote mitochondrial functions. In addition, because mitochondria are also involved in fatty acid oxidation, MTHEmo could promote conversion of fatty acids to glucose in coral cells under bleaching stress (Figure 5). Accordingly, E. montiporae might be able to enhance host's fitness to environment changes by protecting mitochondria from the oxidative injury and/or altering carbon flows of host metabolism.

T3SS effectors might also involve in modulating host's signaling pathways

Two eukaryotic signal pathway proteins, serine/threonine protein kinase (STPK; EZMO1_1618 and EZMO1_3446), were predicted as T3SS effectors in E. montiporae. Both EZMO1_1618 and EZMO1_3446 were similar to STPK genes from the protozoa Trypanosoma congolense (TCIL3000_11_7730; identities: 33%) and the fungus Phytophthora nicotianae (L916_19948; identities: 34%), respectively. The two proteins were unlikely to be membrane-bound receptors, due to a lack of transmembrane regions in their amino acid sequences. Furthermore, predicted subcellular locations of EZMO1_1618 and EZMO1_3446 were all in the host nucleus (Figure 4). Bacterial types of STPKs can have roles during infection of their hosts (Whitmore and Lamont, 2012; Canova and Molle, 2014). For example, the human pathogen Yersinia species can secrete STPK protein (i.e., YopO) into host cells via T3SS and can interfere with the actin cytoskeleton of macrophage to prevent phagocytosis (Black et al., 2000; Juris et al., 2000; Grosdent et al., 2002). In contrast, the T3SS-secreted STPKs in E. montiporae were not similar to bacterial counterparts regarding protein orthology and predicted subcellular locations. Perhaps the two STPKs in E. montiporae were involved in altering certain gene expression in their host's nucleus.

Endozoicomonas montiporae could be a facultative coral endosymbiont

Based on niche specificity, special genomic characteristics, and growth features, we inferred that E. montiporae could be a facultative endosymbiotic bacterium in corals.

One of the indirect albeit rational evidences for niche specificity is that Endozoicomonas were mostly detected or isolated from tissues of marine invertebrates in reefs (Kurahashi and Yokota, 2007; Yang et al., 2010; Bourne et al., 2013; Nishijima et al., 2013; Pike et al., 2013; Hyun et al., 2014), implying their niches were likely associated with reef marine invertebrates. Moreover, the Endozoicomonas strains in this study all had genes of the host-dependent secretion system, carrying several putative animal genes and were able to degrade testosterone, strongly suggesting that the natural niche of those bacteria should be spatially close to or indeed within their animal hosts. We further deduced that the bacterial niches being intimate with host cells might be an important factor for gene transfer from hosts to bacteria.

Highly represented repeat sequences and pseudogenes in a bacterial genome could be a sign of genome erosion, a common feature in host-restricted symbionts or pathogens (McCutcheon and Moran, 2011). Occurrence frequencies of the IS and pseudogenes among the three endozoicomonal genomes could be grouped into various stages of genome erosion: E. elysicola, free-living; E. numazuensis, intermediate between free-living and host-restricted symbiont; and E. montiporae, a recently host-restricted symbiont/pathogen (Supplementary Figure S1) according to the category proposed by McCutcheon and Moran (2011). Such genomic features have been reported from an E. elysicola-like species, which may cause epitheliocystis in sharpsnout seabream larva (Katharios et al., 2015). Even though the genetic basis which can cause fish diseases remains unclear, increasing genome plasticity by accumulating a high proportion of insertion sequences or improving fitness to hosts by eliminating redundant genes through a pseudogenization process seemed to be a common strategy to survive inside host tissue. Furthermore, this was consistent with an evolutionary change toward tightly host-associated lifestyles that could occur in certain Endozoicomonas species.

Ecological roles of bacteria for mutualism or parasitism are commonly dynamic. Based on genomic analysis, we inferred that E. montiporae could be considered as a bacterium beneficial to corals. Not only does this bacterium have genes that could help its hosts, it was noteworthy that no known bacterial toxin homologs were present in its genome. However, potential parasitic roles of other Endozoicomonas species could not be excluded. For example, there are two reports that E. elysicola-like species could be associated with epitheliocystis in fish larvae (Mendoza et al., 2013; Katharios et al., 2015). However, these observations were different from previous surveys of marine invertebrates. Due to limited knowledge, it is unknown whether such different life styles of Endozoicomonas are due to genetic variations (e.g., presence or absence of virulence genes), physiological variations of their hosts, or both.

Some enzymes of E. montiporae could be secreted into host cytoplasm and promote glucose production in host cells. Glucose could be a critical factor in host cells and a feedback from the host cell, which can support growth of E. montiporae. In the present study, compared to other strains, E. montiporae did not achieve high-density growth in rich medium without glucose (Supplementary Figure S5). Therefore, supplying glucose could be a key to maintain a viable E. montiporae population inside corals.

To survive inside corals, it would be necessary for E. montiporae to communicate with its coral hosts as well as with endosymbiotic zooxanthellae, a key symbiotic partner inside coral cells. Glucose supply might be essential to maintain a symbiotic relationship between E. montiporae and coral. Therefore, it is also important to detect interactions between bacterial and algal endosymbionts. Therefore, we searched for potential proteins in the bacterium that could be involved in interactions, although none was detected. Unfortunately, our bioinformatics analyses lacked suitable models for prediction. In that regard, current prediction models of subcellular locations have not been established from corals or animals with algal endosymbionts. Therefore, for example, a protein could be translocated to specific organelles or compartments in which zooxanthellae reside, i.e., the symbiosome (Trench, 1979; Birkeland, 1997). It was noteworthy that eukaryotic domain protein coding genes present in E. montiporae, especially STPKs, were similar to proteins in invertebrates or fungi, but not plants. Therefore, we concluded that exchanges of genetic material between E. montiporae and zooxanthellae was unlikely. Although, our current results suggested the two endosymbionts might have no direct contact, we cannot assume that the E. montiporae does not communicate with algal symbionts, for example, by metabolites (Glick, 2012).

In conclusion, bacterial consortiums within corals are extremely complex and there has been a lack of detailed knowledge regarding interactions between bacteria and their coral host. Therefore, studying interactions between one bacterium and one coral species may provide more direct evidence and details. In this research, we focused on the Endozoicomonas, a group of bacteria present mainly on healthy coral, and how they may assist their hosts. We provided an abundance of genomic information, and discussed how Endozoicomonas, particularly E. montiporae could interact with their hosts. Our findings provided a valuable guide for future in-depth molecular or physiological studies for coral microbiology.

Author contributions

ST, WC, and JS conceived of the work. JD and JS conducted whole genome sequencing and were involved in genome assemble. JC contributed to comparative genome analysis, functional gene prediction, and other bioinformatics analysis. YC and JD designed the experiment of testosterone degradation and was carried out by JD. JD and ST contributed to writing the manuscript and JD elaborated the figures and tables. All authors critically reviewed, revised and ultimately approved this final version.

Conflict of interest statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.