Whole-Genome Transcriptome Analyses of Native Symbionts Reveal Host Coral Genomic Novelties for Establishing Coral-Algae Symbioses.

Reef-building corals and photosynthetic, endosymbiotic algae of the family Symbiodiniaceae establish mutualistic relationships that are fundamental to coral biology, enabling coral reefs to support a vast diversity of marine species. Although numerous types of Symbiodiniaceae occur in coral reef environments, Acropora corals select specific types in early life stages. In order to study molecular mechanisms of coral-algal symbioses occurring in nature, we performed whole-genome transcriptomic analyses of Acropora tenuis larvae inoculated with Symbiodinium microadriaticum strains isolated from an Acropora recruit. In order to identify genes specifically involved in symbioses with native symbionts in early life stages, we also investigated transcriptomic responses of Acropora larvae exposed to closely related, nonsymbiotic, and occasionally symbiotic Symbiodinium strains. We found that the number of differentially expressed genes was largest when larvae acquired native symbionts. Repertoires of differentially expressed genes indicated that corals reduced amino acid, sugar, and lipid metabolism, such that metabolic enzymes performing these functions were derived primarily from S. microadriaticum rather than from A. tenuis. Upregulated gene expression of transporters for those metabolites occurred only when coral larvae acquired their natural symbionts, suggesting active utilization of native symbionts by host corals. We also discovered that in Acropora, genes for sugar and amino acid transporters, prosaposin-like, and Notch ligand-like, were upregulated only in response to native symbionts, and included tandemly duplicated genes. Gene duplications in coral genomes may have been essential to establish genomic novelties for coral-algae symbiosis.

Significance

In the past, to study molecular mechanisms of symbiosis between corals and algae (dinoflagellates), transcriptomic analyses of coral larval inoculations with various algae have been conducted; however, those studies used algal strains that were not isolated from corals. In this study, we performed whole-genome transcriptome analyses of coral larvae inoculated with native algal strains (Symbiodinium) isolated from the coral, and compared the transcriptomic responses of corals with two closely related, non- and occasionally symbiotic algal strains. We identified genes specifically involved in symbiosis with native symbionts, and found that active utilization of symbionts by host corals occurred exclusively when native symbionts were inoculated. We also found that some of these genes originated by gene duplication in the coral genome.

Introduction

Coral reef ecosystems, built on a foundation of reef-building corals (anthozoan cnidarians), harbor ∼30% of all marine species and are regarded as the most diverse marine ecosystems on earth (Roberts 2002; Fisher et al. 2015). Photosynthetic algae (Symbiodiniaceae) live inside host coral cells (endosymbionts) in a mutualistic relationship. These algae are believed to provide up to 90% of the nutrients required by the host corals (Muscatine and Porter 1977; Muscatine et al. 1981; Falkowski et al. 1984; Yellowlees et al. 2008; Melo Clavijo et al. 2018), and host corals, in turn, provide shelter and CO2 for the symbiotic algae. These mutualistic relationships are fundamental to maintenance of coral reefs. Collapsed symbioses induced by environmental stressors, such as increased seawater temperature or ocean acidification, often result in coral bleaching (Brown 1997), which ultimately leads to loss of habitat for numerous marine species and can collapse entire coral reef ecosystems (Bellwood et al. 2004).

The family Symbiodiniaceae is divided into nine clades (e.g., Clade A-I), seven of which are formally described as genera (LaJeunesse et al. 2018). Each clade or genus can be further subdivided genetically into numerous types or species (LaJeunesse 2005). Although Acropora, the most common and widespread coral genus in the Indo-Pacific (Wallace 1999), generally harbors algae of the genus Cladocopium (formerly Clade C) in its adult stage (LaJeunesse 2005; Lien et al. 2012), Acropora recruits often harbor other genera, especially Symbiodinium (formerly Clade A) or Durusdinium (formerly Clade D) (Abrego et al. 2009; Yamashita et al. 2013, 2014). Moreover, Symbiodinium microadriaticum (Clade A, ITS2 type A1) was predominant among Symbiodinium types in Acropora recruits, although numerous types of Symbiodiniaceae, including Cladocopium were commonly observed in coral reef sea water at Ishigaki Island, Okinawa, Japan (Yamashita et al. 2014), indicating that S. microadriaticum is a native symbiont in early life stages of Acropora in Okinawa. The genus Symbiodinium is highly diverse and includes nonsymbiotic and symbiotic types or species (Yamashita and Koike 2013). A nonsymbiotic type, S. natans (Clade A, ITS2 type A2-relative) cells can be artificially introduced to Acropora tenuis larvae under high cell densities that are unlikely to occur in nature, but these are not maintained by the corals for very long (Yamashita et al. 2018). Thus, mechanisms for acquisition and maintenance of only appropriate symbiont types in host corals must exist, but molecular mechanisms enabling this specificity are unknown.

In order to study molecular mechanisms of symbiosis between corals and algae, various transcriptomic analyses of coral larval inoculation with Symbiodiniaceae have been conducted (Meyer and Weis 2012; Mohamed et al. 2016) however, those studies used Symbiodiniaceae strains that were not isolated from corals, possibly due to difficulties of culturing symbiotic strains. Recently, two studies using Symbiodiniaceae strains isolated from anthozoans for inoculation experiments were reported. One inoculated A. tenuis larvae with algal strains isolated from other coral species or from anemones (Yuyama et al. 2018), and another used corals (corallimorpharians) artificially bleached with menthol to inoculate strains isolated from the same or other coral species (Lin et al. 2019). Cultured strains of major symbionts from adult Acropora in Okinawa, Cladocopium, are not available to date. Furthermore, no methods exist for removing symbiont algae completely from adult corals without stress, for maintaining adult corals without symbionts, or for reintroducing arbitrary algae artificially into adult corals. Thus, previous studies may not reflect natural conditions of scleractinian corals in coral reefs.

On the other hand, Acropora larvae intrinsically have no symbiotic algae, and a method for artificially introducing native and nonnative symbiont algae into coral larvae has been established (Yamashita et al. 2018). In this study, in order to clarify molecular mechanisms of symbioses occurring in nature, we sequenced the whole genome of a native symbiont (S. microadriaticum) in early life stages of Acropora (fig. 1), isolated from an Acropora recruit collected at Ishigaki Island, Okinawa, Japan (Yamashita and Koike 2013). We compared its genome with the whole genome of A. tenuis to investigate possible genomic complementarity between the host coral and its symbiotic alga. We performed inoculation experiments with the strain using A. tenuis planula larvae without algal cells (fig. 1). Then we performed whole-genome transcriptomic analysis of the planula larvae to identify genes involved in symbiosis occurring during early life stages of Acropora in nature. Recent studies comparing transcriptomic responses of corals or anthozoans to different Symbiodiniaceae strains used various clades (different genera) (Yuyama et al. 2018; Lin et al. 2019). Here, we also compared transcriptomic responses of A. tenuis larvae with two closely related Symbiodinium strains (same genus), S. natans or type A2-relative (herein, type A2), which has never been detected in wild corals, and S. tridacnidorum, or type A3, which has occasionally been detected from natural corals (fig. 1; Yamashita et al. 2014, 2018). This was done in order to identify coral genes specifically involved in establishing and maintaining the symbiosis with native symbionts.

