Composition, uniqueness and connectivity across tropical coastal lagoon habitats in the Red Sea

Background Tropical habitats and their associated environmental characteristics play a critical role in shaping macroinvertebrate communities. Assessing patterns of diversity over space and time and investigating the factors that control and generate those patterns is critical for conservation efforts. However, these factors are still poorly understood in sub-tropical and tropical regions. The present study applied a combination of uni- and multivariate techniques to test whether patterns of biodiversity, composition, and structure of macrobenthic assemblages change across different lagoon habitats (two mangrove sites; two seagrass meadows with varying levels of vegetation cover; and an unvegetated subtidal area) and between seasons and years. Results In total, 4771 invertebrates were identified belonging to 272 operational taxonomic units (OTUs). We observed that macrobenthic lagoon assemblages are diverse, heterogeneous and that the most evident biological pattern was spatial rather than temporal. To investigate whether macrofaunal patterns within the lagoon habitats (mangrove, seagrass, unvegetated area) changed through the time, we analysed each habitat separately. The results showed high seasonal and inter-annual variability in the macrofaunal patterns. However, the seagrass beds that are characterized by variable vegetation cover, through time, showed comparatively higher stability (with the lowest values of inter-annual variability and a high number of resident taxa). These results support the theory that seagrass habitat complexity promotes diversity and density of macrobenthic assemblages. Despite the structural and functional importance of seagrass beds documented in this study, the results also highlighted the small-scale heterogeneity of tropical habitats that may serve as biodiversity repositories. Conclusions Comprehensive approaches at the “seascape” level are required for improved ecosystem management and to maintain connectivity patterns amongst habitats. This is particularly true along the Saudi Arabian coast of the Red Sea, which is currently experiencing rapid coastal development. Also, considering the high temporal variability (seasonal and inter-annual) of tropical shallow-water habitats, monitoring and management plans must include temporal scales.

