Strong Genetic Structure and Limited Gene Flow among Populations of the Tropical Seagrass Thalassia hemprichii in the Philippines

: Seagrasses are marine angiosperms, and seagrass beds maintain the species diversity of tropical and subtropical coastal ecosystems. For proper understanding, management and conservation of coastal ecosystems, it is essential to understand seagrass population dynamics. Population genetic studies can cover large geographic scales and contribute to a comprehensive understanding of reproductive dynamics and potential dispersal among locations. The clonal and genetic diversity and genetic connectivity of Thalassia hemprichii in the Philippines were estimated by a population genetics approach. The geographic scale of this study has a direct distance of approximately 1600 km. Although high clonal diversity was found in some sites ( R = 0.07–1.00), both sexual and asexual reproduction generally maintains separate populations. Genetic diversity is not deﬁnitely correlated with latitude, and genetic differentiation is signiﬁcant in all pairs of sites ( F ST = 0.026–0.744). Complex genetic structure was found in some regions, even at a ﬁne geographic scale. The migration of fruits and seedlings was elucidated as an infrequent and stochastic event. These results suggest the necessity for the conservation of this species due to a deﬁciency in migrants from external regions.


Introduction
Recently, rapid anthropogenic disturbances of coastal ecosystems have been occurring at global and local scales [1,2].Such disturbances could cause the extinction of coastal organisms and affect many ecosystem services for humans, such as maintenance of tourism, fisheries, erosion control, water purification and carbon sequestration [3,4].Seagrasses are coastal marine angiosperms, and the species richness of seagrasses in Southeast Asia is the highest in the world [5].Estimating the geographic population structure and connectivity based on biogeographic research is essential for the proper understanding, management and conservation of coastal ecosystems.Estimates of the genetic composition and diversity at each location, as well as the genetic structure and connectivity among populations, are crucial for understanding the ecological characteristics of coastal organisms, including seagrasses.
Population genetic studies can encompass large geographic scales and contribute to a comprehensive understanding of the genetic structure and connectivity of seagrasses.Genetic studies have been conducted in some seagrass species, and significant genetic differentiation in seagrasses has been detected among regions based on large geographic analyses at distances of greater than 1000 km (e.g., Cymodocea nodosa [6], Posidonia oceanica [7], Thalassia testudinum [8], Zostera marina [9][10][11]).These analyses have suggested potential long-distance dispersal among regions.Long-distance dispersal is stochastic but important for maintaining connectivity among populations and intraspecific homogeneity.Nevertheless, population genetic studies of seagrasses have seldom been conducted in the Indo-West Pacific region, despite the high species richness in these regions.Recently, the isolation of microsatellites and population genetic studies have been conducted for seagrass species in the region (Enhalus acoroides [12][13][14], Cymodocea serrulata [15,16], Cymodocea rotundata [14,17,18] and Syringodium isoetifolium [19,20]).Overall, those population genetic studies of seagrasses showed high genetic differentiation among populations.
The genus Thalassia is composed of two species, Thalassia hemprichii and T. testudinum.Thalassia hemprichii, known as turtle grass, is found in the Indo-West Pacific region in the tropical and subtropical zones and is one of the most widely distributed and dominant seagrass species in the Northwest Pacific [5].Thalassia hemprichii is a representative species in Blue Carbon ecosystems for carbon storage in tropical and subtropical regions of the Indo-West Pacific [21,22].This genus is distributed over vast shallow coastal areas, and can spread by extending its underground rhizomes and forming dense meadows.Although flowering events of T. hemprichii occur at a local scale, meadows of this species are primarily maintained by vegetative propagation [23,24].van Dijk and van Tussenbroek [25] reported a high rate of clonal expansion by T. testudinum in the Caribbean Sea, with the longest genet estimated to reach up to 230 m.Nevertheless, there is no evidence of inbreeding via sexual reproduction due to the high extent of kinship within local populations of T. testudinum [26].Thalassia hemprichii fruits could theoretically travel 73.5 km over a maximum of 8 days [27].Successful germination of fruits floating for 19 days was reported in T. hemprichii [28].Potential long-distance dispersal of fruits in T. testudinum ranged from approximately 190 to 720 km over 10 days, and this long dispersal distance appears to maintain the gene flow among sites (approximately 350 km) [8].The dispersal period is shorter than E. acoroides [27] and longer than Halophila ovalis [28].However, the dispersal cannot be observed and verified in actual seagrass meadows.
In this study, microsatellite markers were used to investigate the clonal diversity, genetic diversity, genetic structure and connectivity of populations of T. hemprichii in the Philippines to determine the reproductive dynamics and potential dispersal mechanisms among sites.Long-distance dispersal of fruits by ocean currents may contribute to maintain the genetic connectivity among sites in T. hemprichii.However, if long-distance dispersal is rare and stochastic, genetic differentiation and structure will be confirmed among sites.In this study, we verify whether T. hemprichii maintains connectivity by gene flow owing to frequent dispersal or is genetically separated among sites due to rare and stochastic dispersal in that region.