. 1

Experimental design and numbers of differentially expressed genes of Acropora tenuis inoculated with three Symbiodinium strains. (A) A colony of A. tenuis and Symbiodinium microadriaticum. (B) A schematic illustration of the experimental design. (C) Numbers of acquired Symbiodinium cells per planula larva (modified from Yamashita et al. [2018]; raw data are shown in supplementary table S3, Supplementary Material online). (D) Numbers of differentially expressed genes (DEGs) of A. tenuis. Circular areas (right) are proportional to numbers of DEGs. Upward arrows indicate upregulated genes, and downward arrows indicate downregulated genes compared with controls (without Symbiodinium inoculation). Numbers in the legend of the Venn diagram indicate each group (1, specific to A1; 2, A1–A3 common; 3, specific to A3; 4, A1–A2 common; 5, common; 6, A2–A3 common; and 7, specific to A2).

Results

Metabolic Complementation Analysis of the A. tenuis and S. microadriaticum Genomes

We obtained 692 Mb of S. microadriaticum genome assembly with an N50 size of 151 kb (supplementary table S1, Supplementary Material online). In total, 42,529 protein-coding genes (gene models) were identified in the genome assembly of S. microadriaticum (supplementary table S1, Supplementary Material online). Statistics of the genome assembly indicate that it is of comparable or better quality than other published genomes (supplementary table S1, Supplementary Material online). As reported in Aranda et al. (2016) and González-Pech et al. (2017), we confirmed that the GC content of the S. microadriaticum (AJIS2-C2) genome was higher than that of other Symbiodiniaceae genomes (Breviolum minutum and Fugacium kawagutii, supplementary table S1, Supplementary Material online), suggesting that this feature is probably common in S. microadriaticum, regardless of the strain (AJIS2-C2 isolated from Acropora sp. in Okinawa, Japan and CCMP2467 isolated from Stylophora in the Gulf of Aqaba; Aranda et al. 2016). In addition, the mean transcript length and numbers of exons per gene in S. microadriaticum (AJIS2-C2) are larger than those of the B. minutum genome (supplementary table S1, Supplementary Material online).

Metabolic dependencies between host and symbionts sometimes cause genome reduction and gene loss. For example, Buchnera, a bacterial endosymbiont in aphids, lacks some genes involved in essential amino acid pathways (Shigenobu et al. 2000; Hansen and Moran 2011), and Chlorella, a symbiotic alga in the cnidarian Hydra, lack some genes essential for nitrogen assimilation (Hamada et al. 2018). These studies indicate that complementarity and syntrophy between hosts and symbionts reflect their genome components. Complementation analysis of metabolic pathways between A. tenuis and S. microadriaticum genomes revealed that almost all enzymes involved in “Genetic information processing” were detected in both genomes (supplementary fig. S1, Supplementary Material online). Modules involved in “carbohydrate and lipid metabolism” (Lipopolysaccharide metabolism and Other terpenoid biosynthesis), in “nucleotide and amino acid metabolism” (Arginine and proline metabolism, Serine and threonine metabolism, Lysine metabolism, Histidine metabolism, and Aromatic amino acid metabolism), and in “energy metabolism” (Methane metabolism, Sulfur metabolism, Nitrogen metabolism, and Photosynthesis) were principally observed in the S. microadriaticum genome (>50% of all enzymes, supplementary fig. S1, Supplementary Material online). Most of these genes were not detected in the A. tenuis genome (supplementary fig. S1, Supplementary Material online). Comparative genomics revealed that Acropora corals may lack an enzyme essential for cysteine biosynthesis (Shinzato et al. 2011, 2020) and previous transcriptomic analyses suggested probable genetic complementarity between host corals and their symbionts (Shinzato et al. 2014; Lin et al. 2017). Taken together, A. tenuis may depend on its algal symbionts for all of these metabolites.

More Coral Genes Were Differentially Expressed When Inoculated with Native Algal Symbionts

In order to reveal the transcriptomic responses of Acropora corals to their native symbionts, we introduced three different Symbiodinium strains, native, occasionally symbiotic, and nonsymbiotic algae, to A. tenuis larvae (fig. 1). An average of 21.9 million RNA-seq reads per sample was retained after quality and adaptor trimming, 80% of which were mapped to A. tenuis gene models (supplementary table S2, Supplementary Material online). Hierarchical clustering analysis based on gene expression levels clustered our RNA-seq samples into three major clusters, according to days post fertilization (supplementary fig. S2, Supplementary Material online), indicating that overall gene expression patterns were more affected by time after fertilization or by planula developmental stage than by inoculation with Symbiodinium strains. Thus, in order to highlight genes involved in symbiosis and to ignore genes specific to larval development, we performed differentially expressed gene (DEG) identification for each of the sampling times.

As several genes were differentially expressed at two or three time points or days post fertilization, numbers of unique genes are also provided in parentheses (table 1). A total of 72 genes were assigned as DEGs across all three Symbiodinium types used to inoculate naïve coral larvae (fig. 1 and table 1, “common” group). Those DEGs included five essential, circadian rhythm-regulated genes (supplementary fig. S3, Supplementary Material online), suggesting that circadian rhythm-regulated genes are differentially expressed, regardless of Symbiodinium type.

Table 1

Exact Numbers of DEGs Shown in the Venn Diagram ( fig. 1D )

Number of DEGs	6 dpf	10 dpf	14 dpf	Total (Unique)
In A1-inoculated samples	216	175	135	526 (419)
Upregulated	85	73	58	216 (165)
Downregulated	131	102	77	310 (256)
Specific to A1	97	113	70	280 (258)
Upregulated	34	49	39	122 (108)
Downregulated	63	64	31	158 (151)
In A2-inoculated samples	28	34	41	103 (86)
Upregulated	18	19	16	53 (44)
Downregulated	10	15	25	50 (42)
Specific to A2	11	3	9	23 (23)
Upregulated	10	1	8	19 (19)
Downregulated	1	2	1	4 (3)
In A3-inoculated samples	122	130	69	321 (266)
Upregulated	52	70	22	144 (117)
Downregulated	70	60	47	177 (149)
Specific to A3	2	63	3	68 (68)
Upregulated	0	42	2	44 (44)
Downregulated	2	21	1	24 (24)
In common section	17	24	31	72 (55)
Upregulated	7	13	7	27 (18)
Downregulated	10	11	24	45 (37)

Note.—As several genes were overlapped at two or three time points or days post fertilization, numbers of unique genes are provided in parentheses.

We identified 526 DEGs from all three time points in larvae inoculated with type A1 Symbiodinium, 103 DEGs from all time points in those exposed to type A2, and 321 DEGs from all time points in those cultured with type A3 (fig. 1 and table 1). Although some authors have reported that smaller numbers of genes were differentially expressed when corals were cultured with competent algal symbionts (Voolstra et al. 2009; Schnitzler and Weis 2010). In our study, more genes were differentially expressed when A. tenuis larvae were cultured with native algal symbionts (type A1) (fig. 1). In 6-dpf samples, although the number of Symbiodinium cells in type A2-inoculated planula larvae (avg. five algae cells per planula) was larger than in planula larvae inoculated with type A3 (avg. one algae cell per planula) (supplementary table S3, Supplementary Material online) (Yamashita et al. 2018). The number of DEGs in type A3-inoculated samples (122 genes) was larger than that in type A2-inoculated samples (28 genes) (table 1), indicating that degrees of gene expression of A. tenuis were not simply determined by the number of algal cells per planula larva. Due to the small numbers (<0.01%) of sequencing reads mapped to the S. microadriaticum genome and gene models (mRNA), or previously reported S. microadriaticum gene models (supplementary table S4, Supplementary Material online), we were unable to perform transcriptomic analysis of Symbiodinium in planula larvae.

