The study area was located in the provinces Khanh Hoa (Thuy Trieu and Van Phong Bay) and Phu Yen (Xuan Day and Cu Mong) (Figure 1A). Studies determined seagrass distribution over time and identified possible drivers for seagrass loss in Khanh Hoa which revealed a total loss of 74.2% of seagrass area between 2008 and 2018 and various human activities in the coastal zone (e.g., tourism, shellfish collection and aquaculture) appeared as main causes of seagrass decline (Chen et al., 2016; Khanh Ni et al., 2020). In Van Phong Bay, a total decrease of 35.8% over 30 years (1988–2019) with an approximate loss of 46.5–87.3 Ha per decade has been shown (Vo et al., 2020). Consequently, seagrass meadows in Thuy Thrieu and Van Phong Bay are fragmented and restricted to shallow waters with a maximum depth of 3.5–5 m near the coast (Khanh Ni et al., 2020; Vo et al., 2020). Thieu Thrieu is characterized by high turbidity with spatio-temporal variations driven by rainfall and resuspension of the sediment with maximum values for NTU exceeding 13.1 (Quang et al., 2017). Additionally, studies confirmed enrichments of Cu and other metals (As, Cr, Ni, Pb, Zn) due to shipping activities in Khanh Hoa ranging from moderate to severe (Dung et al., 2014) which could affect seagrass beds (Nguyen et al., 2017). The rapid expansion of aquaculture area (especially lobster farming) in recent years has caused eutrophication along the coast of Phu Yen (Hoang et al., 2018), however, no information on seagrass distribution and change is currently available for this province.

A total of 384 individual shoots of E. acoroides were collected in April 2019 at 8 populations in four different lagoons along the South Central Coast of Vietnam (Figure 1A and Table 1). We sampled along transects with an interval of 1 meter and collected approximately 50 individual shoots for each site. This sampling method ensures an adequate resolution to determine clonal characteristics and fine-scale structure within populations and should provide sufficient unique multilocus genotypes (MLG’s). Samples were dried in envelopes and brought to Vrije Universiteit Brussel (Belgium) for further analysis. GPS coordinates at the starting point of each transect were taken. Sampling depth varied between 50 and 140 cm. Turbidity, salinity, pH and temperature were measured (Horiba U-50 Series multi-parameter) in all sites and ranged from 2.5 NTU to 17.5 NTU, 34.1 ppt to 37.2 ppt, 8.62 to 9.02, and 31.3 to 33.75°C. The surface layer of sediment consisted of fine sand and mud. High abundances of algae and Cnidaria (jellyfish and anemones) were present in the water and on seagrass shoots. All sampled meadows were monospecific except for VP3 where the seagrass species Thalassia hemprichii was found intermingled with E. acoroides and TT1 where a monospecific patch of Halophila ovalis existed at 2 m distance from the sampled transect.

Genomic DNA was extracted from approximately 20 mg of dried leaf tissue using the E.Z.N.A. SP plant DNA Mini kit (Omega bio-tek, Norcross, GA, United States). Eleven polymorphic microsatellite markers (EA447, EA1461, Eaco_001, Eaco_002, Eaco_009, Eaco_048, Eaco_050, Eaco_051, Eaco_052, Eaco_054, Eaco_055) developed for E. acoroides (Gao et al., 2012; Nakajima et al., 2012) were selected for genotyping (Supplementary Table 1). Primers were fluorescence-labeled with 4 different dye-labels (6FAM/VIC/NED/PET) and a primer mix was made by mixing 0.2 μM of each primer. Multiplex PCR reactions consisted of 6.25 μl master mix (Qiagen Multiplex PCR kit), 1.25 μl primer mix, 2.5 μl H₂O and 2.5 μl of genomic DNA. PCR was performed in a thermal cycler (Bio-Rad MyCycler) with the following conditions: an initial denaturation of 95°C for 15 min followed by 35 cycles of: 30 s denaturation at 95°C, 90 s annealing at 57°C and 80 s elongation at 72°C followed by a final extension of 30 min at 60°C. PCR products were separated on an ABI3730XL sequencer (Macrogen, Seoul, Korea) and allele sizes were determined with GeneMarker (v.2.60) (Holland and Parson, 2011).

