Strong horizontal and vertical connectivity in the coral Pocillopora verrucosa from Ludao, Taiwan, a small oceanic island

Mesophotic habitats could be sheltered from natural and anthropogenic disturbances and act as reproductive refuges, providing propagules to replenish shallower populations. Molecular markers can be used as proxies evaluating the connectivity and inferring population structure and larval dispersal. This study characterizes population structure as well as horizontal and vertical genetic connectivity of the broadcasting coral Pocillopora verrucosa from Ludao, a small oceanic island off the eastern coast of Taiwan. We genotyped 75 P. verrucosa specimens from three sites (Gongguan, Dabaisha, and Guiwan) at three depth ranges (Shallow: 7–15 m, Mid-depth: 23–30 m, and Deep: 38–45 m), spanning shallow to upper mesophotic coral reefs, with eight microsatellite markers. F-statistics showed a moderate differentiation (FST = 0.106, p<0.05) between two adjacent locations (Dabaisha 23–30 and Dabaisha 38–45 m), but no differentiation elsewhere, suggesting high levels of connectivity among sites and depths. STRUCTURE analysis showed no genetic clustering among sites or depths, indicating that all Pocillopora individuals could be drawn from a single panmictic population. Simulations of recent migration assigned 30 individuals (40%) to a different location from where they were collected. Among them, 1/3 were assigned to deeper locations, 1/3 to shallower populations and 1/3 were assigned to the right depth but a different site. These results suggest high levels of vertical and horizontal connectivity, which could enhance the recovery of P. verrucosa following disturbances around Ludao, a feature that agrees with demographic studies portraying this species as an opportunistic scleractinian.

