A preliminary survey of zoantharian endosymbionts shows high genetic variation over small geographic scales on Okinawa-jima Island, Japan

Symbiotic dinoflagellates (genus Symbiodinium) shape the responses of their host reef organisms to environmental variability and climate change. To date, the biogeography of Symbiodinium has been investigated primarily through phylogenetic analyses of the ribosomal internal transcribed spacer 2 region. Although the marker can approximate species-level diversity, recent work has demonstrated that faster-evolving genes can resolve otherwise hidden species and population lineages, and that this diversity is often distributed over much finer geographical and environmental scales than previously recognized. Here, we use the noncoding region of the chloroplast psbA gene (psbAncr) to examine genetic diversity among clade C Symbiodinium associating with the common reef zoantharian Palythoa tuberculosa on Okinawa-jima Island, Japan. We identify four closely related Symbiodinium psbAncr lineages including one common generalist and two potential specialists that appear to be associated with particular microhabitats. The sea surface temperature differences that distinguish these habitats are smaller than those usually investigated, suggesting that future biogeographic surveys of Symbiodinium should incorporate fine scale environmental information as well as fine scale molecular data to accurately determine species diversity and their distributions.


INTRODUCTION
Symbiodinium is an important genus of dinoflagellate photosymbionts found in tropical and subtropical marine ecosystems. These "zooxanthellae" transfer energy to their invertebrate hosts in nutrient-poor environments, enhancing the growth of reef-building corals and other reef organisms such as zoantharians (Muscatine & Cernichiari, 1969;Baker, 2003). They also serve a key role in establishing the thermal tolerance of coral colonies and shaping the adaptive response of reef organisms to climate change (Sampayo et al., 2008;Thornhill et al., 2014). As climate change intensifies and diversity patterns are expected to alter, it is increasingly important to map marine species distributions on scales both large (McClanahan et al., 2014;Reimer et al., 2017) and small (Mieszkowska & Lundquist, 2011).
The biogeography of Symbiodinium is determined by many factors, including host species distributions and environmental parameters (LaJeunesse et al., 2004). To determine Symbiodinium biogeography, it is critical to use DNA markers with the resolving power to delineate between and among species (LaJeunesse & Thornhill, 2011). Our ability to discern different lineages of Symbiodinium has improved in tandem with the resolution of the molecular markers used to identify them. Originally, highly divergent lineages or "clades" designated A-I were confirmed via the sequencing of 18S and 28S ribosomal DNA (Rowan & Powers, 1991;Pochon & Gates, 2010), followed by the delineation of numerous subclades via phylogenetic analyses of the internal transcribed spacer regions (e.g., internal transcribed spacer 2 (ITS2); LaJeunesse, 2001). More recently, sequences from domain V of the chloroplast large subunit (cp23S) ribosomal DNA (Santos, Gutierrez-Rodriguez & Coffroth, 2003;Kirk et al., 2009;LaJeunesse & Thornhill, 2011;LaJeunesse, Parkinson & Reimer, 2012) and the non-coding region of the plastid minicircle (psbA ncr ; Takishita et al., 2003;LaJeunesse & Thornhill, 2011) have provided even greater resolution, often corresponding to the species level (LaJeunesse, Parkinson & Reimer, 2012). Microsatellite markers have also been developed to investigate fine-scale diversity in Symbiodinium (Pettay & LaJeunesse, 2007;Grupstra et al., 2017).
As a result of the development and implementation of finer resolution molecular markers, it is apparent that the extent of Symbiodinium diversity may be much greater than previously recognized based on studies using only internal transcribed spacer regions. For example, in examinations of Symbiodinium within the common reef zoantharian Palythoa tuberculosa (Esper, 1805) across a latitudinal gradient of 800 km in the Northern Red Sea, psbA ncr sequences revealed up to four unique lineages, despite ITS2 results showing only one single "subclade" lineage (Reimer et al., 2017). The distribution of these individual psbA ncr lineages strongly correlated with sea surface temperature (SST) differences of approximately 1 C. Such results demonstrate a need to re-examine previously reported Symbiodinium diversity with high-resolution markers, as fine scale, ecologically important differences may have been missed.
Okinawa-jima Island is ∼1,200 km 2 in area, just over 100 km in length, and 3-30 km in width. Although it is not "large" when compared to oceanic or regional scales, the surrounding marine environment encompasses a variety of ecosystems including large bays, muddy tidal flats, patch reefs within lagoons, and fringing reefs (Fujita et al., 2015). Additionally, recent population genetic studies on scleractinian corals (Shinzato et al., 2015), sea cucumbers (Soliman, Fernandez-Silva & Reimer, 2016), and amphipods (White, Reimer & Lorion, 2016) have described unexpected genetic structure among locations around Okinawa-jima Island, particularly between the Kuroshio-influenced west coast and the more isolated east coast.
Based on ITS2 sequence analyses of Symbiodinium within P. tuberculosa in Southern Japan (including Okinawa), we previously reported that subclade C1 or closely related types are dominant (Reimer, Takishita & Maruyama, 2006). Here we use psbA ncr sequences to re-examine how much variation exists within Symbiodinium C1 associating with P. tuberculosa, specifically focusing on Okinawa-jima Island's shallow coral reef environments. We additionally explore potential associations between observed diversity and the environment. This study is intended to complement our recent work in the Red Sea (Reimer et al., 2017) by focusing on the same host species and describing the extent of symbiont diversity on a much smaller geographical and environmental scale.