We suggest that the area of human-induced land use is an indirect but integrative proxy for human-induced disturbance especially in habitats (i.e., lagoons) characterized by high environmental fluctuations (water depth, turbidity, salinity). Therefore, we quantified land use in lagoons using remote sensing data. Satellite imagery was achieved by Sentinel-2 and downloaded from ESA Copernicus Open Access Hub¹. We targeted Sentinel-2 L2A images with atmospheric correction sensed between February and June 2019 with a maximum cloud cover of 5%. Twelve spectral bands were resampled to a pixel size of 10m, merged and exported as raster files using Sentinel Toolbox SNAP (v.2.0.2) (ESA Sentinel Application Platform). Satellite raster files were merged with the software QGIS (v.3.10.3) (QGIS org, 2021) and consequently, we obtained one raster file with spectral data covering all four lagoons. Based on the satellite-2 image we identified training samples and classified six land use classes (forest, urbanization, agriculture, bare soil, unvegetated area, water) with a random forest classifier using EnMap-Box in QGIS (Supplementary Figure 3). We carried out an accuracy assessment analysis based on validation samples including total accuracy (90.48%), kappa accuracy (84.74%) and user’s and producer’s accuracy for each class accompanied with the confusion matrix (Supplementary Tables 4, 5). Additionally, we classified aquaculture manually as an additional land use class. Catchment area boundaries of the four sampled lagoons were determined using Digital Elevation Models in QGIS. For each lagoon, we calculated the catchment area excluding water and the absolute and relative area of land use for each class within the catchment. To assess the effect of land use at different spatial scales we calculated catchment and land use area around populations within spatial buffers of four different radiuses (0.5, 1, 5, and 10 km) and within the catchment area of the total lagoon (Supplementary Table 6).

We performed a Principal Component Analysis (PCA) on genetic parameters (AR, H_e, β and Sp_gen) using R (R Core Team, 2021). Populations were organized along the first axis (62,39%) with an increase in Pareto index (β), a decrease in genetic diversity (AR, He) and a decrease in fine-scale genetic structure (Sp_gen) with increasing PC1 values (Figure 2A) and transformed the genetic variables based on the first axis of the genetic PCA (PC1.GENETIC). Based on the acquired land use data, we selected 5 factors that were highly consistent with human-induced land use (catchment, urban, agriculture, aquaculture and forest area) and tested their effect and the effect of spatial scale on the transformed genetic variables. Because the area of urbanization, agriculture and forest was highly correlated with the catchment area we used the relative areas of these land use factors in the analyses to decouple the effect of catchment area from land use. We determined the effect for each land use class separately by linear regressions over different spatial scales and corrected for multiple tests within each scale (Holm-Bonferroni method). The strongest relationship between human-induced land use and the transformed genetic variables was found on the 10 km spatial scale; consequently, we performed an additional PCA on land use variables on this scale. Two groups of populations along the first axis (81.24%) were identified with TT1, VP1 and VP2 characterized by large catchment sizes, large areas of aquaculture, agriculture and urbanization, and low areas of forest (Figure 2B). We transformed land use variables based on the first axis (PC1.LANDUSE) and performed an additional linear regression between the transformed variables of both PCA’s (PC1.GENETIC, PC1.LANDUSE) to test the overall effect of land use factors within a radius of 10km on clonal diversity, genetic diversity and fine-scale genetic structure.

A total of 247 unique multilocus genotypes were identified for 367 E. acoroides samples. Clonal diversity differed between populations with large variation in clonal richness (R = 0.67) ranging from 0.36 to 0.89 and Pareto index (β = 1.26) from 0.47 to 2.26 (Table 2). Clonal evenness (V = 0.86) was overall high and similar between populations indicating that lineages of multilocus genotypes are present in more or less equal abundances (Table 2). The maximum number of ramets per genet (N_Gmax = 4.5) and clonal subrange (CS = 5.9 m) varied with a range from 2 to 8 and 1 to 13 m, respectively (Table 2). Remarkably, XD2 which has the second-largest clonal subrange only has a small number of maximum ramets per genet which is explained by one clone harboring a space up to 11 m but which only consists of four ramets and could be an indication of small scale dispersal by clonal fragments or fragmentation of ramets. β was strongly correlated with R (R² = 0.94), negatively correlated with N_Gmax (R² = –0.73) but not significantly correlated with V (R² = –0.65) showing that populations with a low Pareto index have higher proportions of clones and contain larger genets. For all populations, the aggregation coefficient (A_C) was significant and ranged from 0.47 to 0.77 which infers high spatial aggregation of ramets (Table 2).