In most benthic marine species, recruitment and settlement of new individuals are achieved through pelagic larval movement [9]. While larval dispersal and the effective recruitment at shallow and mesophotic depths are critical to the reef dynamic [33], quantifying them remains a challenge [34] considering that benthic larvae have limited mobility and small recruits. Although biophysical models predict that such movements across depth significantly contribute to demography and resilience in reef coral populations [15], genetic connectivity approaches offer a valuable proxy for assessing population structure and estimating larval flux across depths [6,12,[35][36][37][38][39][40][41][42][43][44][45]. To date, genetic approaches have revealed contrasting patterns of connectivity between shallow and deep populations of reef corals, suggesting that connectivity is taxon-and location-specific [1,46]. However, despite a growing body of literature on this topic, research evaluating connectivity between shallow and mesophotic coral populations is still considered to be preliminary [1,47].
In Taiwan, well-developed fringing reefs have been observed around Ludao, a small volcanic island 33 km off the coast of Taitung. Along the bathymetric gradient, Pocillopora verrucosa is one of the most dominant scleractinian hard corals, distributed from the near-surface to as deep as 55 m [48]. Despite variation in corallum micro and macro morphologies with depth [49], sequencing of the mitochondrial Open Reading Frame marker (mtORF, from a gene with unknown function) has confirmed this species' distribution across a large bathymetric gradient in Ludao [50]. Several microsatellite markers (also called Simple Sequence Repeats, SSRs) have been developed for Pocillopora species [51][52][53][54] and were used to investigate species diversity and population connectivity across their biogeographic distribution [55][56][57][58]. However, to date little has been done to assess connectivity across depths, population structures, or larval fluxes between shallow and deep populations of Pocillopora verrucosa. The present study aims to make these assessments.

Sampling and DNA purification
In 2016 and 2017, Pocillopora verrucosa specimens were collected at three sites around Ludao -Gongguan, (22.683530˚N, 121.496631˚E), Guiwan, (22.639212˚N, 121.483529˚E), and Dabaisha (22.637706˚N, 121.490844˚E)-using compressed air or Trimix diving, with a minimum distance of 3 m between samples to avoid collecting colonies that might be clones (Fig  1). A collection permit was obtained from the Taitung County Government (1040000285) and all specimens were registered at the Biodiversity Research Museum (Academia Sinica, Taiwan). Distances between sites are < 6 km (Fig 1).
DNA extraction process followed De Palmas et al. [50]. Additionally, DNA extracts were purified and concentrated using a sodium acetate (3M, pH 5) precipitation method. For each DNA sample, 1/10 volume of sodium acetate was added, and mixtures were gently inverted for 10 min. Three volumes of pure ethanol were added, and solutions were kept at -20˚C overnight. Samples were centrifuged at 14,000 g for 30 minutes and DNA pellets were rinsed three times with -20˚C ethanol. The DNA pellets were re-suspended in 30 μL TE buffer (1X, USB Corporation, Cleveland, OH, USA). DNA concentration and quality were measured on a Nanodrop 2000 (Thermo Fisher Scientific) and diluted aliquots were prepared to homogenize the DNA concentration of each sample before Polymerase Chain Reaction (PCR).

Microsatellite marker selection
Microsatellite marker effectiveness was determined by amplifying eight samples haphazardly picked with 30 microsatellite markers, accounting for most microsatellites developed for Pocillopora species. Each PCR mixture contained 15 μL of Master Mix RED (Ampliqon, Odense M, Denmark), 15 μL of ddH2O, 2 μL of each pair of primers (2.5 μM), and~50 ng of DNA template. PCR conditions were as follows: 60 s of denaturation step at 94˚C; 40 cycles of 60 s at 94˚C, 60 s at the annealing temperature (Ta), and 60 s at 72˚C; and an extension step at 72˚C for 5 min. PCR amplifications were checked on a 2% agar gel to confirm amplification success rate and product size. Successful amplifications were diluted 22X before processing on a 5400 Fragment Analyzer™ Automated CE System (Agilent, CA, USA) using a dsDNA 905 reagent kit (1 to 500 bp). Any markers displaying low PCR amplification success or non-specific PCR amplification (multiple additional PCR products) were removed from subsequent analyses.

Microsatellite screening and allele scoring
Loci that passed the selection process were amplified to genotype the 85 P. verrucosa. Each PCR product was run three times on the fragment analyzer (as described above) to account for the accuracy of the fragment analyzer (i.e., 2-3 bp, per the manufacturer's instructions). The three repeated electropherograms were scored using Prosize 3.0 (Agilent, CA, USA) and allele repeat size was then estimated using GenAlEx 6.5 [59] based on the three scored electropherograms. The estimation of repeat size and repeated allele scoring were used to solve allele scoring ambiguity. Micro-Checker 2.2.3 [60] was used to detect and correct eventual scoring inconsistencies and scoring errors due to stuttering and large allele drop-out. Samples showing low standard scoring quality were removed from the analysis.

Computational analysis
Arlequin 3.5.2.2 [61] was used to calculate the total number of alleles (Na), the observed (Ho) and expected (He) heterozygosities, and the exact probability of departure from Hardy-Weinberg equilibrium for each locus (pHWE) with 10 6 steps in a Markov chain and 10 5 dememorization steps. Arlequin was also used to investigate linkage disequilibrium between all pairs of loci with 5 x 10 5 permutations. Null allele frequency was inferred using MicroChecker 2.2.3 [60] with a 95% confidence interval and 10 5 randomizations.
Genetic differentiation among all locations was estimated using a pairwise F ST comparison at all loci and their statistical significance with 5 x 10 5 permutations in Arlequin. Pairwise F ST were also calculated using FreeNa software [62] with the number of replicates fixed at 50,000. Unlike Arlequin, FreeNa takes into consideration the presence of null alleles and provides an accurate estimation of F ST values [62]. To account for the low number of samples per location, pairwise F ST were also computed at the site level (pooling 3 depth ranges at one site together) and at the depth level (pooling similar depth ranges from all 3 sites together).
Genetic structuring was assessed using STRUCTURE 2.3.4 [63] with 10 6 Markov Chain Monte Carlo simulations, a burn-in period of 10 5 and three iterations for each K from k = 1 to k = 12. The LOCPRIOR option was used (to account for the location characteristics) under the admixture model and assuming independent allele frequency. The web-based program STRUCTURE HARVESTER 0.6.94 [64] was used to estimate the K number of populations and compare them to the web-based program CLUMPAK [65]. A principal coordinate analysis (PCoA), which summarizes genetic distances between individual multi-locus genotypes to cluster individuals/locations relative to each other, was performed at the location level using GenAlEx 6.5 [59].
In order to detect individuals with a recent migration history, we use the "Detection of first-generation migrants" option from Geneclass 2 [66] which assigns potential new migrant to their most likely source population [67]. This option calculates the likelihood L home /L max (with L home the likelihood computed from the population where the individual was sampled and L max the highest likelihood value among all population including the population where the individual was sampled) using the frequencies-based method from Paetkau et al. [68] and Monte-Carlo resampling algorithm with 10,000 simulated individuals and a type I error threshold α = 0.05 [67]. To account for the low number of samples per locations, this simulation was also computed at the site (pooling 3 depth ranges from one site together) and at the depth level (pooling similar depth range from all 3 sites together).

Results
Of 85 P. verrucosa specimens genotyped in this study, 75 produced satisfactory amplification and scoring results (� 4% of data missing) for eight microsatellite markers, Psp_16, Psp_23, Psp_29, Psp_32, Psp_35, Psp_39, Psp_41, and Psp_48 (Table 1). Full dataset with allele scoring, mtORF haplotype and museum numbers is provided in S1 Table. Loci were polymorphic, but not for all combinations of location and marker (48 combinations in total): Psp_29 was monomorphic at Dabaisha D and Gongguan S; Psp_35 was monomorphic at Dabaisha M, Gongguan S, and Gongguan M; and Psp_39 was monomorphic at Gongguan M ( Table 2). The number of alleles per location and locus ranged from one to 19 (Table 2). Twelve combinations of location and loci deviated from Hardy-Weinberg equilibrium (Table 2), representing 16% of the total combinations. Linkage disequilibrium was found in Psp_29/Psp_39, and Psp_32/Psp_48. Null alleles were detected in Psp_29, Psp_35, Psp_39, Psp_41, and Psp_48. No clone was detected.
Except for the Dabaisha M and Dabaisha D comparison, which displayed significant differentiation (F ST = 0.106, p < 0.05, Fig 2), all F ST values were low (0.000 < F ST < 0.058) and not significant. Pairwise F ST computed at the site and depth levels were also low (0.000 < F ST < 0.007) and non-significant (S1 File). In addition, F ST values that were corrected for the presence of null alleles were not significantly different from non-corrected ones. The PCoA did not reveal any site-specific or depth-related pattern (Fig 3).
Assignment tests yielded 40% of our specimens (30/75) as potential new migrants: 10 came from deeper locations, 10 came from shallower locations and the remainder relocated from similar depth but from a different site (Table 4). This result was consistent when specimens were grouped by site (10 potential migrants) and depths (10 potential migrants, including 3 detected as coming from shallower locations and 7 detected as coming from deeper locations; S1 File).

Discussion
Pocillopora verrucosa colonies originating from three sites and three depth ranges were genotyped at eight microsatellite loci to evaluate genetic differentiation and genetic structure, and to estimate migration flux between sites and across depths. Overall, P. verrucosa exhibits a low level of genetic differentiation, with the exception of Dabaisha M and Dabaisha D, which showed a mild genetic differentiation despite being separated by a short horizontal distance. Pairwise F ST values generally support high genetic connectivity between sites and across depths (Fig 2 and S1 File), and the absence of clustering suggested that P. verrucosa specimens were collected from the same, well-mixed population unit (Table 3). We further identified possible Table 1. Characteristics of microsatellite markers tested in this study. Initial PCR amplification success was evaluated for a subset of eight samples. Loci displaying PCR amplification success >75% and expected fragments size for the initial subset (in bold) were used to genotype the 85 Pocillopora verrucosa or otherwise excluded from this study.

PLOS ONE
migrants originating from sites/depths different from those where they were collected (Table 4 and S1 File), highlighting that MCEs could act as a source of propagules for nearby shallow water communities (and vice versa). Our data indicate that populations of P. verrucosa around the oceanic island of Ludao may approximate panmixia, thus we hypothesized that the mesophotic zone around this island could 1) act as an ecological refuge against disturbances for several other depth generalist species inhabiting this area and 2) benefit from the Kuroshio Current for dispersion of coral larvae and increase the overall connectivity along its path.