Study Sites and Sample Collection
Seagrass sampling sites were located in the Philippines (33 sites; Figure 1, Table 1).Initial letters in the population codes indicate the region: B, Batanes; N, North Luzon; W, West Luzon; D, Mindoro; E, East Luzon; P, Palawan; V, Visayas; M, Mindanao (see Table 1).The maximum distance between sites, from Chanarian (B-CNR) to San Agustin (M-SAG), was approximately 1600 km.Major surface ocean currents are also shown in Figure 1, following the map by Ravago-Gotanco and Juinio-Meñez [29].A strong westward current, the North Equatorial Current, flows from the Pacific Ocean to the eastern coast of the Philippines, through the study region.The current bifurcates into two currents, and the northward-flowing branch, the Kuroshio Current, is a particularly strong current.The southward current, the Mindanao Current, flows toward the southern Philippines through the eastern coastal area.Other currents flow into the inland sea within the Philippine Islands, with southward currents flowing from the South China Sea to the Sulu Sea.We collected samples from approximately 200 m × 300 m sites in the intertidal or subtidal zones, which were less than approximately 3 m in depth.We maintained a distance of more than 10 m among shoot samples collected within a site to unify the sampling strategy with other studies in this region [13,16,18,20].A young leaf from each shoot was collected, dried with silica gel and preserved at room temperature until use.).The maximum distance between sites, from Chanarian (B-CNR) to San Agustin (M-SAG), was approximately 1600 km.Major surface ocean currents are also shown in Figure 1, following the map by Ravago-Gotanco and Juinio-Meñez [29].A strong westward current, the North Equatorial Current, flows from the Pacific Ocean to the eastern coast of the Philippines, through the study region.The current bifurcates into two currents, and the northward-flowing branch, the Kuroshio Current, is a particularly strong current.The southward current, the Mindanao Current, flows toward the southern Philippines through the eastern coastal area.Other currents flow into the inland sea within the Philippine Islands, with southward currents flowing from the South China Sea to the Sulu Sea.We collected samples from approximately 200 m × 300 m sites in the intertidal or subtidal zones, which were less than approximately 3 m in depth.We maintained a distance of more than 10 m among shoot samples collected within a site to unify the sampling strategy with other studies in this region [13,16,18,20].A young leaf from each shoot was collected, dried with silica gel and preserved at room temperature until use.Luzon and Mindoro.The map was produced using the Generic Mapping Tools (GMT) ver.4.5.5 [30].

DNA Extraction and Microsatellite Analysis
Silica-gel-dried leaves from each shoot were ground using TissueLyser (Qiagen, Hilden, Germany).Genomic DNA was extracted using a modified cetyltrimethylammonium bromide method [31].After extraction, genomic DNA was stored at −30 • C until Polymerase chain reaction (PCR).We determined the microsatellite genotypes of all samples collected in those regions using nine polymorphic loci (Themp030, Themp041, Themp047, Themp048, Themp049, Themp066, Themp103, Themp107 and Themp193) isolated from T. hemprichii and characterised by Matsuki et al. [32] (Table S2).PCR was performed using a Multiplex PCR Kit (Qiagen) in a 6 µL reaction containing 20-200 ng template genomic DNA, 2× Multiplex PCR Master Mix and 0.2 µM (final concentration) primer that is fluorescently labelled with FAM, VIC, NED or PET (Table S2).The PCR settings were 15 min at 95 • C, followed by 30 cycles of 30 s at 94 • C, 90 s at 56 • C and 60 s at 72 • C, with a final extension of 30 min at 60 • C. The fragment length of PCR products was analysed with an automated capillary DNA sequencer (ABI 3130xl Genetic Analyzer; Thermo Fisher Scientific, Waltham, MA, USA) and GeneMapper software ver.4.1 (Thermo Fisher Scientific).