The among ramet spatial autocorrelation analysis gave a significant mean kinship coefficient for the first distance class (F_i1ram = 0.190) for all populations with the logarithmic regression slope (b) being significant over 50 m except for CM2 (Supplementary Table 2). The among-ramet Sp-statistic (Sp_ram) varied greatly between populations from 0.007 to 0.211 (Table 2). Similar results were obtained for the among genets analysis with lower but still significant kinship values (F_i1ram = 0.091) for all populations except VP3, significant slopes for TT1,VP1,VP2 and XD1 ((Supplementary Table 2) and among-genet Sp-statistics (Sp_gen = 0.041) ranging from 0 to 0.098 (Table 2). Sp_gen was negatively correlated with R (R = –0.78) and positively with N_Gmax (R = 0.89) indicating that clonal structure affects the fine-scale structure of populations even when clone pairs are not considered. TT1, VP1, VP2, and XD2 have a higher fine-scale genetic structure for both analyses with a maximum distance ranging from 5 to 10 m depending on site and analysis (Supplementary Figure 1).

Number of alleles (N_a = 3.05), effective number of alleles, N_e = 1.88) and allelic richness (AR = 2.79) varied between populations and ranged from 2.09 to 3.64, 1.40 to 2.35 and 1.95 to 3.48, respectively (Table 2). Values for expected and observed heterozygosity were similar within each population and the within-population inbreeding (F_is = −0.030) was close to zero and nowhere significant (Table 2). AR and H_e were strongly correlated (R² = 0.94) indicating that those alleles that contribute to a higher allelic richness in populations are generally found in a heterozygous state.

The three-level AMOVA analysis (Table 3) showed that a high proportion of the molecular variance was explained among regions (33%) but only a small proportion among populations within lagoons (5%) with F_rt and F_sr values being, respectively, 0.327 and 0.081. These results, together with high Fst P-values of populations between and low Fst P-values within lagoons (Supplementary Table 3), indicate high differentiation between populations among lagoons but low differentiation and thus high genetic connectivity of populations within lagoons. The highest Fst P-value (0.518) was found between XD1 and CM1 which could indicate a barrier between neighboring lagoons XD and CM (Supplementary Table 3). The AMOVA inbreeding coefficient (F_is = 0.044) was low with a small proportion of the variation explained by individuals within populations (3%). In contrast, the largest proportion was found within individuals (59%) (Table 3). Bayesian clustering analysis revealed k = 3 (Delta K = 1371) as the most likely number of gene pools (Supplementary Figure 2) with individuals clustered almost exclusively to their respective lagoon (blue = VP lagoon, orange = XD lagoon, purple = CM lagoon) (Figure 1B) except for TT1 where most individuals were fully assigned to the VP cluster and a small proportion was admixed between the VP and XD cluster. With k = 4, the largest proportion (92%) of individuals of TT1 were grouped together in an additional cluster (green) (Figure 1C) and verifies the AMOVA results which indicate high connectivity within but low connectivity between lagoons.

Variation in land use variables was found with a consistent pattern on large spatial scales (10 km, total lagoon) with TT1, VP1, VP2 present in areas with larger catchment sizes and higher levels of disturbance characterized by larger areas of urbanization, agriculture and aquaculture compared to VP3, XD1, XD2, CM1, and CM2 (Table 1 and Supplementary Table 6). We identified no similar and uniform pattern at smaller scales (0.5, 1, 5 km) where differences in land use area between populations were scale-dependent (Supplementary Table 6). We found no significant effect of human-induced land use area on the measured environmental parameters (turbidity, salinity, pH, and temperature). Linear regression analyses per land use factor on PC1.GENETIC showed significant negative effects of aquaculture, agriculture and urbanization on large scales (10 km, lagoon) and of catchment size for the total lagoon (Table 4). No significant relationship was found on smaller spatial scales (0.5, 1, 5 km) (Table 4) showing that human-induced land use only has an effect on clonal diversity, genetic diversity and fine-scale genetic structure of populations on large spatial scales. A positive significant relationship between PC1.LANDUSE and PC1.GENETIC (p = 0.002) revealed a decrease in Pareto index (β), an increase in genetic diversity (AR, H_e) and an increase in fine-scale genetic structure (Sp_gen) in populations inhabiting areas with larger human-induced land use areas, larger catchments areas and smaller forest areas on the 10 km spatial scale (Figure 2C).