Co-expressed Gene Clustering Analysis of Each Symbiodinium-Inoculated Sample

To identify co-expressed genes that might be involved in common functions, we clustered genes based on their expression patterns under each type of Symbiodinium inoculation (tightness ten for A1-inoculated samples, tightness five for A2-inoculated samples, and tightness nine for A3-inoculated samples; supplementary table S5, Supplementary Material online). Clustering analyses identified 20 co-expressed gene clusters each in A1-, A2-, and A3-inoculated samples (supplementary figs. S4–S6, Supplementary Material online). In A1-inoculated samples, the number of A1-strain-specific DEGs was significantly enriched in co-expressed gene clusters, C6, C7, and C8 (χ2 test: P < 0.01, supplementary fig. S4, Supplementary Material online). In A3-inoculated samples, the number of A3-strain-specific DEGs was significantly enriched in cluster C1 (χ2 test: P < 0.01, supplementary fig. S5, Supplementary Material online), whereas no enriched clusters were observed in A2-inoculated samples (supplementary fig. S6, Supplementary Material online). Functional enrichment analysis using DAVID revealed that genes having UniProt keywords “Glycoprotein,” “Signal,” “Disulfide bond,” and “Coiled coil” were most commonly enriched (the number of genes containing these functions exceeded 10% of the total gene number in each cluster) in cluster C6, C7, and C8 in A1-inoculated samples (supplementary figs. S7–S9, Supplementary Material online). In A3-inoculated samples, UniProt keywords “Glycoprotein” and “Disulfide bond” were mainly enriched in cluster C1 (supplementary fig. S10, Supplementary Material online). Enriched UniProt keywords differed between A1- and A3-inoculated samples (supplementary figs. S7–S10, Supplementary Material online), possibly reflecting different symbiosis mechanisms of corals harboring A1 and A3 Symbiodinium.

Native Symbiodinium Strain-Specific Transcriptomic Responses of Host Corals

We compared the number of functionally annotated genes in each group of Venn diagram (fig. 2). In planula larvae inoculated with native type A1 Symbiodinium, 122 DEGs from all time points were specifically upregulated (fig. 1 and table 1). Among these 122 DEGs, 94 (77%) were annotated with Pfam clans (Pfam domains; fig. 2). Genes with structural functions (Tetratricopeptide repeat superfamily), catalytic functions (Membrane acyl transferase superfamily, Peptidase clan MA, and Six-hairpin glycosidase superfamily), immune functions (C1q and TNF superfamily), and signaling functions (Hedgehog/intein superfamily) were specifically detected in the A1-specific group (fig. 2), and genes having binding (Helix-turn-helix clan) or receptor functions (Family A G protein-coupled receptor-like superfamily) were continuously detected during the experiment (i.e., at 6-, 10-, and 14-dpf) in the A1-specific group (fig. 2).

. 2

Relative gene expression of functionary annotated differentially expressed genes. (A) Relative gene expression (row Z-score transformed) of upregulated differentially expressed genes with conserved domains during the experimental period (left) and numbers of genes with conserved domains significantly upregulated by three Symbiodinium strains (right). Domain categories in which at least two genes existed in each row were shown. About 6, 10, and 14 indicate the days post fertilization. Each group number in right heatmap below represents genes specific to each group on the Venn diagram in figure 1. (B) Relative gene expression (row Z-score transformed) of downregulated differentially expressed genes with KEGG modules during the experimental period (left) and numbers of genes with KEGG modules significantly downregulated by three Symbiodinium strains (right). Each group number in right heatmap below represents genes specific to each group on the Venn diagram in figure 1.

About 23 DEGs had transporter domains (fig. 2), and the number of genes having sugar transporter domains (Major facilitator superfamily) or amino acid transporter domains (amino acid-polyamine-organo-cation superfamily) was significantly larger than the number of DEGs specific to A2, A3, common, or others (Z-score >1.96, fig. 2, supplementary figs. S11 and S12, Supplementary Material online). These showed a tendency to increase as the period of exposure to Symbiodinium lengthened (fig. 2). Interestingly, some DEGs having sugar transporter and amino acid transporter domains belong to gene clusters in the A. tenuis genome (fig. 3 and supplementary fig. S13, Supplementary Material online). Molecular phylogenetic analysis revealed that these genes can be divided into two clades, Facilitated trehalose transporter 1 (Tret 1) and Solute Carrier Family 2 Member 8 (SLC2A8), in which a tandem duplication in Acropora was observed (fig. 3 and supplementary fig. S14, Supplementary Material online). Two DEGs (aten_s0094.g50 and aten_s0094.g54) were in each of the two clades containing sea anemone genes (fig. 3 and supplementary fig. S14, Supplementary Material online), indicating that both genes originated before the anthozoan diversification. Molecular phylogenetic analysis revealed genes in clade SLC2A8 were expanded in Acropora (or Complexa), and one of the genes (aten_s0094.g54) was involved in symbiosis with the native symbiont in A. tenuis (fig. 3). These two DEGs (aten_s0094.g50 and g54) were classified into C7 in a co-expressed gene clustering analysis based on their expression patterns in A1-inoculated samples (supplementary table S6, Supplementary Material online), suggesting that these serve a common function in symbiosis. Fossil and molecular data suggest most extant shallow-water, reef-building scleractinians fall into two major clades, Robusta and Complexa, to which Acropora belongs (Romano and Palumbi 1996). Molecular phylogenetic analysis of genes having transmembrane amino acid transporter domains revealed that these could be divided into three clades (fig. 3 and supplementary fig. S15, Supplementary Material online), and a DEG (aten_s0058.g54) and the two neighboring genes (g56 and g57) in the scaffold 0058 were clustered into a probable Acropora-specific (or Complexa-specific) transporter cluster (fig. 3 and supplementary fig. S15, Supplementary Material online), suggesting that Acropora acquired these genes after divergence of Complexa/Robusta. In contrast, molecular phylogenetic analysis of genes having an amino acid permease domain revealed these could be divided into distinct three clades, each of which contained sea anemones, Robusta corals, and Acropora (Complexa) corals (fig. 3 and supplementary fig. S16, Supplementary Material online). Interestingly, we found that the syntenic relationship of the A. tenuis gene cluster is partially conserved in sea anemones, and perfectly conserved within Robusta corals, including transcription directions (fig. 3 and supplementary fig. S16, Supplementary Material online), indicating that this gene cluster originated before the divergence of anemones and corals. Taken together, some A. tenuis genes duplicated before the anthozoan diversification, 550 Ma and after divergence of the Complexa/Robusta, 250 Ma (Park et al. 2012) are involved in symbiosis with native algae.

. 3