Genotypic Variation and Clonality
The number of multilocus genotypes (N MLG ) among the nine loci was counted using GenAlEx ver.6.503 [33] for each site.Whether two shoots belong to the same genet or to two different genets was determined based on the probability of identical multilocus genotypes being derived from distinct sexual reproduction events by chance (P SEX ) and calculated using GenClone ver.2.0 [34].P SEX was calculated by considering the inbreeding coefficient (F IS ) estimates in the dataset and was used to test for clonal identity and clonal propagation [34,35].To remove duplicate clonemates from further analysis, identical multi-locus genotypes at a given site were considered clonal replicates (P SEX < 0.01).However, we did not evaluate the probability with a low number of multilocus genotypes, due to the low statistical power under such conditions.The index of clonal diversity R proposed by Dorken and Eckert [36] was calculated with the formula R = (G − 1)/(N − 1) for each site, where G represents N MLG accounting for P SEX and N indicates the number of shoots successfully genotyped.To estimate the resolution level of the microsatellites, the combined probability of identity (P ID ) was calculated for each site using GenAlEx.

Genetic Diversity and Recent Bottleneck Effect
The mean number of alleles per locus (N A ) and allelic richness (A R ) across all the loci at each sampling site were estimated using GenAlEx and FSTAT ver.2.9.3.2 [37], respectively.A R is the mean number of alleles per locus and standardised to the site with the smallest number of genets.The values of observed and expected heterozygosity (H O and H E , respectively) and the number of private alleles (P A ) at each site were also calculated with GenAlEx.To estimate the effects of latitude on genetic diversity, the significance of the correlation between genetic diversity (A R and H E ) and latitude was tested using regression analyses.In addition, we evaluated whether populations had experienced a recent bottleneck event using BOTTLENECK ver.1.2.02 [38].We performed Wilcoxon's signed-rank test for significant deviations from population mutation-drift equilibrium using 1000 replications under assumptions of both the infinite allele mutation model (IAM) and the two-phase model (TPM) (30% IAM and 70% stepwise mutation model).

Genetic Differentiation, Structure and Connectivity among Sites
Genetic differentiation within and among populations was estimated using a hierarchical analysis of molecular variance (AMOVA) [39].The genetic differentiation index (pairwise F ST ) between sites was calculated using GenAlEx with 999 permutations to test the significance of each F ST value.The Mantel test was performed using 999 random permutations in GenAlEx for the significance of the correlation between pairwise F ST /(1 − F ST ) values [40] and direct geographic distances between paired sites.Furthermore, we estimated the patterns of genetic differentiation derived from pairwise F ST values among sites via principal coordinates analysis (PCoA) using GenAlEx.Population genetic structure was inferred from microsatellite data using STRUCTURE ver.2.3.4 [41].STRUCTURE analysis employs a Bayesian clustering algorithm to assign genotypes to clusters while minimising the Hardy-Weinberg equilibrium and linkage disequilibrium.Ten replicate runs without prior information were conducted for each K between 1 and 20 using the admixture model and assuming correlated allele frequencies [42].Each run consisted of 1,000,000 Markov chain Monte Carlo replications after a burn-in period of 100,000 iterations.The optimal K was determined using the ∆K method of Evanno et al. [43] as implemented in STRUCTURE HARVESTER [44].Run data were merged using the CLUMPAK server [45].The Q-matrix of each cluster was plotted onto a map in R ver.4.1.2[46] using the POPS script ver.1.2 [47].Furthermore, discriminant analysis of principal components (DAPC) was conducted in R using the package adegenet ver.2.1.8[48], to represent the distribution patterns for each multilocus genotype as plots.This clustering method does not assume a particular model.Therefore, it is free from assumptions regarding Hardy-Weinberg equilibrium and linkage disequilibrium.We allowed 60 principal components to be retained for discriminant analysis.To infer the extent of gene flow and direction, BayesAss ver.3.0.3[49] was used to calculate the migration rate between sites per generation over the last few generations.The following settings were used: number of iterations, 10,000,000; burn-in, 1,000,000; and sampling frequency, 1000.