Low genetic differentiation and absence of genetic structure across depth in P. verrucosa
In this study, P. verrucosa from Ludao exhibited low genetic differentiation between sites located less than 6 km apart, and across depth zonations occurring at shorter horizontal distances. Only one pairwise comparison showed weak but significant differentiation (F ST = 0.106), which corresponded to two adjacent locations at Dabaisha (M: 23-30 m and D: 38-45 m, Fig 3). Scleractinian coral populations may exhibit strong genetic breaks between shallow and mesophotic depths, even when separated by a short horizontal distance (< 1 km) [12,39,43,70]. High levels of genetic differentiation could be linked to biotic factors such as selection against migrants (or genotypes) in their non-native environment [12], temporal differences in reproductive timing [71], and/or depth-dependent parental effects [72]. Alternatively, high levels of genetic differentiation could also be linked to abiotic factors such as contrasting reef geomorphology and water currents between shallow and mesophotic habitats that create a physical barrier to dispersal [39,70]. In Ludao, there are evidences for shallow and mesophotic reefs to produce recruits in synchrony [73]. For Pocillopora, similar molecular operational taxonomic units (MOTU) are recruiting at both shallow and mesophotic depths [73]. In addition, few scleractinian corals from the mesophotic zone have been observed to have gametes, in synchrony with shallow water corals, affirming that these scleractinians are reproductively active (Fymbriaphyllia ancora: -45 m and Acropora spp.: -40 m, De Palmas S, personal observation).
Despite the short horizontal distance separating Dabaisha M from Dabaisha D (< 100 m), both locations are well-known for their strong currents. The slight difference between these two adjacent locations could be due to diverging currents generating a physical barrier to P. verrucosa larval dispersal. Overall, low genetic differentiation and the absence of any genetic Table 3. Estimates and probabilities of K population structure. Mean LnP (K), Ln' (K), Ln" (K), and Delta K are given using the Evanno method (see [69] for details). Prob (K = k) is the probability for each K given by STRUCTURE. � computed in CLUMPAK. PLOS ONE Table 4. Detection of potential new migrant. For each specimen (Assigned sample) from the location where it was sampled (Origin) the table gives the likelihood (-log (L home /L max ) of this specimen to be a new migrant with its respective probability p < α and the likelihood of originating from another location (-log(L)). Shaded lines highlight potential new migrants their possible location of origin in bold. Number of loci (Nb. of loci) are the loci used for calculations (loci with missing data were not included into the calculation).  structures within and between sites further suggest that P. verrucosa specimens originated from a single population, and that this population was probably panmictic. Previous investigations have found contrasting patterns of connectivity between shallow and mesophotic populations, varying from open to closed populations, with connectivity patterns occurring over thousands of kilometers [39,43,70,[74][75][76]. For instance, populations of the coral Montastraea cavernosa showed genetic differentiation with depth in Florida [39] and Little Cayman Island [41], but not in Bermuda or the US Virgin Islands [39,41]. Similar findings were seen in Seriatopora hystrix from the Great Barrier Reef: the species displayed genetic differentiation with depth at Yonge Reef and Day Reef [36], but not at Scott Reef [40]. Interestingly, while it was proposed that broadcast spawners (i.e., species that release gametes into the water column and undergo external fertilization) may have higher dispersal abilities than brooders (i.e., species that release sperm and larvae into the water column and undergo