directly supporting local populations. These ecosystems are naturally stressed on daily to annual-time scales [4][5][6][7][8] and display high environmental variability (e.g. temperature, salinity, primary productivity, nutrients, dissolved oxygen). Such variability is reflected in the biological patterns that alter in response to the new environmental conditions. Lagoon ecosystems are also being increasingly affected by human disturbances that can compromise their ecological and socio-economic values [5,[9][10][11][12].
Subtropical and tropical coastal lagoons encompass a range of essential soft-substrate habitats, such as mangroves, seagrasses and unvegetated bottoms. These habitats are associated with different environmental conditions, resulting not only from their location along the depth profile but also their structural complexity, and biological assemblages [13][14][15][16]. However, while these habitats contain a diverse range of organisms spatial distribution patterns and connectivity in subtropical and tropical lagoon habitats have mainly been assessed using fish and other mobile marine fauna [17][18][19][20][21][22][23]. Studies describing and comparing macrobenthic distribution patterns and the strength of connectivity linkages across different shallow-water tropical lagoon habitats are particularly limited compared to temperate systems (e.g. [15,[24][25][26][27]). Spatial differences in the community can provide information regarding the ecological requirements of species. For example, species able to colonize multiple habitats will most likely be less sensitive to environmental changes, whereas those more directly associated with a specific habitat may be less tolerant to environmental changes. In general, harsher environmental conditions are observed in the intertidal area, dominated by mangrove trees, with conditions being attenuated with increasing depth, a pattern that is associated with a consistent increase in species richness and abundance [28,29]. Indeed, mangrove habitats are characterized as unfavourable environments influenced by high salinity, high fluctuation of temperature, desiccation, and poor soil condition (depleted oxygen) [30]. On the other hand, if undisturbed, seagrass habitats provide comparatively more stable environmental conditions through time [31][32][33] as well as protection from predators [34].
Furthermore, the knowledge about the role of temporal variability in driving macrobenthic patterns is still scarce [35][36][37][38][39][40][41]. While seasonal changes in tropical regions are comparatively less distinct than in temperate regions [42], temporal variability in benthic patterns exists [39,43,44]. Investigating temporal variability patterns is essential to obtain a deeper knowledge of the dynamics and processes regulating lagoon communities. Indeed, considering the current scenario of global climate change, it is critical to better understand how the distribution patterns of organisms in these habitats are changing and particularly how they respond to changes in temperature and other key environmental drivers [1,45]. Temporal variation patterns in the abundance and composition of macrofaunal invertebrates have been intensively studied in temperate coastal ecosystems in relation to environmental variables [46][47][48][49]. Temporal variability in temperature and food availability, for example, can influence recruitment events with consequences for the structure, distribution, and abundance of the community [50][51][52]. Similarly, sediment composition, organic matter, and vegetation cover, which may vary in time, are also main drivers of observed ecological patterns. However, most of those studies have been conducted in temperate regions and, more recently in polar habitats (e.g. [53][54][55][56][57]). Comparatively, less attention has been dedicated to sub-tropical and tropical areas (e.g. [58][59][60][61]). This is even more striking in regards to the assessment of inter-annual variability (but see [40,62,63]).
Assuming that harsher environmental conditions will occur towards the intertidal area (i.e. mangrove habitats), we hypothesise (i) a decrease in species richness (i.e. the total number of species) and in the number of exclusive species from subtidal to intertidal areas, as less resistant species are progressively excluded along the the environmental gradient. We also hypothesise that (ii) shallow water seagrass meadows will harbour higher numbers of species particularly compared with unvegetated bottoms, as a result of habitat complexity, protection from predators and food availability [64][65][66]. Likewise, we hypothesise (iii) that temporal changes will be less evident in subtidal (vegetated and unvegetated) than intertidal habitats [30,67] and that subtidal seagrasses areas will support more stable communities through time. Ecologically related management decisions require a sound knowledge of the biodiversity of the ecosystem. By assessing the variability in spatial and temporal patterns of macro benthic organisms we expand on the existing knowledge on tropical coastal lagoons which are sensitive as well as ecologically and economically valuable.

Macrobenthic community composition: general characterization and connectivity among habitats
A total of 4771 invertebrates were identified within the different habitats surveyed in the lagoon (Fig. 1a), belonging to 272 operational taxonomic units (OTUs) distributed among 11 phyla, 16 classes, 40 orders, and 80 families. Annelida dominated both in abundance and number of taxa, contributing to, respectively, 51.0% and 42.0% of the total values. Sipuncula (15.0%), Arthropoda (13.0%), Mollusca (12.0%), and Echinodermata (7.0%) also contributed to the overall density. Regarding the number of species, Arthropoda (28.0%) and Mollusca (18.0%) were, along with Annelida, the phyla contributing the most to the total number of species.
Both seagrass habitats showed a higher percentage of resident species (i.e. species present in over 85% of the sampling dates in a certain habitat) compared to mangrove and unvegetated areas ( Table 2). In terms of the number of individuals, those taxa contributed to 45.0% and 34.0% for S1 and S2, respectively, of the site's total abundance. S2 showed a more balanced distribution of the four habitat preference traits analysed (i.e. resident, frequent, occasional, rare) and relatively stable numbers throughout the study period ( Table 2). Regardless of the habitat, occasional species accounted for more than 12.6% of the total number of species.
Macrobenthic patterns of variability across the lagoon seascape show that the community was structured by habitat with limited seascape ecological connectivity across the different habitats (Fig. 2a). The environmental data gathered partially explained the multivariate variability of the biological data with the two first axes of the distance-based redundancy analysis (dbRDA) explaining more than half of the constrained variability but only 19.1% of the total variability of the biological communities. The dbRDA plot reinforces a clear separation of the communities inhabiting mangrove areas, S1, and the unvegetated habitat, whereas S2 presented affinities (i.e. higher connectivity) with either S1 or mangrove stations depending on the sampling period (Fig. 2b). Samples from the unvegetated habitat were associated with depth and percentages of medium and fine sand. Seagrass habitats (particularly S1) were separated based on the higher silt and clay (fine particles) content, whereas mangrove habitats presented a slightly higher percentage of coarse sand. Multivariate patterns suggest that the nature of the biotope itself drives the composition and structure of macrobenthic communities. The investigation of temporal variability was undertaken for each habitat separately.

Temporal variability within habitats
The high variability patterns in the seagrass biomass along the study period ( Fig. 3) was reflected in the biological changes but was not fully aligned with the temporal pattern in sea water temperature (Fig. 1b). When analysing the full dataset and regardless the diversity metric considered, S2 consistently presented the highest number of taxa (155, observed; 184.8-219.7, estimated), whereas M2 was the poorest taxa site. Density was also higher at S2 (801.9 ind.m −2 ) and lowest at the unvegetated area (388.8 ind.m −2 ) ( Table 1).
In general, a higher number of OTUs were observed in the subtidal habitats than the intertidal mangrove areas (Fig. 4a), with M2 showing a consistently depressed number of taxa across all sampling dates. Abundance was also generally higher within seagrass meadows (Fig. 4b).
M2 also presented the lowest Shannon-Wiener diversity whereas, in general, higher values were observed at S2 or at the unvegetated habitat (Fig. 4c).
Biological similarity within each habitat was markedly low, ranging from 14% (M2) to 25% (S1) ( Table 2). Both habitats also showed a higher dominance with only four and six species contributing to over 62% of the habitat's abundance, respectively. In the remaining habitats, a minimum of 13 taxa was needed to reach the same level of abundance (Table 2). Except for S1, where none of the dominant taxa was a polychaete, this group dominated all the other habitats. S1 was dominated by a Temporal variation in the structure of macrobenthic assemblages within each habitat examined on the basis of the Bray-Curtis and Jaccard resemblance measures indicated different patterns depending on the habitat in analysis. Major differences were not detected between metrics and therefore only plots for Bray-Curtis matrices are presented (Fig. 5). The results of the Permutational Multivariate Analysis of Variance (PERMANOVA) confirmed different temporal trajectories in the analysed habitats (Table 3). Both resemblance metrics applied to M1 and S1 datasets showed a significant interaction of the main factors (Year x Season). The pair-wise tests indicated for M1 a significant inter-annual difference both in winter and summer. For S1, inter-annual differences were only detected in winter. With regard to seasonal differences, S1 presented significant variability in both years (except in the composition-Jaccard-for 2015) but in M1 differences were only detected in 2014 (Table 3).
Macrobenthic communities at M2 and S2 showed significant inter-annual variability (except for S2 with presence/ absence) (Table 3). Finally, the unvegetated area showed significant and independent seasonal and inter-annual variability (Table 3).

Discussion
This study investigated the distribution patterns of macrobenthic communities inhabiting adjacent shallowwater habitats in a tropical coastal lagoon with particular focus on how they are connected and how communities within each habitat vary over time. Even though ecological seascape connectivity has been previously demonstrated particularly for fish, information on the benthic dynamics in tropical lagoons is still scarce. The Al Qadimah lagoon, likewise other tropical lagoons, encompasses a wide range of habitats including both hard (not addressed here) and soft-substrates. Within the latter, changes in the vegetation cover result in a mosaic of habitats with different sedimentary properties that will determine the structure of local macrobenthic communities [68]. Here, we observed a clear zonation of the benthic communities, driven by habitat-related factors acting  at varying spatial scales [69]. The present results also provided new insights into the temporal variability (seasonal and inter-annual) of different lagoon shallow-water habitats in a tropical seascape.

Uniqueness of lagoon habitats within the seascape
A clear pattern of habitat-dependent association was observed with the different habitats harbouring distinct macrobenthic assemblages. The high spatial variability of macrofaunal patterns is most likely linked to the heterogeneity of the seascape and to the high contribution of rare species to the overall abundance. Recent studies showed that biological variability is driven by the relative high contribution of rare and common species, with rare species playing a major role in the temporal patterns, as a result of their vulnerability to fluctuations in environmental conditions (e.g. [70,71]). Subtidal habitats harboured 70% of the total number of species. Overall, seagrass habitats showed the highest number of taxa, which agrees with previous studies [65,68,[72][73][74]. Variability was, however, high and significant differences within the subtidal area were not detected. The structural complexity provided by the seagrass canopy and the developed rhizome and root systems that contribute to sediment stability may favour the development of diverse communities [70,75,76]. In the tropics, the canopy can play an additional critical role providing shade that can attenuate the effects of sea water temperature [8] that in the study region can reach over 32 °C in the summer. Yet,   [68,81]. Also, despite the increased aeration within the sediment due to the developed root system [82], the decomposition of the high amounts of organic matter Non-metric multidimensional scaling (nMDS) based on Bray-Curtis dissimilarity matrices based on untransformed data, for temporal variation in the structure of macrobenthic assemblages within each habitat. M1 and M2, mangrove; S1 and S2, seagrass; and unvegetated area (Unv) will require increased oxygen consumption and result in anoxic regions and accumulation of toxic products [83,84]. Therefore, vegetated areas with comparatlively lower cover might harbour higher species numbers as a result of species avoiding toxic anoxic conditions in densely covered areas [85].
Within mangrove habitats species encounter harsh physical environmental conditions (e.g. high salinity, hypoxia, desiccation, high concentration of toxins) and in general nitrogen limitation (C/N ratio often > 100; although mangroves in the Red Sea are carbon limited compared to other locations [86]) due to a low nutritional

Table 3 Two-way PERMANOVA model and pair-wise tests based on Bray-Curtis and Jaccard matrices within habitats among seasons and year (Year and Season interaction; Yr x Se)
M1 and M2, mangrove; S1 and S2, seagrass; and unvegetated area (Unv) value of the main source of organic matter, i.e. leaf litter [25]. Under these consitions, populations of a few tolerant/opportunistic species dominate the macrobenthic communities [25,87]. In the present study, the deepest mangrove area (M2) was dominated by only four species, the polychaetes Simplisetia erythraeensis, Ceratocephale sp. and Paucibranchia adenensis, and the bivalve Paphies angusta contributed to over 60% of the total abundance. In the shallowest mangrove area, despite the dominance of polychaetes, the sipunculid (Phascolion (Phascolion) strombus strombus) and some decapods (Diogenes costatus and Thalamita poissonii) were also co-dominant. Decapods are critical players for the ecosystem functioning of these habitats by processing leaf litter and oxygenating sediment through their burrows [88,89] and therefore their dominance in the habitat is not surprising. As observed elsewhere, mangrove habitats showed the lowest number of species compared to nearby seagrass and unvegetated substrates, as previously found [90,91].

Connectedness and stability at the scale of the seascape
In the present study, nearby seagrass meadows differed in cover and depth location, which might have resulted in limited similarity in faunal communities (both habitats shared 35.0% of total number of species). Higher similarities (~ higher seascape connectivity) were detected among subtidal habitats than between those and mangroves (intertidal habitats). Nevertheless, 62 taxa, representing 22.8% of the gamma diversity, were shared between intertidal and subtidal habitats, suggesting that several species may utilize contrasting yet adjacent habitats within the lagoon seascape. Despite the fact that the overlap of species across the five habitats is lower (eight taxa; 2.9% of the total number of taxa) than previously reported [92,93], the present study suggests the connectivity between intertidal and subtidal areas and the need for integrated management measures. The results obtained may result from the low hydrodynamic conditions present but information on the hydrographic patterns is non-existent. The effect of tides can result in displacement of specimens through water movement [94] and depending on their height can also expose organisms to desiccation for variable periods of time, which may hinder the distribution of most of the species toward the intertidal area. Specially, when analysed together, mangrove habitats contributed to 6.6% (M1, 4.4%; M2, 2.2%) of the gamma diversity, contrasting with the unvegetated subtidal area and the seagrass meadows that supported, respectively, 15.1% and 31.3% (S1, 12.9%; S2, 18.4%). Mangrove forests can produce relatively large amounts of organic matter through the conversion of leaf litter into detritus [64], that are later exported to nearby habitats [95][96][97]. Therefore, the proximity of the mangrove stands to shallow water seagrass meadows will most likely contribute to the higher biodiversity and, particularly, higher density observed within seagrasses. The populations of suspension-feeders, such as Barbatia foliata, which was dominant in the seagrass meadow (S1), supports the idea of higher availability of organic suspended particulate matter derived from, among others, nearby mangrove canopies and this higher availability will also support more resident organisms [68,99]. Despite the high temporal variability observed in all habitats, highlighted by the dissimilarity indices, seagrass habitats showed a comparatively higher stability, with the lowest values of inter-annual variability, similar to previous studies in temperate areas [8,98]. These habitats also supported the highest number of resident species (i.e. those present in over 85% of the sampling periods). At the lagoon entrance, the exclusive presence of Schizaster gibberulus, a sea urchin previously associated with the near shore coastal biotope in the region [16], suggests that the unvegetated area may be located along a corridor connecting offshore and lagoon communities, with patterns likely dependent on the hydrodynamic processes [99]. Its position between the lagoon and the open coastal water may also explain the high number of species observed (121), with a large proportion being exclusively associated with this habitat (33.9%). It is worth noting that given the generally low density observed in the Red Sea [16,100], future studies will require to increase the replication across multiple spatial scale to fully understand the dynamics of benthic macroinvertebrates under low nutrient, high temperature, and high salinity conditions. Therefore, conclusions related to abundance and diversity should be interpreted with caution.
The present findings reinforce the need for an integrated understanding of shallow-water habitats from a seascape perspective, in opposition to a fragmented analysis of the isolated habitats [21,101,102]. Whereas the latter may be relevant when looking at particular species, the contribution of each habitat to the dynamics of the whole macrobenthic assemblages is relevant and should not be disregarded by managers when aiming for marine biodiversity conservation. Indeed, in tropical regions, seagrass beds and mangroves have been reported as key nursery areas for several reef fishes such as parrotfishes (Labridae, Scarini), grunts (Haemulidae) and snappers (Lutjanidae) [103][104][105][106] that rely on the macrobenthos as food resources. Large-scale migrations (over 30 km) by juvenile snappers, between inshore nursery habitats and reefs in the central Red Sea have been reported [22]. Also, mangrove forests have been linked to enhanced biomass and biodiversity of coral reef fishes [18,21,104,107,108]. Sustained connectivity of the habitats may enhance the resilience of coral populations to recover after disturbance [107]. Therefore, disturbing the corridors connecting coral reefs with other inshore habitats may even have consequences for reef conservation at a local scale.

Conclusion
Overall, the present study confirmed a decreasing gradient in the total number of species and number of exclusive species towards the mangrove habitats. It also supports the role of seagrass habitat complexity in promoting diversity and density of organisms. Nevertheless, high and stable seagrass cover does not necessarily result in the highest biodiversity levels. But the presence of these plants plays an essential role in the biodiversity of coastal lagoons. Seagrass habitats in contrast to mangrove forests and the unvegetated area show lower interannual variability and higher number of resident species, suggesting more stable communities.
Current findings highlight habitat-structured patterns and persistent patchiness evidenced by a limited number of overlapping species (dominance of habitat specialists over generalists) within the seascape. This is particularly relevant considering the proximity of the analysed habitats but may result from the low dominance levels compared to temperate regions [92,98,109]. Nevertheless, 22.8% of the gamma diversity was represented by taxa spanning between subtidal and intertidal habitats. Hence, holistic, i.e. interconnected seascape management approaches, rather than those focusing on single habitats should be prioritized to protect biodiversity and fisheries [22,110,111].

Study area and sampling design
The present study was carried out in the Al Qadimah lagoon (22° 22′ 39.3″ N, 39° 07′ 47.2″ E) located in the central region of the Saudi Arabian Red Sea (Fig. 1a). This shallow lagoon (average depth 2.19 m) has an approximate area of 14 km 2 and is not impacted by direct anthropogenic disturbances typical of other coastal lagoons (e.g. freshwater or sewage discharges, fisheries, habitat destruction from coastal development). It is, however, situated between two urbanized areas, which are increasing in size (King Abdullah University of Science and Technology, 7000 inhabitants; King Abdullah Economic City, currently 5000 inhabitants but it is expected to reach 50,000 in the near future) but that are not directly connected with the lagoon. Hence, it offers a rare opportunity to study the natural roles of environmental drivers in shaping macrobenthic communities inhabiting such critical wetlands.
Scattered along the extent of its margins, well-developed mangrove stands of Avicennia marina are observed.
The bottom of the lagoon, particularly in the inner areas is characterized by more or less fragmented seagrass meadows. To depths of approximately 50 cm, Cymodocea rotundata is the dominant species with smaller patches of Cymodocea serrulata also being present. Below this depth, seagrass meadows are mainly characterized by mono-specific stands of Enhalus acoroides down to 2 m depth. Towards the sea, unvegetated bottoms with either sponges mixed with coral rubble or sand progressively replace seagrass meadows.
In the Red Sea, there are two marked seasons (Fig. 1b), winter (November-April) and summer (May-October). In order to investigate inter-annual and seasonal changes in macrobenthic patterns, samples were collected in two different periods in winter (January; March) and summer (June; September) of 2014 and 2015. Five permanent soft-sediment habitats typical of tropical coastal lagoons were selected: 1. upper mangrove area (M1); 2. deeper mangrove area (M2); 3. shallow seagrass meadow (S1, mix meadows of Cymodocea serrulata interspaced with Cymodocea rotundata; relatively high cover all year round); 4. deeper seagrass meadow (S2, monospecific stands of Enhalus acoroides with high variability in the vegetation cover throughout the study period); and 5. unvegetated soft-sediments (Fig. 1a). The unvegetated sandy substrate was located between 8 and 10 m depth. Due to the widespread distribution of seagrasses, mangroves and in order to minimize the direct influence of those habitats on the colonization patterns of unvegetated areas, the site was located at the entrance of the lagoon.

Sampling strategy
At each habitat and sampling period, conductivity, temperature, and depth (CTD) casts were carried out with a multiparameter probe (OCEAN SEVEN 316 Plus and 305 Plus). The CTD casts also recorded oxygen saturation in the water column. Water samples for the analysis of chlorophyll a (chl a) were collected using a Niskin bottle at each station (2 L per station). Sediment samples were collected using a 0.1 m 2 Van Veen grab in the seagrass meadows and the unvegetated area (subtidal stations), whereas in the mangrove habitats (intertidal), samples were collected using hand corers (3 × 10 cm i.d. making one replicate; total area per replicate ~ 0.024 m 2 ). In 2014, two replicates at each site and sampling date were taken for the study of macrobenthic communities, with additional samples being collected for the study of environmental variables (grain particle size distributions and organic matter content). In 2015, the same approach was followed increasing the number of replicates for the study of macrobenthic communities to three. Macrobenthic samples were sieved through 1 mm mesh screens and preserved in 96% ethanol.

Laboratory analyses
In order to estimate the primary production in the sampling area, the concentration of chl a was quantified by fluorescence using the EPA method 445.0 [112]. Water samples were filtered using GF/F filters as soon as we arrived at the laboratory. The filters were then preserved at -80 °C until extraction of the pigments. 10 ml of 90% acetone were used for each extract and left for 24 h in cold and dark conditions to minimize degradation. The procedure was undertaken in low light conditions to minimize degradation. A Turner Trilogy ® fluorometer (Turner Designs) was used to quantify the chl a content using an acidic module. The degradation of the chlorophyll a to phaeophytin was accomplished by acidifying the sample with 60 µl of 0.1 N HCl.
Sediment samples were sorted after all the vegetation associated with sediment was removed. Organisms were whenever possible identified to the species level. Vegetation biomass (seagrass leaves, roots, and mangrove material) was quantified per replicate.
Grain particle-size distribution was quantified after initial wet sieving of the samples (63 μm mesh) to separate the silt and clay fraction from sandy fractions and gravel. The retained fractions were dried at 80 °C for 24-48 h. The dried sandy and gravel sample was then mechanically sieved by using a column of sieves to separate the sandy fractions and the gravel as follows: < 63 μm, siltclay; 63-125 μm, fine sand; 250-500 μm, medium sand; 1000-2000 μm, coarse sand; > 2000 μm, gravel.
The organic content of the sediments was determined by loss on ignition (LOI). Sediments were dried for 24-48 h at 60 °C and then the samples were placed in the muffle furnace at 450 °C for 4 h. After cooling in a desiccator for 30 min, samples were weighed and the LOI was calculated using the following equation [113]: where: LOI = Organic Matter content (%), W i = Initial weight of the dried sediment subsample; W f = Final weight after ignition.

Data analysis General patterns
Macrobenthic patterns were analysed through a combination of univariate and multivariate techniques. Several univariate metrics were calculated including the total number of taxa (S, species richness), density (ind. m −2 ), and Shannon-Wiener (H′). Considering the different sampling methods, and the dependency of species richness on sample size [114], estimates of species diversity were also calculated and compared with S. The nonparametric species richness estimators used: Chao 1, Jacknife 1 order and Bootstrap all follow an asymptotic approach to estimate the number of undetected species richness. These estimators are commonly used in ecological studies because they are simple, intuitive, relatively easy to use and perform reasonably well [115]. The biased corrected form of Chao 1 estimator [114,116] uses the number of singletons and doubletons to estimate the lower bound of species richness. The first order Jacknife estimator [117] assumes that the number of species that are missed equals the ones that were seen once (singletons). The Bootstrap estimator is based on the assumption that if the same data is resampled with replacement the number of missing species after resampling will be similar to those missed originally [117]. All estimators were calculated using the open source software R [118] using function "specpool" from "vegan" package [119]. Abundance data was used for the calculations of all estimators. In order to have a balanced number of replicates, the analyses were conducted for two replicates, with those collected in 2015 being randomly selected. Preliminary analysis showed that the same general patterns in composition and alpha-diversity were obtained for 2014 and 2015 datasets.
To visualize multivariate patterns of abundance in macrobenthic communities within the seascape, non-metric multidimensional scaling (nMDS) was applied based on the Bray-Curtis dissimilarities. Given the differences among habitats for some dominant species, when comparing habitats (i.e. full dataset), Bray-Curtis dissimilarities matrices were calculated using untransformed abundance data. Separate nMDS plots were generated for each one of the sites for a better visualization of the temporal variability. These analyses were also based on untransformed data. Within each site, significant variability in the multivariate patterns over time was analysed initially according to a three-factor design (Year; Season; Date, nested within Season) using Permutational Multivariate Analysis of Variance (PERMANOVA). As the factor "Date" was found not significant, and to increase the power of the analysis, a two-factor PERMANOVA was applied. Whenever significant differences in the interaction term were detected (i.e. Year × Season), pair-wise tests were conducted.

Connectedness within the seascape and stability patterns over time
A preliminary investigation of the patterns of variability across the seascape was carried out to identify generalist versus specialist taxa, i.e. those that span across multiple habitats versus those that are particularly associated with a specific habitat, respectively. We aimed to characterize the main differences in the community patterns in terms of shared and exclusive species that could determine the cause of the connectivity across the lagoon. This analysis was conducted based on the whole dataset, disregarding the seasonal and annual changes, as our main question was related to the constancy of spatial changes in different habitats.
Finally, we analysed the frequency of occurrence of species in each habitat during the study period. Species were classified based on Habitat Preference Trait as follows: (i) resident, present in over 85% of the sampling dates (i.e. eight events); (ii) frequent, observed between 50% and 85% of the dates; (iii) occasional, presence registered in between 25% and 50% of the sampling occasions; (iv) rare, observed in less than or equal to 25% of the sampling dates; (v) seasonal, only observed in one season but in both years.
Community stability was also examined over the sampling period within each habitat based on the indices Bray-Curtis (community structure) and Jaccard (presence/absence; composition). Within each habitat, variability between all pairwise comparisons among terms of interest (e.g. within and between seasons; within and between years) was analysed. We established that low levels of similarity are related to high variability in the macrobenthic communities over time, whereas high similarity is indicative of more stable communities.

Relationships between environmental variables and assemblage structure
Distance-based redundancy analysis (dbRDA) was used to assess the relationship between each environmental variable and the variation in the community structure (given by the direction and length of vectors for each variable). The variables used for the analysis were salinity, temperature, depth, grain size fractions, organic matter content (% LOI), and chl a. Three of the variables were transformed to reduce skewness, namely the fines and coarse sand fractions of the sediment (square root) and organic matter content (natural log). Marginal tests are used to show the significance of each variable individually to the model and sequential tests show the best subset of explanatory variables that explain the biological patterns.