Scoring of Genotypes and Common Genotypes between/among Sites
Of 1339 shoots sampled, 1261 (94.2%) were successfully genotyped at nine loci.Of these 1261 shoots, 1033 (81.9%) multilocus genotypes were detected at each site and 1058 (83.9%) were identified as genets that were raised from different sexual events based on P SEX values calculated for each site.The combined P ID of each site varied between 9.3 × 10 −6 (P-HND) and 1.5 × 10 −2 (N-TPL).The clonal diversity ranged from 0.07 to 1.00, with a mean of 0.84 ± 0.04 SE (Table 1).For one site (P-PPB) with the lowest R value for clonal diversity, P SEX correction was not applied due to the low resolution of this value with a low number of multilocus genotypes (see Materials and Methods).Five multilocus genotypes were shared among multiple sites: one between W-MOR and D-MCH (147.7 km), one between E-PNT and E-CGB (169.3 km), one between E-PNT and E-SSB (440.4 km) and two between E-CGB and E-BAC (275.6 km).Clonal replicates within each site were removed from further analysis, while replicates between sites were preserved for subsequent analysis.

Genetic Diversity and Bottleneck Event
Allelic richness was standardised based on the minimum number of G, except for site P-PPB.The mean number of alleles (N A ) and allelic richness (A R ) per site ranged from 1.67 (V-TCL) to 5.78 (P-HND) and from 1.62 (V-TCL) to 4.49 (P-HND), respectively.Nineteen private alleles (P A ) were detected from eight sites, and nine of these nineteen alleles were detected in P-HND.The observed and expected heterozygosities (H O and H E ) ranged from 0.172 to 0.496 and from 0.145 to 0.484, respectively.Inbreeding coefficients, F IS , varied between -0.127 and 0.242 (Table 1).The values of N A , A R , P A , H O , H E and F IS per locus for each site are shown in Table S1.Correlation analysis showed that genetic diversity index values had no definite correlation with an increase in latitude (A R : R 2 = 0.044, p = 0.248; H E : R 2 = 0.151, p = 0.028).Some sites appear to have undergone a recent bottleneck, resulting in a significant excess of heterozygosity compared to that expected at mutation drift equilibrium (V-GUI and M-BAT) under both IAM and TPM assumptions (p < 0.05; Table 1).

Genetic Differentiation, Structure and Connectivity among Sites
Hierarchical AMOVA detected considerable genetic differentiation among regions and among sites (Table 2).The proportion of genetic variance was significantly partitioned among sites (F ST = 0.291, p = 0.001).Pairwise F ST values ranged from 0.026 (between W-AND and W-SAN) to 0.744 (between N-TPL and V-TCL) (mean F ST = 0.289 ± 0.005 SE), and all combinations showed significantly genetic differentiation (all Ps = 0.001) (Table S3).These results indicate high levels of genetic differentiation in T. hemprichii among sites.Moreover, analysis of isolation-by-distance revealed a significant correlation between geographic distance and genetic differentiation (Mantel test: p = 0.014) (Figure 2a).PCoA plots were not completely separated by the regions; unexpected complex distribution patterns were found in F ST -based data (Figure 2b).The results of STRUCTURE analysis indicated the existence of a strong population genetic structure (Figure 3 and Figure S1).The ∆K value was greatest at K = 2 (∆K = 2782.92,mean LnP(D) = −20,796.93).When K = 2, Batanes, North Luzon, West Luzon, Mindoro (excluding D-HNR) and Palawan were mainly assigned as the same cluster.Furthermore, V-SII, V-BTY, V-TCL and V-GUI in the Visayas formed a separate cluster at K = 3 (Figure 3).In regions such as Mindoro and the Visayas, the main clusters at each site differed despite the geographic proximity due to the delimited borders of genetic clusters (Figure 3).DAPC also identified continuous genetic clusters, indicating genetic differentiation among sites (Figure 4).In this analysis, V-SII, V-BTY, V-TCL and V-GUI in the Visayas composed isolated clusters.V-GUI was particularly far removed from the clusters of the three sites.More specific clusters were detected visually by removing the clusters of these four populations.In that case, D-HNR was closer to the Visayas compared to the other three populations in Mindoro.Most sites showed a migration rate of less than 0.01 (Table S4).A relatively high migration rate (greater than 0.1) between sites was detected in eight pairs of sites, and the geographic distance ranged from 13.7 km (W-AND <-W-SAN) to 402.2 km (E-BAC <-E-PNT).The highest rate between sites was 0.197, and the distance was 169.3 km (E-CGB <-E-PNT).tected visually by removing the clusters of these four populations.In that case, D-HNR was closer to the Visayas compared to the other three populations in Mindoro.Most sites showed a migration rate of less than 0.01 (Table S4).A relatively high migration rate (greater than 0.1) between sites was detected in eight pairs of sites, and the geographic distance ranged from 13.7 km (W-AND <-W-SAN) to 402.2 km (E-BAC <-E-PNT).The highest rate between sites was 0.197, and the distance was 169.3 km (E-CGB <-E-PNT).S3.S3.