PLOS ONE internal fertilization) [36], low genetic differentiation was observed in the broadcasting species
Stephanocoenia intersepta [12]. In addition, divergence by depth has also been reported in both brooders and spawners [36, [39][40][41]43]. Pocillopora verrucosa is a broadcasting spawner, which spawns in the early morning [77], 2-3 days after the full moon in April and/or in May in Taiwan [78]. Mulla et al. [78] found larvae of P. verrucosa to display positive photo-movements, i.e., phototaxis or photokinesis. They further suggested that P. verrucosa larvae could disperse further and benefit from stronger surface currents as they move upward into the water column. This larval behaviour could explain the absence of genetic structure observed in the present study. In general, the ability of mesophotic populations to reseed shallow water reefs should be regarded as species-and location-specific [1] with local hydrology playing a key role in disseminating coral larvae [15]. In Ludao, mesophotic and shallow zonations have large scleractinian diversity overlap (37%, [79]), so that many other reef coral species with large bathymetric distribution may contribute to the reproductive refuge.

Possible new migrants and their origins
Effective ecological connectivity between shallow and deep zonations is achieved under three prerequisite criteria: 1) suitable habitats supporting robust populations, 2) habitats that are physically connected via water or propagule movements, and 3) organisms with life-history traits that enable successful replenishment of both their own populations and those connected to them [45]. In Ludao, Denis et al. [48] and De Palmas et al. [50] reported that P. verrucosa is one of the dominant species found across the depth gradient. Our data strongly suggest that larval movement creates connectivity between shallow and deep zonations. Our simulation test detected potential migrants between sites and across depths, which reinforces the notion that mesophotic populations can support themselves as well as adjacent populations. About 40% of our specimens were assigned to different populations from where they were collected (Table 3). In this study, loci and sample numbers were relatively low but our methodology produced a similar outcome when individuals were grouped by sites and depths (S1 File). Among the 30 individuals detected as potential migrant, 1/3 were assigned coming from deeper locations, 1/3 were assigned coming from shallower locations and the rest were assigned to similar depth range from different sites. This result suggests strong exchange of individual between sites and depths, showing that both vertical and horizontal larval exchange shape the populations of P. verrucosa at Ludao. Moreover, this is congruent with the low observed genetic differentiation and absence of population structure between sites and across depths, and suggests that P. verrucosa disperses effectively around Ludao. Assignment methods based on the analysis of genetic data should still be taken with caution due to the relative limitations induced by the presence of null alleles and the sensitivity of such methods to imperfect data (missing data and/or an unbalanced number of loci) [66]; however, concordant evidence suggests that P. verrucosa shows effective ecological connectivity around the island. Our study design includes adult colonies that are probably different ages, so the connectivity pattern is likely the result of multi-generational larval exchange, rather than associated with a specific P. verrucosa cohort. Population genetics comparing adult colonies to newly recruited P. verrucosa could tackle this issue and allow researchers to evaluate the current populations' dynamics.