Molecular phylogenetic analyses of sugar and amino acid transporter clusters in the Acropora tenuis genome, including DEGs specific to native symbiont inoculation. (A) The gene-cluster (upper) and molecular phylogenetic analysis (lower) of sugar transporters having sugar (and other) transporter domains on scaffold 0094. Arrows indicate the transcription direction of genes. Possible gene names are also shown in boxes. Maximum likelihood phylogenetic relationships of genes having sugar and other transporter domains. About 438 gap-trimmed, aligned amino acids were used for phylogenetic tree reconstruction (CPREV model). Genes from the Acropora (Complexa), Robusta, sea anemone, and other bilaterians are colored red, blue, green, and purple, respectively. Circles on branches indicate bootstrap values higher than 80%. The bar indicates expected substitutions per site in aligned regions. The complete phylogenetic tree is also shown in supplementary figure S14, Supplementary Material online. (B) The gene cluster (upper) and molecular phylogenetic analysis (lower) of amino acid transporters having transmembrane amino acid transporter domains on scaffold 0058. Maximum likelihood phylogenetic relationship of genes having a transmembrane amino acid transporter domain. About 760 gap-trimmed, aligned amino acids were used for phylogenetic tree reconstruction (LG model). The complete phylogenetic tree is also shown in supplementary figure S15, Supplementary Material online. (C) The gene cluster (upper) and molecular phylogenetic analysis (lower) of amino acid transporters having amino acid permease domain on scaffold 0139. Maximum likelihood phylogenetic relationship of genes having amino acid permease domains. About 508 gap-trimmed, aligned amino acids were used for phylogenetic tree reconstruction (LG model). Information about the gene clusters of Nematostella vectensis, Exaiptasia pallida, and Orbicella faveolata were retrieved from the NCBI Genome Data Viewer. The complete phylogenetic tree is also shown in supplementary figure S16, Supplementary Material online.

We also compared numbers of DEGs downregulated by Symbiodinium inoculation, grouped into KEGG modules for each function. Downregulation of genes involved in amino acid metabolism (Aromatic amino acid metabolism and Cofactor and vitamin biosynthesis) was commonly detected in A1- and A3-inoculated samples (fig. 2). Eleven genes involved in carbohydrate metabolism (Other carbohydrate metabolism), lipid metabolism (Glycosaminoglycan metabolism), energy metabolism (Carbon fixation and Sulfur metabolism), nucleotide metabolism (Pyrimidine metabolism), and genetic information process (Proteasome, Ribosome, and Ubiquitin system) were specifically downregulated in the A1-inoculated sample (Z-score >1.96, fig. 2). Taken together, these gene classifications indicate that corals reduced amino acid, sugar, and lipid metabolism, and compensated for those metabolites by upregulating transporter genes after exposure to appropriate algal symbionts.

Genes for Which Expression Levels Were Dramatically Changed by Symbiodinium Inoculation

In order to identify genes for which expression was dramatically changed by Symbiodinium inoculation, we focused on those in which expression levels changed more than 2-fold (log2 ≥ |1|). In DEGs shared among the three samples cultured with different Symbiodinium types (72 genes; supplementary table S7, Supplementary Material online), expression levels of nine genes increased more than 2-fold (four genes were upregulated, and five genes were downregulated; supplementary table S7, Supplementary Material online). These included genes with circadian rhythm-related domains (e.g., Clock interacting protein circadian and PAS domain) and transcription factor domains (e.g., Basic region leucine zipper and Zinc finger domain) (supplementary table S7, Supplementary Material online). No genes in which expression levels changed more than 2-fold were found among DEGs specific to type A2- and A3-inoculated samples (supplementary table S7, Supplementary Material online).

We found that a DEG (aten_s0011.g21) showing similarity to a Ras-related and estrogen-regulated growth inhibitor (RERG), predicted to be involved in symbiosome formation (Mohamed et al. 2016), was significantly upregulated by inoculations involving all three Symbiodinium types (supplementary table S7 and fig. S17A, Supplementary Material online), and another RERG-like gene (aten_s0011.g24) was also significantly downregulated by all three Symbiodinium inoculations (supplementary fig. S17A, Supplementary Material online), suggesting that these two genes may possess different functions. Furthermore, we investigated expression patterns of genes putatively involved in symbiosis with a member of Symbiodiniaceae, B. minutum (Mohamed et al. 2016). Histone H2A (H2A-like; aten_s0063.g72), a gene reportedly involved in transcriptional control and cell cycle regulation, was significantly downregulated at three time points in A1-inoculated samples and was also significantly downregulated in A3-inoculated samples, 14-dpf (supplementary fig. S19B, Supplementary Material online). In addition, three genes involved in immunity, apoptosis, and stress responses (aten_s0002.g104, aten_s0002.g115, and aten_s0083.g44), Neuropeptide FF Receptor 2-like (NPFF2-like), were upregulated in A1-inoculated samples, 6-dpf (supplementary fig. S17C, Supplementary Material online). Moreover, expression levels of two genes similar to transient receptor potential cation channel subfamily A member 1 (TRPA1-like) were significantly changed in both type A1- and A3-inoculated samples at 6-, 10-, and 14-dpf, respectively (supplementary fig. S17C, Supplementary Material online).

Among 280 DEGs specific to type A1-inoculated samples at all time points (table 1), expression levels of 15 DEGs changed more than 2-fold (13 were upregulated and two were downregulated; supplementary table S7, Supplementary Material online). A sulfate ion transporter showed the highest expression change at log2 4.22 (14-dpf, supplementary table S7, Supplementary Material online). It has been also reported that another sulfate ion transporter gene was upregulated when A. tenuis juveniles were cultured with algal symbionts (Yuyama et al. 2016). Among 122 DEGs with expression levels upregulated by type A1 symbionts, lipid transporter (aten_s0225.g34, 14-dpf) showed the 6th largest change. Sugar transporters (aten_s0094.g54, samples at 14 dpf), which belong to the sugar transporter cluster (fig. 3), showed the 12th largest expression level changes (supplementary table S7, Supplementary Material online). These results reflect types of nutrients that corals obtain from native symbionts, using sulfate, sugar, and lipid transporters.

Genes with Expression Levels That Changed Continuously during the Experiments

Genes differentially expressed during the experimental period (planula larvae 6-, 10-, and 14-dpf) could be involved in general coral responses to Symbiodinium. Two genes similar to circadian locomotor output cycle protein kaput (CLOCK), a circadian rhythm-regulated gene, and a transcription factor, Thyrotroph embryonic factor (Tef), changed expression levels throughout the experimental period (supplementary table S8, Supplementary Material online), in a manner independent of Symbiodinium type. DEGs for which expression level constantly changed throughout experiments were not specifically detected in type A2- or A3-inoculated samples (supplementary table S8, Supplementary Material online).

In DEGs specific to type A1-inoculated samples, only one gene (aten_s0093.g61) similar to prosaposin (PSAP), a glycoprotein involved in antioxidant responses in human mitochondria (Kashiba et al. 2014), was constantly upregulated throughout the entire period (fig. 4 and supplementary table S8, Supplementary Material online). Interestingly, another PSAP-like gene is located in the same scaffold in the A. tenuis genome (scaffold 0093) (fig. 4). Molecular phylogenetic analysis showed that the two neighboring PSAP-like genes were clustered with other anthozoan genes (Bootstrap value = 94, supplementary fig. S18, Supplementary Material online) and that the node containing the DEG (aten_s0093.g61) includes exclusively Acropora genes (fig. 4), suggesting that these PSAP-like genes in the node are specifically expanded or are rapidly evolving in Acropora and may be involved in Acropora-specific responses to Symbiodinium. Actually, Acropora orthologs of PSAP-like genes including the A. tenuis gene (aten_s0093.g61) were under positive selection in Acropora (Shinzato et al. 2020).