Discussion
We concluded that populations of T. hemprichii have been maintained by both sexual and asexual reproduction, as the index of clonal diversity R ranged from 0.07 to 1.00 with a mean value of 0.84 ± 0.04 SE.The index of clonal diversity R was never low, even in high-latitude sites, where clonal plants could have lower seed production and dominant asexual reproduction [36].The value of clonal diversity was slightly lower than that for E.

Discussion
We concluded that populations of T. hemprichii have been maintained by both sexual and asexual reproduction, as the index of clonal diversity R ranged from 0.07 to 1.00 with a mean value of 0.84 ± 0.04 SE.The index of clonal diversity R was never low, even in highlatitude sites, where clonal plants could have lower seed production and dominant asexual reproduction [36].The value of clonal diversity was slightly lower than that for E. acoroides in the Philippines (mean R = 0.88) [13], which was collected using the same procedure employed in this study.On the other hand, it is higher than values for C. serrulata (mean R = 0.62) [16], C. rotundata (mean R = 0.66) [18] and S. isoetifolium (mean R = 0.38) [20].Large genets in some populations with low clonal diversity (e.g., P-PPB) suggest that seagrass meadows formed as the result of ecological processes over a long time scale.Although meadows of T. hemprichii are primarily maintained by vegetative propagation [23,24], sexual reproduction also contributes to the persistence of populations, even in high-latitude areas.Along with sexual reproduction, asexual reproduction methods, such as clonal expansion of T. hemprichii through elongation of the horizontal rhizome within a meadow or drifting of shoots between sites, may contribute to population maintenance.Three and seven genotypes were detected in E. acoroides and S. isoetifolium, respectively, which were common among sites in this region [13,20].Cymodocea serrulata also exhibited shared multilocus genotypes, with patterns found in the Ryukyu Archipelago [16].In addition, five common multilocus genotypes were found among the sites in this study for T. hemprichii.There was no clonal replication between distant sites in an extremely large-scale geographic analysis of T. hemprichii (see Hernawan et al. [50] and the corresponding microsatellite genotypic data deposited in DRYAD).A monoclonal seagrass meadow of this species does not appear to be present in our focal area.Because the clonal growth of seagrasses is also influenced by environmental conditions such as nutrients, temperature, water salinity, light availability and other physical factors [25,51,52], variation might reflect differences among environments, even in high-latitude areas.
Genetic diversity did not clearly decrease in high-latitude areas, and the central region of this species' habitat distribution has been considered a hotspot of diversity for various organisms including seagrass species [5].This result is similar to the result for E. acoroides in the Philippines [13].This result implies that seagrass populations in the Philippines are close to the centre of the distribution area, and genetic diversity does not distinctly decrease with latitude in relation to the temperature.Meanwhile, genetic diversity appears to decrease in high-latitude areas and populations at the periphery of the distribution ranges of some tropical coastal organisms [53,54].In the case of seagrass species, the Coral Triangle represents a hotspot.The Coral Triangle region is important from the viewpoint of species and genetic diversity.In T. hemprichii, high genetic diversity has been maintained in Indonesia, in the southern part of the Coral Triangle, although it decreases toward the southwest [50].Further collection and analyses that are standardised among sites will provide detailed data on the clonal distribution and genetic diversity among sites in subtropical areas of the Northwest Pacific region.Spatial genetic information is needed to estimate the genetic diversity and allele distribution for local or fine-scale geographic analysis.In addition, the genetic diversity of seagrasses may play a role in the maintenance of species diversity in organisms that live in seagrass beds [55].
The genetic differentiation of T. hemprichii among sites is generally large, even between nearby sites.Pairwise F ST values and migration rates indicated that gene flow is highly restricted and successful recruitment occurred primarily to the natal habitat.Previous studies have found evidence of lower genetic differentiation in the Northeast Pacific reef region for broadcast-spawning invertebrates [56][57][58][59].In broadcast-spawning invertebrate species, mean larval duration ranged from several days to a few weeks and some of the larvae can survive for a few months [60].Meanwhile, seagrass species usually maintain high genetic differentiation among sites, even at local or regional geographic scales [61,62].The pattern of genetic structure appears to have been influenced by the North Equatorial Current (see Figure 1).Such complex genetic structures have been confirmed in other seagrass species in the Philippines.For example, a complex and irregular genetic structure was confirmed in a regional study of E. acoroides in Visayas [63]; populations in central Visayas (V-SII and V-BTY) were assigned to the cluster found in eastern Visayas (V-TCL and V-GUI).Meanwhile, a complex genetic structure at the local scale may be caused by differences in environment or historical fluctuation of populations between sites.As an example of variation in environment, differences in sea bottom (soft bottom and hard bottom) may be influenced in populations of T. hemprichii [64].However, the physiological mechanism of habitat selection has not yet been elucidated.On the other hand, historical fluctuation of populations may influence genetic structure.T. hemprichii populations in Lesser Sunda and the Indian Ocean exhibited irregular genetic structures unrelated to location [50], as did T. testudinum along the Mexican Coast [8].This disagreement may be influenced by historical fluctuations in coastlines and basin boundaries since the Pleistocene.Sea levels during the Pleistocene were ~120 m lower than at present in our focal regions [65].Repeated historical events may have caused complex genetic differentiation and structures.Adaptive reproductive strategies and biological characteristics may also lead to genetic background differences in some areas.
Although most of pairwise migration rates were less than 0.01, the higher migration rate (greater than 0.1) was confirmed in eight pairs of sites which widely ranged from 13.7 km to 402.2 km.Such long-distance dispersal of larvae and repeated stepping-stone dispersal events will contribute to the maintenance of genetic connectivity among populations at ecological time scales (a few generations).Due to ecological and historical migration, the long-distance dispersal of this species may be particularly stochastic, as it has a short maximum floated-fruits-dispersal period of eight days [27].Similarly, the common multilocus genotypes among sites (147.7 to 440.4 km) may be caused by the rare and stochastic dispersal of floated seedlings, which can survive for over 12 weeks in T. hemprichii [28].Seagrass species appear to have the potential for long-distance dispersal of clonal propagules at the regional scale [66].The habitat range of T. testudinum has extended over a long time frame due to stochastic flow [8].The range of T. hemprichii in the Northwest Pacific could have been extended by repeated fragmentation events and historical stepping-stone dispersal of clonal propagules and seeds or fruits from sexual reproduction, as suggested by continuous clusters observed with DAPC.Moreover, clarification of genetic connectivity at the local scale within this region is needed for appropriate population management [63].Proper management depends upon a comprehensive understanding of the relationship between genetic connectivity and the reproductive and biological characteristics of seagrasses species.The specific floating processes of seagrass species need to be investigated further with a focus on hydrodynamic approaches for a better understanding of population connectivity.Successful dispersal rate would need to be analysed using an ocean circulation model due to the extraordinary low probabilities involved.In addition, not only ocean currents but also biotic vectors must be considered simultaneously, because some propagules, seeds or fruits of seagrasses can be transported by manatees, dugongs, turtles, ducks [61,62] or other birds [67,68].Whether propagules, seeds or fruits of T. hemprichii in this region can be carried or swallowed by these vertebrates should be determined.