Environmental data
Satellite-derived SST and chl-a data for the waters around Okinawa-jima Island and nearby Amami Oshima Island were acquired from National Aeronautic and Space Administration Giovanni database (https://giovanni.gsfc.nasa.gov/giovanni/; Acker & Leptoukh, 2007), developed and maintained by the NASA Goddard Earth Sciences Data and Information Services Center. Error ranges of the moderate resolution imaging spectroradiometer (MODIS) Aqua data were approximately ±0.25 C for SST, and ±40% for chl-a. Yearly average SST (SST avg ) and chl-a data and maps used in this study were derived from 4 km resolution data from the MODIS Aqua database. These generated maps provided the basis for estimating annual average SST and chl-a at each sampling location. Data from February 2000 to May 2015 were used for SST avg analyses, and from July 2002 to May 2015 for chl-a. To examine yearly winter minimum SST (SST min ) and summer maximum SSTs (SST max ), we averaged monthly data from February and from August, respectively (2000-2014, n = 14 each).

Specimen collection
The zooxanthellate cnidarian species P. tuberculosa is the most common zoantharian on coral reefs surrounding Okinawa-jima Island (Irei, Nozawa & Reimer, 2011), and is easily identifiable (Hibino et al., 2014). Previous work in Southern Japan has also ascertained that the taxon does not contain cryptic species (Reimer et al., 2007), making it a simple and reliable species for research and citizen science (Parkinson et al., 2016).
Specimens of P. tuberculosa were collected from January 2012 to November 2015 from eight locations around Okinawa-jima Island and one location on Amami Oshima Island, Kagoshima, to the north of Okinawa in the Middle Ryukyus (Table 1). All specimens were collected from the low intertidal zone (0-2 m depths depending on tides) via snorkeling. Small portions of colonies were collected and fixed in 70-99.5% ethanol for further molecular analyses. The collections were limited by the low numbers of P. tuberculosa present at some sites (n = 3-11).
DNA extraction, PCR, and phylogenetic analyses DNA was extracted from preserved colony samples using a DNeasy Blood and Tissue Kit (Qiagen, Tokyo, Japan) following the manufacturer's protocol. We amplified two Symbiodinium DNA markers; the ITS2 in the ribosomal DNA array, and a portion of the plastid minicircle non-coding region (psbA ncr ). ITS2 sequences were obtained to place our new specimens within the phylogenetic framework of this well reported marker, and with past research on Symbiodinium within P. tuberculosa in Southern Japan (Reimer, Takishita & Maruyama, 2006), while psbA ncr sequences were obtained to examine finer scale phylogenetic patterns (Reimer et al., 2017). ITS2 was amplified using the primers zITSf (5′-CCG GTG AAT TAT TCG GAC TGA CGC AGT-3′) and ITS4 (5′-TCC TCC GCT TAT TGA TAT GC-3′) (White et al., 1990;Rowan & Powers, 1992;Hunter, Morden & Smith, 1997). psbA ncr was amplified using the primers 7.4-Forw (5′-GCA TGA AAG AAA TGC ACA CAA CTT CCC-3′) and 7.8-Rev (5′-GGT TCT CTT ATT CCA TCA ATA TCT ACT G-3′) (LaJeunesse & Thornhill, 2011). Reaction mixes contained 1.0 ml of genomic DNA, 7.0 m1 of Milli-Q water, 10.0 ml of HotStarTaq Plus Master Mix, and 1.0 ml of each primer (10 pmol). Thermocycler conditions were as follows: for ITS2: 95.0 C for Table 1 Sites from which P. tuberculosa specimens were collected in this study to examine Symbiodinium spp., and information on numbers of specimens, sea surface temperature (SST), and chlorophyll-a (chl-a) concentrations. The nucleotide sequences for ITS2 and psbA ncr were separately aligned within Geneious v9.1.3 (Biomatters Limited, Auckland, New Zealand). Alignments were inspected manually, and primer regions and short sequences were excluded. Because the long plastid non-coding region rarely sequenced completely, we used only the forward psbA ncr reads. The ITS2 alignment contained 24 sequences of 216 bp, while the psbA ncr forward alignment contained 61 sequences of 300 bp. Previously reported sequences from GenBank were incorporated into the ITS2 alignment for reference (DQ480631, DQ480639, DQ889741, DQ889743-all Symbiodinium subclade C1 or closely related from P. tuberculosa from Southern Japan; and AB207184-Symbiodinium subclade C15 related from Zoanthus sp. in Southern Japan), while the psbA ncr alignment contained only novel sequences (no previously reported sequences bore strong similarity).
Both alignments were analyzed using maximum likelihood (ML), neighbor-joining (NJ), maximum parsimony (MP) and Bayesian inference (BI) methods. ML analyses for both datasets were performed using PhyML (Guindon et al., 2010) with input trees generated by NJPlot (Perriere & Gouy, 1996) under automatic model selection by smart model selection with Akaike Information Criterion. Both datasets were analyzed under the HKY85 model (Hasegawa, Kishino & Yano, 1985) with the transition/transversion ratio estimated, the proportion of invariable sites fixed at 0.0, and the number of substitution rate categories as 1. PhyML bootstrap trees were made using the same parameters as the individual ML tree. The distances were calculated using a Kimura's twoparameter model (Kimura, 1980). NJ analyses for both alignments were run within Geneious on default settings under the HKY85 model. MP analyses were performed in Paup Ã 4.0a147 (Swofford, 2000) with indels included as a fifth character state. All trees were run with 1,000 bootstraps. Bayesian posterior probabilities were calculated with the software Mr. Bayes (Huelsenbeck & Ronquist, 2001) using the HKY85 substitution model and default parameters (chain length = 1,000,000; burn-in = 250,000). Genetic distances between and within lineages were calculated in MEGA6 (Tamura et al., 2013) using the Maximum Composite Likelihood Model (Tamura, Nei & Kumar, 2004).

Environmental data
Yearly average SST showed southern sites to be warmer than northern sites, with a difference of 0.95 C between Wase on Amami Oshima Island (24.6 C) compared to Kyan and Odo on the southern tip of Okinawa-jima Island (25.35 C) (Table 1). For SST max , Wase was lowest (28.88 ± 0.93 C) while Uken on the east coast of Okinawa-jima Island was highest (29.43 ± 0.76 C) (Table 1). For SST min , Wase (20.35 ± 0.43 C) was coldest, with highest SST min at Odo (21.42 ± 0.60 C). The highest observed SST in any year was at Uken (30.9 C in 2001), and the lowest was at Bise (18.7 C in 2015) and Oku (18.7 C in 2008) on the northwest and north coasts of Okinawa-jima Island, respectively (Table 1). Yearly average chl-a concentration values were generally low at all sampling sites, ranging from a low of 0.08 mg/m 3 at Odo to a high of 0.50 mg/m 3 at Uken (Table 1).
Finally, we examined the range of environments in which each lineage could be found (Fig. 3). Symbiodinium lineage 1 appeared at all sites and thus all environments in this study (SST avg = 24.6-25.35 C; SST max = 29.43 C; SST min = 20.35 C; chl-a 0.08-0.50 mg/m 3 ). Lineage 2 was only observed at Oku, Teniya, and Mizugama, where SST avg ranged between 24.95 and 25.15 C, with SST min of 20.69 C (Oku) and SST max of 29.16 C (Teniya, Mizugama), and chl-a ranged between 0.08 and 0.25 mg/m 3 . Lineage 3 was only found at Bise and Nerome, with SST avg of 24.85-24.95 C, SST min of 20.47 C, SST max of 29.33 C (both Nerome), and chl-a of 0.15-0.30 mg/m 3 . Lineage 4 members were only observed once each at three locations, but these stretched across the geographic range of this study (Wase, Mizugama, Odo).

Notes:
Shaded diagonal values represent within-lineage distances. * Lineage 4 was represented by two identical sequences.

DISCUSSION
Using the high-resolution psbA ncr marker, we identify four Symbiodinium lineages associated with P. tuberculosa on Okinawa-jima Island. The lineages feature surprisingly unique distributions over a small geographic scale, and would be considered at most two entities based on lower-resolution ITS2 data (Fig. 1). Because the between-lineage molecular distances for psbA ncr (0.105-0.256; Table 2) are greater than those reported between the psbA ncr sequences of two divergent Symbiodinium ITS2 types (∼0.045 for C26a vs. C31; LaJeunesse & Thornhill, 2011), these lineages likely represent reproductively isolated species rather than populations within a species. Lineage 1 appears to be a widely distributed generalist, at least over the range of this study in the Central Ryukyus (Fig. 2). It is found at all sites and often occupies the majority of colonies at a given site. Lineage 3 is observed only at the two sites on the northwestern coast of Okinawa-jima Island, and thus appears to have narrow geographical and environmental components to its distribution (SST avg = 24.85-24.95 C; chl-a = 0.15-0.30 mg/m 3 ) (Fig. 3). Lineage 4 is present in very low numbers (n = 3) across the latitudinal range of the study, making it difficult to infer its environmental niche. A B Figure 2 Map of Amami Oshima Island (A) and Okinawa-jima Island (B) with average sea surface temperature (SST avg ) and Symbiodinium psbA ncr lineage ratios at each site investigated. Note thermal distortions near coastlines were ignored in all SST analyses as these are generated by influence of terrestrial portions of islands within the 4 km resolution of satellite data.
Lineage 2 has a somewhat restricted range, found at only three locations (where it also occupied the majority of colonies): Oku, Teniya, and Mizugama. These three sites are located near river outflows, suggesting this Symbiodinium lineage may be able to tolerate changes in salinity more effectively than the others. Unfortunately, fine-scale salinity data are not yet available to investigate this trend further. Chl-a levels at these sites are generally low (<0.25 mg/m 3 ), while the SST avg range is intermediate (24.95-25.15 C). Lineage 2 is absent at other locations within this SST avg range; for example, at Bise and Nerome, where only lineages 1 and 3 are present, and at Uken, where only lineage 1 is detected. Although Uken's SST avg is firmly in the middle of the SST range investigated in this study, Uken's SST min and SST max are generally more extreme than those at other sites (Table 1) due to Uken's position within shallow Southern Kin Bay, which is isolated from the stabilizing temperature effects of the open ocean (Montani, 1996). Additionally, chl-a levels at Uken are higher than those of all other locations (0.50 mg/m 3 ).  By focusing on a very small area of the Northwest Pacific, and by using a rapidly evolving molecular marker, we could resolve a much finer scale of Symbiodinium biogeography than has previously been recognized in the region. Some psbA ncr lineages appear to be partitioned on the basis of temperature (SST max ) differences on the scale of 0.1-0.3 C, much lower than the 0.7-1.0 C observed in previous studies (Baums, Devlin-Durante & LaJeunesse, 2014;Reimer et al., 2017). It also appears SST stability (e.g., differences between SST max and SST min ), as well as fine scale salinity differences (not measured here) may play an important role in the distribution of Symbiodinium diversity. Chl-a levels are generally low at all locations, although it should be noted that only generalist lineage 1 is found in Uken (although n = 3), the site with the highest chl-a levels.
Overall, the contribution of chl-a (a proxy for turbidity) is not clear, as all P. tuberculosa specimens were collected from very shallow waters (0-2 m).
It is surprising that such biogeographic patterns could be uncovered given the low numbers of host P. tuberculosa at some locations. However, these low numbers also caution that the symbiont species' distributions are unlikely to have been completely resolved. Another issue is that the fine-scale environmental variation among sites falls below the margin of error in the satellite datasets; future research of this nature must utilize more precise instrumentation. We did not examine host P. tuberculosa population genetics, and it remains to be seen if fine-scale host structure may also play a role in the observed patterns. Previous studies indicate cnidarian host and symbiont genetic structure can be mismatched (Baums, Devlin-Durante & LaJeunesse, 2014;Leydet & Hellberg, 2016), although clear cases of matching genetic structure have also been observed (Bongaerts et al., 2010;Prada et al., 2014).
As our ability to discern between different lineages of Symbiodinium has increased, so too has our understanding of the complexity and nuances of Symbiodinium diversity and distribution. The present work supports previous studies that show Symbiodinium evolution is driven largely by specialization to different environmental niches (Frade et al., 2008b;LaJeunesse et al., 2010;Kamezaki et al., 2013;Reimer et al., 2017), and that specialization may occur on much finer micro-environmental scales than usually addressed. Further characterization of the Symbiodinium-P. tuberculosa symbiosis within the Ryukyus or comparable island chains should be carried out to confirm similar symbiont structuring based on fine-scale environmental heterogeneity.
The Symbiodinium diversity patterns on Okinawa-jima Island highlight three major considerations for future investigations of this kind. First, it is advantageous for researchers to obtain fine-scale environmental data of the study area if available, so as to better characterize niches that might otherwise be deemed homogenous. Second, researchers should refrain from considering pooled specimens or results from different nearby locations as representative of a larger area without thoroughly checking for fine scale environmental patterns. Third, as has been suggested elsewhere (LaJeunesse & Thornhill, 2011;Reimer et al., 2017), researchers should incorporate both ITS2 data (to tie to past work) as well as rapidly evolving markers (to uncover hidden diversity) when investigating Symbiodinium biogeography, as a failure to address fine-scale niche adaptation could lead to a misinterpretation of results. These suggestions