. 4

Gene expression patterns and molecular phylogenetic analyses of differentially expressed prosaposin (PSAP)-like and Gretchenhangen-3 (GH3)-like genes. (A) Expression patterns (normalized counts) calculated by SLEUTH (upper) and molecular phylogenetic analysis (lower) of differentially expressed PSAP-like genes. Arrows on the box indicate directions of genes in the Acropora tenuis genome. Tandemly duplicated PASP-like genes in scaffold 0093 (middle). Maximum likelihood phylogenetic relationships of PSAP-like genes. About 385 gap-trimmed, aligned amino acids were used for phylogenetic tree reconstruction (FLU model). Genes from the Acropora (Complexa), Robusta, sea anemone, and other bilaterians are colored red, blue, green, and purple, respectively. Circles on branches indicate bootstrap values higher than 80%. The bar indicates expected substitutions per site in aligned regions. The sequence of surfactant-associated protein B (Sftpb) from Mus musculus was used as an outgroup. The complete phylogenetic tree is shown in supplementary figure S20, Supplementary Material online. (B) Expression patterns (upper) and molecular phylogenetic analysis (lower) of differentially expressed GH3-like gene. Maximum likelihood phylogenetic relationship of prosaposin-like genes. About 559 gap-trimmed, aligned amino acids were used for phylogenetic tree reconstruction (JTT model). Genes from the Acropora (Complexa), Robusta, other bilaterians, and land plants are colored red, blue, green, purple, and orange, respectively. The complete phylogenetic tree is shown in supplementary figure S21, Supplementary Material online.

Among DEGs shared by type A1- and A3-inoculated samples, one gene without homology in the Swiss-Prot database and one similar to an auxin-responsive gene, Gretchenhangen-3 (GH3) were constantly downregulated throughout the experiments (fig. 4 and supplementary table S8, Supplementary Material online). One of these (aten_s0217.g29) in A. tenuis was downregulated by coculturing with type A1 and A3 algae (fig. 4 and supplementary fig. S19, Supplementary Material online). Molecular phylogenetic analysis revealed that GH3-like genes are also observed in bilaterians, that coral GH3-like genes originated from an ancestral gene within the scleractinia and that most corals possess two GH3-like genes. Since other corals have genes phylogenetically close to the differentially expressed GH3-like gene in A. tenuis, functions of the GH3-like genes may be conserved among reef-building corals.

Genes Possibly Controlling Coral–Algal Symbiosis

In order to identify genes that may govern coral–algal symbiosis, we focused on signal molecules and transcription factors among DEGs. No downregulated DEGs with signaling domains were detected (supplementary fig. S20, Supplementary Material online). Three DEGs (both up- and downregulated) having transcription factor domains (Helix-loop-helix DNA-binding domain, Hairy Orange, and Basic region leucine zipper) were identified in the three Symbiodinium-inoculated samples (supplementary fig. S20, Supplementary Material online), indicating that these transcription factors responded to the presence of Symbiodinium, regardless of type. Although no transcription factor was differentially expressed commonly among type A1- and A2-cultured samples, ten genes with transcription factor domains were differentially expressed in both type A1- and A3-samples (supplementary fig. S20, Supplementary Material online).

Six DEGs having signaling and transcription factor domains (Hint domain, P53 DNA binding domain, Transcription factor DP, Homeobox domain, E2F/DP family winged-helix DNA-binding domain, and DM DNA binding domain) were exclusively detected in type A1-inoculated samples (supplementary fig. S20, Supplementary Material online), suggesting that these genes may control gene expression associated with native symbionts. In particular, two DEGs (aten_s0035.g9 and aten_s0035.g10) showing similarities to a ligand for the Notch signaling pathway, Delta-like protein 4 (Dll 4), were constantly upregulated in 10-, and 14-dpf samples (fig. 5). In addition, a neighboring gene (aten_s0035.g8) also showed similarities to another Notch ligand protein, jagged-1b (jag-1b) (fig. 5). Although expression patterns of these DEGs look similar (fig. 5), aten_s0035.g9 was classified into C7 and aten_s0035.g10 was classified into C8 in a co-expressed gene clustering analysis of A1-inoculated samples (supplementary table S6, Supplementary Material online), suggesting that these two genes may serve different functions in symbiosis with A1 symbionts. Interestingly, molecular phylogenetic analysis showed that these genes belong to neither canonical Dll nor Jagged (fig. 5 and supplementary fig. S21, Supplementary Material online). Although a Repeated Delta Serrate Ligand (RDSL) with unique domain composition has been reported from a sea anemone, Nematostella vectensis, (Gazave et al. 2009), these genes on scaffold 0035, including the two DEGs were not clustered with RDSL (fig. 5 and supplementary fig. S21, Supplementary Material online), and clustered into a novel clade; thus, we named it Anthozoan Notch ligand-like (AnNLL) (fig. 5 and supplementary fig. S21, Supplementary Material online). In addition, conserved protein domains of AnNLL were not detected, except for a Hint domain (fig. 5 and supplementary fig. S21, Supplementary Material online), suggesting that AnNLL genes may possess novel functions that have not been reported. In addition, clustering each of the three A. tenuis genes with Acropora digitifera and Acropora millepora counterparts (fig. 5 and supplementary fig. S21, Supplementary Material online) demonstrated that these genes were duplicated after divergence of the Complexa/Robusta and before the Acropora divergence. Novel Notch ligand-like genes, AnNLL, may have evolved to regulate gene expression required in the Acropora (or Complexa) clade for mutualism with algal symbionts.

. 5

Gene expression patterns, molecular phylogenetic analysis, and sequence feature of Notch ligand-like genes. (A) Gene expression patterns (normalized counts) by SLEUTH of Anthozoan Notch ligand-like genes located on scaffold 0035 of the Acropora tenuis genome assembly. Arrows on the box indicate transcription directions of the genes. (B) Molecular phylogenetic analysis and protein domain organizations of Anthozoan Notch ligand-like genes. The maximum likelihood phylogenetic tree was constructed using the WAG model based on 1,780 gap-trimmed, aligned amino acids. The branch length is the number of substitutions per site. Circles on branches indicate bootstrap values higher than 80%. Genes from the Acropora (Complexa), Robusta, sea anemone, and other bilaterians are colored red, blue, green, and purple, respectively. Sequence features are shown in the box as follows: MNNL, N terminus of Notch ligand; DSL, Delta serrate ligand; EGF, EGF-like domain; Cb-EGF, calcium-binding EGF domain; Hint, Hint domain. Note that these sequence features do not reflect the sequence length. The complete phylogenetic tree is shown in supplementary figure S23, Supplementary Material online.

Discussion

Symbiodinium Symbionts Affect Circadian Rhythms of Host Anthozoans

