Comparison of Symbiodiniaceae diversities in different members of a Palythoa species complex (Cnidaria: Anthozoa: Zoantharia)—implications for ecological adaptations to different microhabitats

Specimens collection

Eighty-two colonies of three Palythoa species (P. tuberculosa, P. sp. yoron, and P. mutuki) were collected from a shallow fringing reef of Tokunoshima Island, Kagoshima, Japan (Figs. 1 and 2). Specimens of these three Palythoa species were collected in four different areas (Table 1, Fig. 2A): reef edge (Fig. 2B, 27.76998333N, 129.03988611E) for P. tuberculosa (Fig. 2C); reef flat 1 (Fig. 2D, 27.76997777N, 129.03925000E) for P. tuberculosa (Fig. 2E) and P. mutuki; reef flat 2 (Fig. 2F, 27.77195277N, 129.03843611E) for P. mutuki (Fig. 2G); and backreef moat (Fig. 2H, 27.76990833N, 129.03855833E) for P. tuberculosa and P. sp. yoron (Fig. 2I). To avoid collecting clones, we collected individuals from clearly different colonies while maintaining a set distance from each other of at least 1 m. In a previous study, even when closer to each other (within approximately 50 × 50 cm), no clones were observed in Zoanthus (Cnidaria: Anthozoa: Zoantharia) colonies (Albinsky et al., 2018). In addition, eighteen previously collected specimens of Palythoa species including 10 P. aff. mutuki specimens from Mizuyama, Masucci & Reimer (2018) were also examined in this study (Table 1).

DNA extraction and PCR amplification