Conclusions
We conclude that T. hemprichii is genetically separated among sites due to its uncommon and stochastic dispersal in the Philippines.These results suggest the necessity of conservation due to a deficiency in constant migrants from external regions at the contemporary time scale and a difficulty in rapid recovery of the population with high genetic diversity after disturbances.The relationship between environmental factors and reproductive characteristics of seagrasses is also a question that warrants future study.Data from the field of genetics could be used to maximise the effects of conservation efforts for seagrass meadows by providing useful information for retaining the function of carbon storage in Blue Carbon ecosystems [69], contributing to ecological monitoring of the health of seagrasses [1,61].

Figure 1 .
Figure 1.Map of the 33 collection sites in the Philippines.A-C are enlarged views at Batanes, WestLuzon and Mindoro.The map was produced using the Generic Mapping Tools (GMT) ver.4.5.5[30].Major surface currents are drawn following the map by Ravago-Gotanco and Juinio-Meñez[29].

P
-PPB was removed due to low sample size with a large number of clonal replicates (NMLG = 4, R = 0.07).

Figure 2 .
Figure 2. (a) Isolation-by-distance for all sites except for P-PPB based on pairwise FST/(1 − FST) values plotted against geographic distance (km) using the Mantel test.(b) Principal coordinates analysis (PCoA) from the covariance matrix based on pairwise FST values except for P-PPB.The first two axes explained 35.97% of the variation (the first axis explained 24.14%, the second axis, 11.83%).The values of pairwise FST are shown in TableS3.