Circadian rhythm-regulating genes were differentially expressed regardless of the type of Symbiodinium employed in this study (supplementary fig. S3, Supplementary Material online). The intracellular pH of giant clams (Tridacna), which also support algal symbionts in extracellular space, is more alkaline in sunlight and acidic in the dark (Fitt et al. 1995), due to the effects of carbon dioxide, which decreases during photosynthesis and increases with nocturnal respiration (Miller and Yellowlees 1989). Likewise, intracellular pH of a coral (Stylophora) also tends to become more alkaline, depending on photosynthetic activities of intracellular symbionts (Laurent et al. 2013). Thus, animals harboring algal symbionts need to coordinate their circadian rhythms to adapt to intracellular environmental changes induced by Symbiodinium. Circadian rhythm-regulating genes exist in coral genomes (Levy et al. 2007), and these genes have been suggested to regulate circadian rhythms in corals (Levy et al. 2007). Moreover, it is reported that circadian rhythms of sea anemones (Aiptasia) harboring Symbiodiniaceae change from 12- to 24-h cycles as a result of inoculation with Symbiodiniaceae (Sorek et al. 2018). Thus, changing circadian rhythms of hosts could be a common feature in animals that harbor algal symbionts.

Is Symbiodinium Proliferation Regulated by Host Corals?

In this study, downregulation of the GH3-like gene was constantly observed in both type A1- and A3-inoculated samples (fig. 4). The GH3 gene is known as a phytohormone-regulating gene that regulates growth and development in land plants (Westfall et al. 2010). Although humans and mice possess a GH3-like gene containing an auxin-responsive promoter domain (fig. 4), no functional analysis has been reported yet. One of the most famous phytohormones, auxin, can be found in marine algae, including cyanobacteria, red-, green-, and brown-algae (Tarakhovskaya et al. 2007). The role of auxin is controversial, but it may be associated with growth in algae (Fürst-Jansen et al. 2020; Vosolsobě et al. 2020). In corals and Hydra, symbionts are potentially under control inside the animal host cell in which they are surrounded by a membrane, called a symbiosome or perialgal vacuole (Melo Clavijo et al. 2018). Although corals maintain an appropriate number of algal symbionts by expelling excess algal cells (Baghdasarian and Muscatine 2000; Yamashita et al. 2011), details of how algal symbionts in corals grow and proliferate, and how numbers of symbionts within corals are regulated remain enigmatic. It has been reported that plant and insect hosts control the growth of intracellular symbiotic bacteria by inhibiting bacterial cell differentiation (Van de Velde et al. 2010; Login et al. 2011). Expression changes of the GH3-like gene were observed when Symbiodinium types occurring in nature (A1: native symbiont, A3: occasional symbiont) are present. It has been reported that a nonsymbiotic sea anemone (Nematostella) does not possess a GH3-like gene in its genome (Ying et al. 2018) and we could not find it in the genome of a symbiotic anemone (Aiptasia, TBlastN, e-value cutoff: 10) (Baumgarten et al. 2015), suggesting that the GH3-like gene was lost in the anemone lineage. Therefore, it is predicted that the GH3-like gene may be involved in phytohormone regulation, controlling growth, and proliferation of Symbiodinium in corals. Functional analyses of coral GH3-like genes will be needed to confirm that these genes are actually involved in controlling growth and proliferation of algal symbionts.

Active Utilization of Native Symbionts by Host Corals

In photosymbiotic animals, the host shelters the symbiont against predators and environmental fluctuations, and supplies inorganic compounds, such as CO2, for photosynthesis. In turn, the host receives photosynthetic products, mainly in the form of fixed carbon, which can meet at least 50% of their overall nutritional requirements (Melo Clavijo et al. 2018). Metabolic exchanges between coral and algal symbionts are well known (Yellowlees et al. 2008; Davy et al. 2012). Upregulation of some transporter genes under illuminated conditions suggests nutrient uptake from algal symbionts in adult corals (Bertucci et al. 2015), although this was thought to be minimal in planula larvae (Kopp et al. 2016). In contrast to the previous report, downregulation of genes involved in sugar, amino acid, and lipid metabolism and upregulation of transporters for those compounds were observed exclusively in planula larvae inoculated with native algal symbionts in this study (fig. 2). Harii et al. (2010) suggested that planula larvae may use photosynthetic products from algal symbionts, given the high survival rate of planula larvae with algal symbionts and increased lipid levels under light conditions compared with dark conditions. In addition, loss of a cysteine biosynthesis enzyme from Acropora genomes suggests that cysteine is imported from algal symbionts in Acropora corals (Shinzato et al. 2011, 2020), and a possible cysteine transporter (SLC7A9) was specifically upregulated by inoculation with native algal symbionts (supplementary fig. S13, Supplementary Material online). Therefore, it is predicted that transporters identified in this study may be involved in nutrient uptake from algal symbionts in planula larvae. Upregulation of transporters and downregulation of genes involved in sugar, amino acid, and lipid metabolism can occur only when corals acquire appropriate algal symbionts.

Modulation of Host Transcriptomes by Inoculation with Different Symbiodinium Strains

In this study, A. tenuis showed clearly different transcriptome responses for each genetically closely related Symbiodinium strain (fig. 2). In the repertoires of DEGs in A1-inoculated samples (A1; S. microadriaticum, native symbiont in early life stages of Acropora), genes reportedly involved in symbiosis with a Symbiodiniaceae species, B. minutum (Mohamed et al. 2016) (RERG-, NPFF2-, and TRPA1-like genes) were present (supplementary fig. S17, Supplementary Material online). These results suggest that symbiosis mechanisms, such as arrested phagosomes and modulation of immune responses, proposed by Mohamed et al. (2016) may occur in the A. tenuis planula larvae inoculated with native symbionts. Whereas, DEGs such as RERG-, NPFF2-, and TRPA1-like genes were also detected in A2- and/or A3-inoculated samples (supplementary fig. S17, Supplementary Material online), indicating that arrested phagosomes and modulated immune responses are not specific to symbiosis with the native symbiont in A. tenuis. In addition, fewer DEGs were involved in those mechanisms than were reported by Mohamed et al. (2016) (five genes in this study and 52 genes in Mohamed et al. [2016]), indicating that other mechanisms, for example, upregulation of transporter genes and metabolic suppression, are important for symbiosis with native symbionts.