Deviation from Hardy-Weinberg, linkage disequilibrium, and null alleles
Among all the combinations of sites, depths, and loci, 12 (out of 96) deviated from Hardy-Weinberg equilibrium (Table 1). In addition, linkage disequilibrium-i.e., non-random association of alleles-was found in Psp_29/Psp_39, and Psp_32/Psp_48. These markers were not reported to have linkage disequilibrium in Pocillopora populations from the Ryukyus for which they were developed [53]. Linkage disequilibrium can be a result of evolutionary forces such as natural selection or genetic drift [80]. However, null alleles could also slightly affect linkage disequilibrium estimators [81]. Microsatellite null alleles occur when flanking regions have accumulated enough mutations to prevent primers from annealing during PCR amplification [62]; they could also be the result of preferential amplification of short alleles (i.e., Short Allele Dominance, SAD) [82]. Consequently, some alleles simply do not amplify during PCR, in which case individuals that are heterozygous will appear to be homozygous following fragment analysis [83]. In this study, five loci (Psp_29, Psp_35, Psp_39, Psp_41, and Psp_48) were detected as null alleles, but no SAD effect was found at these loci, suggesting that the flanking regions of those microsatellite loci could have accumulated enough variation to prevent marker annealing during PCR. As a deviation from Hardy-Weinberg equilibrium results from a shortage in heterozygotes (or excess in homozygotes), we hypothesize that the high number of potential null alleles found in this study could influence both linkage disequilibrium and deviation from Hardy-Weinberg equilibrium. Null alleles can also influence population differentiation estimators [68], so they require correction for the potential bias they induce, but in our case, no differences were found between biased and corrected F ST . Null alleles are common in a large variety of taxa and, interestingly, are often found in taxa with a large population [62]. P. verrucosa is indeed locally abundant in Ludao, as it is one of the major components of the scleractinian community, residing from very shallow to at least 55 m deep [48,50]. Population censuses for scleractinian corals inhabiting the mesophotic zone could reveal important elements for evaluating the dynamic of their populations and in fine the role of the mesophotic zone as an ecological refuge.