Human activities consisting of agriculture, urbanization and aquaculture were broadly present in the sampling area though the area of human-induced land use differed between lagoons and over spatial scales. Water transparency was low at sites with high abundances of epiphytic algae, anemones and jellyfish in the water and on seagrass shoots (Figure 1D) which indicate that the sampled lagoons were highly eutrophic (Purcell et al., 2007); however, no relationship between land use and turbidity and no effect of turbidity on the genetic variables of E. acoroides populations could be found. Nevertheless, extreme reductions in seagrass coverage in our study area have been shown and are most likely associated with extensive human activities in lagoons and high levels of turbidity (Khanh Ni et al., 2020; Vo et al., 2020). No evidence of inbreeding or loss of genetic diversity was found for E. acoroides and demonstrates the high level of resilience of E. acoroides populations toward environmental perturbations. Still, it remains challenging to infer demographic responses directly from genetic diversity because no knowledge is available on the genetic diversity of populations before, during and after disturbances (Arnaud-Haond et al., 2010). The level of genetic diversity in our study is comparable to populations of Hainan Island, China (Yu et al., 2018, 2019), using a nearly similar set of loci, though lower compared to populations in the Philippines based on comparisons among analogous loci (Putra et al., 2018). Populations in Central Vietnam and South China are located on the margins of the species distribution range of E. acoroides (Short and Waycott, 2010) and lower levels of genetic diversity are most likely caused by range edge effects (Arnaud-Haond et al., 2006; Nakajima et al., 2014). A similar study on E. acoroides along the South Central Coast of Vietnam revealed higher levels of genetic diversity though nearly all sites were significantly inbred (Nguyen and Jutta, 2019). This presumably can be explained from a lower number (five) though highly polymorphic microsatellite markers used and their lower sample size (N = 10) which usually does not allow to detect a sufficient number of heterozygous genotypic combinations.

Both AMOVA and STRUCTURE analyses show that lagoons are strongly differentiated and confirm the low levels of genetic connectivity found between populations in lagoons with the open sea (Nguyen and Jutta, 2019). Lagoons are isolated by land barriers with often small inlets and thus may act as barriers to gene flow. In contrast, low differentiation and the absence of genetic structure among populations within lagoons indicate that local recruitment by seedlings, originated within the meadow or within sites in the same lagoon, contributes greatly to the maintenance of E. acoroides populations. Olesen et al. (2004) studied the recovery of available spaces after disturbance in seagrass meadows and showed that sexual reproduction is highly important for the colonization of empty gaps. Direct mortality due to disturbance and the subsequent recruitment of new seedlings in available spaces from various sources within the lagoon could therefore positively affect genetic diversity which may explain the increased levels of allelic richness and heterozygosity under high levels of disturbance in our study.

A higher relatedness of individuals over small distances within meadows (5–10 m) was found in more disturbed sites on both ramet level, due to the aggregation of clones over several meters, and genet level. The aggregated distribution of pollen in male floating flowers and of seeds that are released simultaneously from the floating fruit could explain the fine-scale genetic structure found for E. acoroides among genets rather than pollen and seed dispersal restriction (Yu et al., 2019). The effect of the aggregated dispersal of propagules on the fine-scale structure may be elevated in disturbed sites when population density is low as a result of high mortality and decreased reproductive output. The successful settling and germination of seeds from the same fruit which are resistant to environmental stressors may be an additional explanation for the genetic structure within disturbed meadows. The influence of clonal structure on the fine-scale genetic structure among genets is less clear and understudied. Possible mechanisms are successful mating between a few clonal lineages, flower synchronization among ramets (Jahnke et al., 2015) or the aggregation of clones with slightly different but related genotypes as a result of somatic mutations (Yu et al., 2020).

No shared multilocus genotypes were found among populations, even within the same lagoon, hence genetic connectivity between populations of E. acoroides is only realized by sexual propagules and confirms the assumption that successful dispersal of vegetative fragments of large-sized seagrass species is low to non-existent (Lai et al., 2018). The genetic structure among lagoons is not explained by geographic distance with distant lagoons Thieu Thrieu and Van Phong Bay (~70 km distance) being less differentiated than neighboring lagoons Xuan Day Bay and Cu Mong (~14 km distance). Populations among lagoons with high levels of disturbance (TT and VP) appear to be less differentiated in our study. Human mediated dispersal of fruits and seeds between areas with higher levels of human activity may increase genetic connectivity among disturbed sites (TT and VP) (Orth et al., 2006; Bullock et al., 2018) though no evidence of direct human-mediated transport of propagules is available for E. acoroides. Clonal growth contributes to the persistence of populations under disturbance.