In contrast to A1-inoculated samples, relatively few DEGs were detected in A2-inoculated samples (A2; S. natans, has never been detected in natural corals) (fig. 1), and DEGs specific to A2 were also limited (fig. 1 and table 1). Interestingly, two genes in the galactose-binding domain-like superfamily (aten_s0030.g1 and aten_s0138.g31) were specifically upregulated in A2-inoculated samples (fig. 2 and supplementary fig. S11, Supplementary Material online). Many members of this superfamily have motifs that are involved in carbohydrate recognition on cell surfaces (Murzin and Bateman 1997). Algal cell surfaces possess glycan ligands, such as α-mannose/α-glucose and α-galactose, which are recognized by corals as microbe-associated molecular patterns during initial contact (Koike et al. 2004; Wood-Charlson et al. 2006; Jimbo et al. 2013). Lectins specific to each carbohydrate were isolated from corals (Jimbo et al. 2000, 2005, 2010; Kvennefors et al. 2008, 2010; Vidal-Dupiol et al. 2009; Kuniya et al. 2015). In addition, symbiotic phytoplankton, including Symbiodiniaceae, possesses a range of oligosaccharide molecules on their surfaces that vary within and among species (Markell et al. 1992; Tien et al. 2005). In land plants, lectin motifs enable plants to recognize invading organisms, such as pathogens, and to resist pathogen colonization (Lannoo and Van Damme 2014). Although Kuniya et al. (2015) reported that Tachylectin-2-like lectin is involved in acquisition of another type A3 algal symbiont (NBRC102920), during the juvenile polyp stage of A. tenuis, no Tachylectin-2-like genes were differentially expressed in this study. These observations suggest that two genes in the Galactose-binding domain-like superfamily (aten_s0030.g1 and aten_s0138.g31) specifically upregulated in A2-inoculated samples may enable A. tenuis to recognize inappropriate algae, like type A2, rather than facilitating acquisition of algal symbionts.

In the occasional symbiont A3-inoculated samples (A3; S. tridacnidorum, is occasionally detected in natural Acropora corals). Almost all DEGs, including genes in the common group, were shared with A1-inoculated samples (119 genes [98%] at 6-dpf, 61 genes [47%] at 10-dpf, and 65 genes [94%] at 14-dpf; fig. 1), suggesting that they share some molecular mechanisms for symbiosis. For example, upregulation of sugar and amino acid transporter genes and downregulation of genes involved in amino acid metabolism in the A1–A3 common group (fig. 2) indicate that nutrient uptake by corals from algal symbionts normally occurs when A. tenuis harbors these symbionts.

Gene Duplications May Have Enabled Co-option of Symbiosis Genes in Corals

We showed that some genes in tandemly duplicated gene clusters (e.g., Notch ligand-like, sugar and amino acid transporters, PSAP-like genes) are upregulated when native algal symbionts are present (figs. 3–5). The Notch signaling pathway is a conserved signaling pathway in metazoans (Artavanis-Tsakonas et al. 1999) and it regulates cell fates during early development in humans (Kopan and Ilagan 2009; Marlow et al. 2012). It was suggested that the Notch signaling pathway interacts with the NF-kB signaling pathway (Ang and Tergaonkar 2007), which is involved in regulation of innate and adaptive immunity in various animals, ranging from insects to humans (Gilmore and Wolenski 2012). It has been suggested that corals recognize their algal symbionts using their immune systems (Wood-Charlson et al. 2006). Therefore, immune systems originating from the Notch signaling pathway may be involved in specific recognition of native symbionts in corals.

Stable and persistent photosymbioses in the animal kingdom are thought to be dependent on cellular response mechanisms coping with elevated levels of reactive oxygen species (ROS) (Melo Clavijo et al. 2018). Although the functions of PSAP genes in cnidarians are unclear, knockdown of a PSAP gene in humans caused a decrease in coenzyme Q10, an antioxidant, suggesting that the PSAP gene may regulate coenzyme Q10 (Kashiba et al. 2014). The presence of coenzyme Q10 in corals has been documented and it functions as an antioxidant during oxidative stress (Lutz et al. 2014). ROS, which cause oxidative stress and DNA damage, are constantly produced in coral mitochondria and/or plastids of algal symbionts (Weis 2008; Buxton et al. 2012). Thus, preventing oxidative stress is important for corals harboring algal symbionts, and the PSAP-like gene may serve this function.

It has been reported that many amino acid transporter genes were tandemly duplicated in the genomes of aphids that harbor endo-symbiotic bacteria (Price et al. 2011), and that these genes are especially expressed in bacteriocytes or guts, suggesting that they participate in symbiotic mechanisms (Price et al. 2011; Dahan et al. 2015). In addition, tandemly duplicated amino acid transporter genes occur in genomes of other sap-feeding insects, citrus mealybugs, cicadas, and potato psyllids (Duncan et al. 2014). Gene duplication is a common phenomenon in many organisms and contributes to acquisition of new gene functions during evolution (Ohno 1970). As in symbioses between insects and bacteria, corals may acquire genes involved in symbiosis with algal symbionts via gene duplication. In addition, Shoji and Hashimoto (2011) reported that a transcription factor recognizes promotor regions of several duplicated nicotine biosynthesis genes in the tobacco genome and controls their expression levels; thus, gene duplication may enhance specificity and responsiveness. Such effective transcription systems may exist in Acropora and may be used for efficient control of gene transcription during establishment of mutualism with native algal symbionts.

Conclusions

In this study, we discovered molecular mechanisms of symbioses during early life stages of corals and Symbiodinium in nature by comparing coral transcriptomic responses to native symbionts and to algae that are phylogenetically close, but that involve different symbiotic processes. When inoculated with native algal symbionts, corals slow metabolic rates and actively acquire photosynthetic products from algal symbionts. Corals may have acquired genes involved in these processes by gene duplication and we were able to clearly identify some of these among tandem duplicated genes. Further investigations of genes, such as duplicated sugar and amino acid transporters, novel Notch ligand-like genes that may control symbiosis, and GH3-like genes, which may control growth and proliferation of algal symbionts, should further illuminate molecular mechanisms of coral–algal symbiosis.

Materials and Methods

Sampling, Genomic DNA Isolation, and Whole-Genome Assembly of S. microadriaticum

The culture strain AJIS2-C2 (type A1; S. microadriaticum) was originally isolated from a coral (Acropora sp.) at Ishigaki Island, Okinawa, Japan (Yamashita and Koike 2013). Genomic DNA isolation, genome assembly, and gene prediction were performed as previously described (Shoguchi et al. 2013, 2018). Briefly, genomic DNA was isolated using the phenol–chloroform method. DNA was fragmented into ∼600-bp lengths and 200 ng of DNA were used for PCR-free shotgun library preparation. For mate-pair libraries, DNA of different sizes (∼3, 7, 10, and 15 kb) was separated using SageEFL, whereas a Nextera Mate Pair Library Prep Kit (Illumina) was used for library preparation following manufacturer instructions. The library was sequenced using a HiSeq 2000 (Illumina). Sequencing data were assembled de novo with the Platanus genome assembler v1.2.440. Subsequent scaffolding was performed with SSPACE v3.0 (Boetzer et al. 2011) using Illumina mate-pair information, and gaps inside scaffolds were closed with paired-end data using Gapcloser (Luo et al. 2012). To identify erroneous scaffolds created by genome assembly errors or contamination, filtered raw reads were mapped against the preassembly to examine scaffold coverage. Then, scaffolds with exceptionally high or low coverage (Kolmogorov–Smirnov test, P < 0.05) were removed from subsequent analyses. Genome assembly statistics were compared with those of other Symbiodiniaceae genomes, S. microadriaticum (Aranda et al. 2016), B. minutum (Shoguchi et al. 2013), and F. kawagutii (Lin et al. 2015). As intron–exon boundaries are uncommon in genomes of dinoflagellates, and three major dinucleotides (GT/GC/GA) of introns were found in Symbiodiniaceae genomes (Shoguchi et al. 2013, 2018; Aranda et al. 2016), in order to allow intron boundaries, a set of gene models was predicted using AUGUSTUS v3.2.3 using the “–allow_hinted_splicesites = GAAG” option (GTAG and GCAG splice sites are already allowed in the default setting of AUGUSTUS) following Shoguchi et al. (2013) and Shoguchi et al. (2018).