Figure 2 .
Figure 2. Figure 2. (a) Isolation-by-distance for all sites except for P-PPB based on pairwise F ST /(1 − F ST ) values plotted against geographic distance (km) using the Mantel test.(b) Principal coordinates analysis (PCoA) from the covariance matrix based on pairwise F ST values except for P-PPB.The first two axes explained 35.97% of the variation (the first axis explained 24.14%, the second axis, 11.83%).The values of pairwise F ST are shown in TableS3.

Figure 3 .
Figure 3. Map showing the frequency of Q-matrices from Bayesian clustering analysis implemented in STRUCTURE using clustering without prior information under the admixture model and assuming correlated allele frequencies.K values of these panels are K = 2 to K = 5.

Figure 3 .
Figure 3. Map showing the frequency of Q-matrices from Bayesian clustering analysis implemented in STRUCTURE using clustering without prior information under the admixture model and assuming correlated allele frequencies.K values of these panels are K = 2 to K = 5.

Figure 3 .
Figure 3. Map showing the frequency of Q-matrices from Bayesian clustering analysis implemented in STRUCTURE using clustering without prior information under the admixture model and assuming correlated allele frequencies.K values of these panels are K = 2 to K = 5.

Figure 4 .
Figure 4. Scatter plots of discriminant analysis of principal components (DAPC) results using the R package adegenet.(a) All populations.(b) Twenty-nine populations after removal of four populations (V-SII, V-BTY, V-TCL and V-GUI).

Figure 4 .
Figure 4. Scatter plots of discriminant analysis of principal components (DAPC) results using the R package adegenet.(a) All populations.(b) Twenty-nine populations after removal of four populations (V-SII, V-BTY, V-TCL and V-GUI).

Table 1 .
Population genetic parameters for each site sampled.

Table 1 .
Population genetic parameters for each site sampled.
Initial letter in the code indicates the region: B, Batanes; N, North Luzon; W, West Luzon; D, Mindoro; E, East Luzon; P, Palawan; V, Visayas; and M, Mindanao.N is the number of samples analysed.N MLG is the number of multilocus genotypes.G is the number of multilocus genotypes, including identical genotypes resulting from sexual reproduction by chance, as estimated by P SEX values for each site calculated using GenClone and an index of clonal diversity R = (G − 1)/(N − 1).N A is mean number of alleles for nine microsatellite loci.A R is allelic richness, the mean allele number for nine loci using FSTAT.P A is the number of private alleles for all loci.H O and H E are mean observed and expected heterozygosities for nine loci, respectively.F IS is mean inbreeding coefficients for nine loci.N A , P A , H O , H E and F IS were calculated using GenAlEx.Bottleneck IAM and TPM are probability values of the infinite allele mutation model and two-phase model in the analysis by BOTTLENECK, respectively (bold: p < 0.05).The values of N A , A R , P A , H O , H E and F IS per locus for each site are shown in Table
P-PPB was removed due to low sample size with a large number of clonal replicates (N MLG = 4, R = 0.07).