Populations of E. acoroides in our study are sustained by clonal growth and local seed recruitment which confirms the mixed reproduction strategy in populations of Hainan Island affected by physical disturbance and catchment effluents (Yu et al., 2019). Becheler et al. (2014) studied temporal dynamics of clonal structure for Zostera marina and identified similar patterns of clonal richness and structure whereas populations persisted by both initial recruitment, with a stable core of clones which dominance increased over time, and repeated (annual) recruitment of seedlings. Our results show a high proportion of clones and clonal size in sites with large areas of human-induced land use. Arnaud-Haond et al. (2010) found a correlation between mortality and clonal richness in Posidonia meadows with resilient sites having a higher number of large-sized clones due to successful competition under disturbance. E. acoroides has a very low rhizome growth rate (ca. 3 cm/year) (Marbà and Duarte, 1998) which suggests that large genets for this species may have persisted and weathered environmental change over hundreds of years. The selection of resistant genets under strong disturbance and the subsequent reproduction by clonal spread of these persistent genets may thus explain the positive effect of human-induced land use on clonality. Clonal growth and structure also have direct fitness assets including competitive advantages for space, higher reproductive success when multiple ramets have the potential to mate (Pan and Price, 2001) and the translocation of resources between ramets (Marbà et al., 2002; Jahnke et al., 2015) though no evidence for carbon allocation in turbid waters has been found for E. acoroides (Kiswara et al., 2005). Recent studies on Z. marina show that somatic mutations and epigenetic variation originated from asexual reproduction can create phenotypic variation between ramets on which selection can act and may result in the adaptation of genets toward environmental change (Jueterbock et al., 2020; Yu et al., 2020). Decreased levels of successful sexual reproduction may be an additional explanation for the negative effect of disturbance on clonal richness. Light availability has been identified as an important driver for flowering in E. acoroides with a lower proportion of female flowers and fruits in deeper and turbid sites due to decreased levels of photosynthesis (Rollón et al., 2003). Fragmentation of populations and reductions in shoot density which decreases pollen trapping potential of leaves and flowers can affect pollination and thus also fruit set (Vermaat et al., 2004). A negative association between clonal and sexual reproduction due to resource availability and allocation in individuals may exist but no evidence for this trade-off is currently present for E. acoroides (Jahnke et al., 2015).

Our findings highlight the importance of both sexual and asexual reproduction for the persistence of populations of E. acoroides in areas with large-scale human-induced land use. Meadows are sustained by clonal growth and local recruitment of seedlings with lagoons acting as barriers to gene flow. Large persistent genets contribute to the resilience of meadows of E. acoroides due to their resistance toward environmental pressures. However, recruitment by seedlings remains crucial for the colonization and recovery of distant populations, the maintenance of genetic diversity and the adaptation toward environmental change and should be prioritized in management though restricted connectivity among lagoons and local adaptation may complicate conservation measures. To fully understand the response and resilience of populations of E. acoroides toward disturbance, spatial analyses combined with temporal dynamics of clonal richness, structure and genetic diversity should be investigated in future research. Moreover, it is important to identify the direct local environmental drivers caused by human-induced land use and their relative effect on populations of E. acoroides. The studied meadows are most likely remnants of large, connected multi-species seagrass beds which disappeared with increasing human activities (Chen et al., 2016; Khanh Ni et al., 2020) and E. acoroides being the only remaining seagrass species persisting under high levels of human-induced disturbance. We hypothesize that seagrass ecosystems with E. acoroides as key species will eventually shift to systems dominated by macro-algae, epiphytic algae and phytoplankton when environmental conditions continue to diminish, combined with the effects of climate change and potential large disturbance events in the nearby future. Crucial steps for protection and recovery is to improve water quality by managing human-induced land use and related activities and maximizing ecosystem resilience taking into account all key actors (e.g., biological features of seagrass species and communities, functional important species trophic interactions, connectivity with other ecosystems) (Unsworth et al., 2015). Only then will it be effective to invest in practical management measures (e.g., planting seeds or the translocation of plant fragments) for the conservation of seagrass ecosystems.