Gene Annotations of A. tenuis and S. microadriaticum Gene Models

Gene models of A. tenuis were obtained from Shinzato et al. (2020). Both A. tenuis and S. microadriaticum gene models were annotated as follows: 1) BlastX v2.2.28+ (Camacho et al. 2009) homology search with an e-value cut off of 1e−5 against the Swiss-Prot database (April 9, 2018), 2) searching for conserved protein domains using the Pfam database (Finn et al. 2016) with InterProScan v5.27-66.0 (Jones et al. 2014) with default settings, and 3) searching for metabolic pathways using KAAS (Kyoto Encyclopedia of Genes and Genomics Automatic Annotation Service) (Kanehisa and Goto 2000; Moriya et al. 2007), with the one-directional best hit method. For each KEGG metabolic pathway, percentages of enzymes in the A. tenuis or S. microadriaticum genomes were calculated with a custom Perl script, as reported in (Guzman et al. 2018). Annotated protein domains were further classified into clans, using the hierarchy of related Pfam entries (Pfam 30.0, posted July 1, 2016).

Inoculation Experiments

The culture strain ISS-C2-Sy (type A2-relative; S. natans) was originally isolated from coral reef sand at Ishigaki Island, but has never been detected in wild corals (Yamashita and Koike 2013). The culture strain CS-161 (type A3; S. tridacnidorum), which has occasionally been detected in wild corals, was purchased from the Australian National Algae Culture Collection, Australia. In this study, we designated S. microadriaticum as type A1, S. natans as type A2, and S. tridacnidorum as type A3. Inoculation experiments using these three strains were conducted as in Yamashita et al. (2018), and planula larvae collected by Yamashita et al. (2018) were used in this study. Briefly, we added each Symbiodinium strain to planula larvae every day from 2- to 5-dpf (more details are available in Yamashita et al. [2018]), and planula larvae were kept in an incubator (27 °C; 80–120 µmol photons m−2 s−1; 12 h light:12 h dark cycle). Then, 6-, 10-, and 14-dpf planula larvae that were cultured with each Symbiodinium strain, and planula larvae without Symbiodinium exposure were collected (12 conditions corresponding to the four conditions [A1, A2, A3, and apo-sym] time the three sampling times). Numbers of acquired Symbiodinium cells per planula larva are shown in figure 1, supplementary table S3, Supplementary Material online, and Yamashita et al. (2018). Collected planula larvae were snap frozen using acetone/dry ice and stored at −80 °C until use.

RNA Sequencing and Transcriptomic Analyses for Identification of DEGs

Total RNA was extracted using an RNeasy Plant Mini Kit (QIAGEN). A TruSeq RNA Sample Preparation Kit (Illumina) was used for sequencing library preparation, and each library was sequenced from 100-bp paired-end libraries using a HiSeq 2000 (Illumina). Low-quality reads (quality score <20 and length <20 bp) and Illumina sequence adaptors were trimmed with CUTADAPT v1.16 (Martin 2011), and then mapped to the A. tenuis gene models (mRNA) using KALLISTO v0.44.0 (Bray et al. 2016) with 100 bootstrap replicates. Normalized gene expression levels (normalized counts) of each Symbiodinium-inoculated sample (supplementary data 1, Supplementary Material online), which were calculated with SLEUTH v0.30.0 (Pimentel et al. 2017) in R v3.4.1 (https://www.r-project.org/, accessed December 4, 2020), were compared with control samples (without Symbiodinium), and samples inoculated with other Symbiodinium types at each time point (6-, 10-, and 14-days post fertilization). Then statistical tests to identify DEGs were performed using SLEUTH in R as well. P values were adjusted using the Benjamini–Hochberg method in SLEUTH. Genes that satisfied the following two criteria were designated as DEGs: 1) gene expression level was significantly different (q < 0.05) from control samples and 2) gene expression level was significantly different (q < 0.05) compared with samples cultured with other Symbiodinium strains in at least one comparison. DEGs were classified into seven groups in a Venn diagram (fig. 1) and were named “specific to A1,” “A1–A3 common,” “specific to A3,” “A1–A2 common,” “common,” “A2–A3 common,” and “specific to A2.”

For detecting Symbiodinium transcripts in planula larvae, the above RNA-seq reads were mapped to the S. microadriaticum genome using HISAT v2.1.0 (Kim et al. 2015) without penalties for pairs of noncanonical splice sites (–pen-noncansplice 0) due to a large number of noncanonical intron boundaries in Symbiodiniaceae genomes (Shoguchi et al. 2013, 2018). We also mapped reads to gene models in this study and gene models from a previously reported S. microadriaticum (CCMP2467) genome (Aranda et al. 2016) using KALLISTO.

Hierarchical Clustering Analysis Based on Gene Expression Levels

Pearson’s correlation coefficients between all samples were calculated based on their gene expression levels calculated with SLEUTH after removing low-expression genes (fewer than ten normalized counts) using the “genefilter” package in R software. Then these were hierarchically clustered using “hclust” with the “ward.D2” method in R software.

Clustering Co-expressed Genes Based on Gene Expression Patterns of Each Symbiodinium-Inoculated Sample

In order to identify co-expressed gene clusters, genes were clustered using Clust (Abu-Jamous and Kelly 2018), based on expression patterns induced by each Symbiodinium type inoculation with different tightness (1–10). Tightness means the weight of cluster tightness versus cluster size. Tightness values smaller than 1 produce larger and more diffuse clusters, and tightness values larger than 1 produce smaller and tighter clusters. The tightness that produced the largest number of clustered genes was selected. We named each cluster C1 to C20, as shown in supplementary figures S4–S6, Supplementary Material online. Then, numbers of strains A1-, A2-, and A3-specific DEGs in each cluster were examined. A χ2 test was performed to identify differences in the percentage of DEGs in each cluster (“the number of DEGs”/“the number of total predicted genes” and “the number of DEGs in each cluster“/”the number of assigned genes in each cluster”). Clusters in which DEGs were significantly enriched were functionally annotated using UniProt Keyword via web platform DAVID v6.8 (Huang et al. 2009) with a background that comprised all of the UniProt IDs from gene models of A. tenuis.

Molecular Phylogenetic Analyses

Amino acid sequences of homologous genes from cnidarian genomes were retrieved using the BlastP server (https://blast.ncbi.nlm.nih.gov/Blast.cgi, accessed December 4, 2020) with default settings against the NCBI RefSeq database on February 27, 2020. Homologous amino acid sequences of bilaterians and land plants were retrieved using BlastP homology searches against the Swiss-Prot database with an e-value cutoff 1e−5. Amino acid sequences were aligned using MAFFT v7.310 (Katoh 2002; Katoh and Standley 2013) and gaps in alignments were removed using trimAl v1.2 (Capella-Gutiérrez et al. 2009) with the “gappyout” option. After removing gaps, maximum likelihood analyses were performed using RAxML v8.2.10 (Stamatakis 2014) with “bootstrap 100” and “protgammaauto” options.

Supplementary Material

Supplementary data are available at Genome Biology and Evolution online.

Click here for additional data file.