From each of these specimens, several polyps were cut with a surgical knife and DNA was extracted using DNeasy Blood and Tissue Kit (QIAGEN). DNA concentrations were checked by Qubit Fluorometer (ThermoFisher, Waltham, USA). Two molecular markers for genotyping symbiotic algae of Palythoa species were examined: nuclear internal transcribed spacer ribosomal DNA (ITS-rDNA) region including partial 18S–ITS1–5.8S–ITS2–partial 28S (primers zITSf: CCG GTG AAT TAT TCG GAC TGA CGC AGT and ITS4: TCC TCC GCT TAT TGA TAT GC, (Baillie, Belda-Baillie & Maruyama, 2000; appx. 700–750 bp) and plastid mini-circle non-coding region DNA (psbA^ncr) (primers 7.4-Forw: GCA TGA AAG AAA TGC ACA CAA CTT CCC and 7.8-Rev: GGT TCT CTT ATT CCA TCA ATA TCT ACT G, (Moore et al., 2003; appx. 800–850 bp). These regions were amplified according to the PCR thermal conditions in Wee, Kurihara & Reimer (2019). Amplified PCR products of symbionts were directly sequenced, and sequence data were manually checked based on the chromatogram files and low quality sites were removed at the 5′ and 3′ ends by BioEdit v.7.0.5.3 (Hall, 1999). Obtained sequences were deposited in the GenBank database (MN654128–MN654306, Table 1).

Haplotype network inference and phylogenetic estimation

Obtained sequences for ITS-rDNA, psbA^ncr forward and reverse regions were aligned, respectively. In order to discriminate taxa of Symbiodiniaceae, we extracted the ITS2 region utilizing SymPortal (Hume et al., 2019; https://symportal.org/) and performed BLASTN search against the nt database using the NCBI website (https://blast.ncbi.nlm.nih.gov/Blast.cgi) for ITS-rDNA sequences. Haplotype network inference was performed for ITS-rDNA sequences using the alignment with TCS networks method (Clement et al., 2002) in PopART (Leigh & Bryant, 2015). Any columns in the alignment with gaps or ambiguous sites were automatically masked in the inference. The phylogenetic analyses were performed by MEGA version X (Kumar et al., 2018) and any loci with ambiguous (double peaks) sites and gaps was automatically deleted completely for calculation in order to avoid over/underestimation of genetic distance among each sequence. Molecular phylogenetic trees of each marker were constructed by maximum likelihood (ML) and neighbor joining (NJ) methods under the JC+G model for ITS-rDNA region and the JC model for psbA^ncr regions adopted by modeltest program within MEGA X. The significance of each node was tested by bootstrap test with 1,000 replications. Bayesian inference was performed using BEAST2 (Bouckaert et al., 2019) under default settings other than the clock model being changed to the relaxed log normal model, which showed the highest likelihood value according to the model comparison program compiled in BEAST2 (Drummond et al., 2006). Posterior probability (PP) on each branch was calculated summarizing four independent 10 million MCMC simulations.

Statistical analyses

To clarify the relationships between (1) symbiont lineages and host species, and (2) symbiont lineages and host microhabitats, Fisher’s exact test was conducted for the compositions of genotype for ITS-rDNA region and monophyletic clades for psbA^ncr forward and reverse regions. It should be noted that host microhabitat was restricted by host species for P. sp. yoron and P. mutuki, and thus we only targeted P. tuberculosa for these analyses (aim 2 above) When significance was detected in Fisher’s exact test, Cramér’s coefficient of association (V) was calculated to evaluate which factors (host species or host microhabitat) were strongly associated with each other.

Sequence alignment

The total number of sequences of Symbiodiniaceae from specimens of the four Palythoa species obtained in this study was 98 sequences for the ITS-rDNA region (513–773 bp), 40 sequences for the psbA^ncr forward region (330–547 bp), and 41 sequences for the psbA^ncr reverse region (352–494 bp). As the primer set for psbA^ncr used in this study did not make a congruent contig, obtained sequences of forward regions and reverse regions were aligned separately (Noda et al., 2017). After alignment, a total of 449 sites with 5 parsimony informative (=PI) sites for the ITS-rDNA region, 260 sites with 94 PI sites for the psbA^ncr forward region, and 293 sites with 40 PI sites for the psbA^ncr reverse region were used for each phylogenetic estimation.

Barcoding, haplotype network and phylogenetic trees

As the result of BLAST searches, all query sequences of the ITS-rDNA region (n = 98) were confirmed as belonging to the genus Cladocopium. Seventeen ITS-rDNA unique sequences (=genotypes) were observed in TCS network, with most of the sequences belonging to one of major three ITS-rDNA genotypes (Fig. 3, Table S1). No significant clade was detected for the ITS-rDNA phylogenetic tree (Fig. S1). Summarizing these ITS-rDNA genotypes from the viewpoint of host species, P. tuberculosa possessed mainly Genotype01 (n = 20) followed by Genotype02 (n = 7), and P. sp. yoron also possessed mainly Genotype01 (n = 20) followed by Genotype03 (n = 8) (see details in Table S1). On the other hand, P. mutuki possessed mainly Genotype02 (n = 13) with a few Genotype01 (n = 3) and Genotype03 (n = 2). Although the number of specimens examined was smaller (n = 6) than those the other species, P. aff. mutuki also possessed mainly Genotype01 (n = 5).

In contrast, phylogenetic trees generated from psbA^ncr regions had a higher resolution. Two monophyletic clades were well supported by bootstrap values and posterior probability in both forward (Fig. 4 clf1, ML = 100, NJ = 100, PP = 1 and clf2, ML = 100, NJ = 100, PP = 1) and reverse trees (Fig. 5 clr1 and clr2, ML = 100, NJ = 100, PP = 1). Summarizing these Symbiodiniaceae lineages from the viewpoint of host species, P. tuberculosa inhabiting the reef edge possessed clf1/clr1 lineage (n = 7∕5) and one specimen inhabiting at the backreef moat possessed clf2/clr2 lineage. Palythoa sp. yoron inhabiting at the backreef moat possessed mainly clf1/clr1 (n = 9∕13), however, approximately one third of specimens (n = 5) possessed other lineages. On the other hand, P. mutuki inhabiting the reef flat possessed mainly clf2/clr2 (n = 8∕8) other than two specimens that possessed clf1/clr1. Unfortunately, as most of P. aff. mutuki specimens were not amplified by this primer set, we could only obtain phylogenetic information on one specimen which possessed the same lineage as P. sp. yoron (C24ToKa-PF) for the forward region and clr1 for the reverse region.

Relationships among symbiont genotype/lineages, host species and host microhabitats

From the results of Fisher’s Exact test, significant differences were detected in all combinations, i.e., ITS-rDNA genotype and host species (p < 0.01), psbA^ncr forward lineages and host species (p < 0.01), psbA^ncr reverse lineages and host species (p < 0.01), and ITS-rDNA genotype and host microhabitats for P. tuberculosa (p < 0.05) (Table 2). In other words, it was shown that Symbiodiniaceae lineages and host species were not independent, nor were Symbiodiniaceae lineages and host microhabitats for P. tuberculosa. The effective dose calculated by Cramér’s coefficient of association (V) was largest between host species and psbA^ncr forward/reverse lineages (V = 0.786, V = 0.682, respectively), and moderate for the other combinations (host species and ITS-rDNA genotypes, V = 0.477; host microhabitats and ITS-rDNA genotypes).

Symbiodiniaceae genotype/lineage and host species

The development of molecular markers such as psbA^ncr that have higher resolution than commonly used 18S or ITS ribosomal DNA markers has helped unveil a more detailed picture of the genetic diversity of Symbiodiniaceae (Takishita et al., 2003; LaJeunesse & Thornhill, 2011; LaJeunesse et al., 2018) (but see also Hume et al., 2019 who utilized intragenomic variation of ITS2 to resolve genetic delineations). Accordingly, host species biodiversity has been discovered from the initial observation of differences of Symbiodiniaceae phylotypes in some cnidarian species (e.g., gorgonian Eunicea flexuosa, Prada et al., 2014; scleractinian coral Seriatopora hystrix, Warner, Van Oppen & Willis, 2015).

From the results of Mizuyama, Masucci & Reimer (2018), none of the four molecular markers utilized could clearly delineate four Palythoa species, although they could delineate two closely related species groups composed of P. tuberculosa—P. sp. yoron and P. mutuki—P. aff. mutuki. These previous results seem to be reflected in the results in the current study of Symbiodiniaceae genotypes of ITS-rDNA and lineages of psbA^ncr regions. Palythoa tuberculosa and P. sp. yoron mostly shared the same symbiont genotype (Genotype01); nevertheless, they also partially shared the other genotypes with P. mutuki (Genotype02 and Genotype03). With regard to psbA^ncr lineages, even though the delineation of species groups between P. tuberculosa—P. sp. yoron and P. mutuki were shown more clearly, they were not divided completely. The situation requires further investigation via obtaining more P. aff. mutuki specimens’ psbA^ncr sequences. Unfortunately, in the current study, despite much searching, we could not find large numbers of P. aff. mutuki on the reef in Tokunoshima Island, even though they were previous sampled for Mizuyama, Masucci & Reimer (2018). We do not know what happened to P. aff. mutuki colonies, but they may have been strongly affected by the bleaching events of 2016 and 2017 observed in southern Japan (Masucci et al., 2019).

Symbiodiniaceae genotype/lineage and microhabitat of host species

From the results of the phylogenetic analyses, three microhabitats were not exclusively allocated in distinct Symbiodiniaceae genotypes or monophyletic clades, but the ratios of different genotypes were significantly different for P. tuberculosa. Regarding P. tuberculosa, Symbiodiniaceae Genotype01 was mostly detected on the reef edge and reef flat, while Genotype02 was mainly observed in the backreef moat. Although there were not enough samples to conduct statistical examinations of P. sp. yoron and P. mutuki due to their habitat specificity, Genotype02 and clf2/clr2 were detected mainly on the reef flat, while Genotype01 and clf1/clr1 were observed from all three environments.

It has been reported that zoantharian species with different symbiotic genotypes show species-specific photosynthetic responses against seawater temperature and p CO₂ (Graham & Sanders, 2016; Reimer et al., 2017b; Wee, Kurihara & Reimer, 2019). Although the four Palythoa species in this study occurred sympatrically on one reef, the environmental conditions in a reef can be quite different according to small-scale geographical features. Seawater temperatures on reef flats frequently reach near 40 °C (Achituv & Dubinsky, 1990). In enclosed reefs, seawater temperatures and p CO₂ show higher variations than those in exposed reefs (Suzuki, Nakamori & Kayanne, 1995; Fitt et al., 2001). Thus, the relationship between Symbiodiniaceae and host Palythoa species may change among different microhabitats in a reef area, facilitating ecological divergence of Palythoa species within a narrow geographic range.

Although a previous molecular study could not distinguish the boundaries among these Palythoa species (Mizuyama, Masucci & Reimer, 2018), it is suggested by our results that these species are ecologically divergent, and physiological differences within Symbiodiniaceae species may contribute to their ecological adaptation. In fact, Howells et al. (2012) reported that Cladocopium C1 in Acropora tenuis showed different physiological responses between northern and southern populations in the Great Barrier Reef. Considering that Cladocopium contains various species distinguished by differences of only a few bp in the ITS2 maker (Thornhill et al., 2014), meta-barcoding analyses via next-generation sequencing would be necessary to further understand the detailed relationship between Symbiodiniaceae and Palythoa species complex.