Cross-specific amplification for Pocillopora microsatellite markers
The microsatellite markers used here were developed for a variety of Pocillopora species (P. damicornis [54,84], P. verrucosa, P. meandrina and P. damicornis [51,52], P. acuta, and P. meandrina [53]), and are known to amplify closely related taxonomic groups [85][86][87]. However, our analyses found that only eight markers produced both satisfactory amplification and scoring results for P. verrucosa from Ludao, and several others did not provide any exploitable results (low amplification success, low specificity, Table 1). Cross-specific amplifications of microsatellites among Pocillopora species display contrasting results when compared to previous work. Nakajima et al. [53] confirmed that 14 previously described Simple Sequence Repeat (SSR) markers can amplify Pocillopora types 1 (P. meandrina and P. eydouxi), 3 (P. verrucosa), 4 (P. damicornis), 5 (P. acuta), and 8 (an undescribed Pocillopora species), while several other SSR markers could not (PV3, PV5, PV6, Pd3-002, Pd2-003, Pd3-010, and Pd13). None of these previously developed markers produced satisfactory results in P. verrucosa from Ludao. For instance, PV3, PV5, PV6, Pd3-002, Pd2-003, and Pd13 had low amplification rates (<50%) and Pd3-010 did not amplify at all (Table 1). Markers PV2, PV5, PV6, and PV7 successfully amplified P. verrucosa specimens from Southeast Africa [88], and most of these markers generally amplify Pocillopora species from the Western Indian Ocean, the Southwestern Pacific, Southeast Polynesia [55,58,89], and the Red Sea [90]. Among the 13 markers developed by Nakajima et al. [53] for populations from Japan and predicted to amplify most Pocillopora species, four (Psp_1, Psp_2, Psp_10, and Psp_18) showed no or poor amplification (<25%), and an additional marker (Psp_33) was found to have major additional products (Table 1). We hypothesize that the low amplification success of most SSR markers on P. verrucosa in Ludao, and differences in cross-specific amplification, could be the result of the accumulation of important mutations in the microsatellite flanking regions [91], which may reflect strong differentiation between populations from Taiwan and those for which these markers were developed and tested. Verifying this hypothesis would require sequencing both microsatellite regions and their flanking regions or applying single nucleotide polymorphisms (SNPs), such as used in Taninaka et al. [92], to reveal population difference of P. verrucosa across large biogeographic scale in the future.

Conclusions
Our results demonstrate that P. verrucosa populations around Ludao from 7 m to as deep as 45 m are close to panmictic. This species is more likely to recover from acute disturbances occurring in shallow waters (typhoons or warm water bleaching events) through the recruitment of larvae originating from adjacent shallow and deep habitats. We think that the reproductive status and reseeding capacity of mesophotic corals in this area require further investigations. Ludao, as an oceanic island could act as an important source of recruits along the Kuroshio Current for populations located on its path [93]. However, to the best of our knowledge, no other study addresses genetic connectivity or larval flux between shallow and deep coral populations in Taiwan [48]. In fact, deeper coral ecosystems around Taiwan are generally not the focus of conservation plans [48], and this study calls for greater consideration to be paid to these ecosystems, as they may constitute important reproductive refuge on the path of the Kuroshio current.
Supporting information S1