The Codevelopment of Mangroves and Infaunal Community Diversity in Response to the Natural Dynamics of Mud Deposition in French Guiana

: The sustainability of mangrove ecosystems requires a knowledge of their spatiotemporal variability as a function of regional properties. The unique coastal ecosystems of the mangrove belt along the coast of the Guianas in South America are inﬂuenced by cycles of a massive accretion of mud supplied by the Amazon River and wave induced erosion. This study characterized, for the ﬁrst time, how benthic infaunal assemblages, as proxies of mechanisms of mangrove resilience, were structured by the natural growth track of Avicennia germinans dominated mangroves in French Guiana. We sampled 4 mobile mud stations and 27 consolidated mud stations distributed over 9 tidal transects from bare to vegetated mudﬂats colonized by young mangroves during the dry season. We collected a complete dataset of sediment and vegetation variables together with the benthic meso-(>0.25 mm) and macrofauna (>1 mm). We used a combination of eigenvector based multivariate analyses and variance partitioning on this multiple set of variables to identify which environmental variables likely drive the benthic diversity patterns. Mangrove early development increased the alpha and beta diversities of the infaunal communities for the two size classes. A total of 20–30% and 7–12% of the beta diversity are explained by linear and nonlinear spatial variables, respectively. However, 7% to 9% of the variance partioning could be determined by other biotic/abiotic variables, biological interactions or neutral processes, not described here. This study has highlighted the necessity of taking into account mangrove dynamics at suitable spatial scales for benthic biodiversity evaluation and mangrove management or restoration plans.


Introduction
Mangroves are one of the most productive ecosystems along subtropical and tropical coastlines.These systems are threatened by a range of anthropogenic pressures, including deforestation, pollution, and urbanization, and by climate changes, including recurrent droughts, sea level rise, cyclones, and coastal erosion [1].It has been estimated that 35% of the world's mangroves are degraded or have been lost [2].Key ecosystem services provided by mangroves (i.e., feeding areas, nurseries, blue carbon, sediments and contaminants retention, bioremediation), and strongly aided by benthic biodiversity [3][4][5], are thus weakened [6].One of the main difficulties for mangrove management and restoration strategies, however, is the lack of knowledge regarding mangrove ecosystem functioning and its spatial and temporal variabilities.This lack of knowledge often leads to a selection of small numbers of sites based on poorly defined criteria/variables, and, ultimately, to monospecific plantations that reduce biotic interactions between animal and plant species, and the diversity of ecological functions [7,8].The efficient sustainability of mangroves requires better characterization of these spatial and temporal variabilities as a function of regional characteristics, and in response to the various environmental and anthropogenic pressures that constraint mangrove functioning [9].
Mangrove forest structure is characterized by specific zonation along the intertidal gradient, each plant belt dominated by a tree species [10].The composition and structure of macrobenthic faunal communities change with this succession of the mangrove forest, in response to variations in substrate elevation, local geomorphology, salinity, predation, food resources and hydroperiod within each plant zone [5,[11][12][13][14][15].This pattern suggests that several benthic species may be good indicators of successional stages [16,17].The main taxa of mangrove benthic assemblages within sediments are annelids, crustaceans, mollusks and insects, but their respective abundances and species richness are highly variable between and within the same mangrove area due to local habitat heterogeneity [3,18].Benthic species richness is known to be lower in mangroves than in adjoining habitats such as seagrass beds or sandy beaches [19].However, benthic infaunal biodiversity is often under-estimated in mangroves and not easily compared both among mangrove zones and with adjoining habitats because of different sampling strategies and taxonomic resolutions [5].Despite low benthic species richness in mangroves, communities are composed of highly abundant small sized individuals [20].Indeed, Dittmann [21,22] emphasized the need to consider mesofauna (0.25 < x < 1 mm; [23]) in addition to macrofauna (x > 1 mm) to describe the whole benthic infaunal assemblage in tropical intertidal mudflats and mangroves.
The French Guiana (FG) coast in northeastern South America is controlled by a mud supply from the Amazon River and it's alongshore transport.In addition to coastal FG, the 1500 kilometer-long muddy tract includes the coast of Brazil west of the mouths of the Amazon River, the coasts of Surinam and Guyana, and the mouths of the Orinoco River in eastern Venezuela, collectively referred to as the Guianas coast.Mangrove distribution and dynamics on the Guianas coast are constantly affected by an alternation of accretion and erosion phases determined by the onshore welding of mud banks and their reworking by waves [24][25][26].Shore attached mud banks periodically migrate along the coast (i.e., mobile muds), become topographically elevated during continuous mud accretion, and consolidate, allowing pioneer mangroves to colonize the substrate and develop extremely rapidly.Vast monospecific Avicennia germinans dominated forests have developed in FG, contrary to other worldwide mangroves that display specific tree succession over short (<km) distances [27,28].Rates of mud consolidation and mangrove development previously measured (tree growth up to ~2 m.yr −1 ) are probably the fastest worldwide [27,29].
Previous studies of the benthic infaunal communities in the mobile muds and newly accreted bare sediments along the north Brazil, FG, and Surinam coasts have shown a low taxonomic richness of macrofauna (>1 mm) [30][31][32][33]).Two recent studies in FG demonstrated that mesofauna (>0.25 mm) was abundant in consolidated bare muds, with bioturbation activity as high as that of macrofauna in young mangroves [34,35].However, the spatial distribution of benthic meso-and macrofaunal communities along the mangrove age gradient and the environmental variables that drive them remain to be determined.This information is essential to better characterize the functioning of this ecosystem and necessary for designing sustainable coastal zone management under the influence of the Amazon River discharge.This knowledge is particularly critical, given the fact that the FG coast is still largely exempt from the impacts of human activities and thus could, provide insight into benthic community responses to natural disturbance cycles.Unlike FG, where mangrove removal by humans has been extremely limited [36], Surinam, and especially Guyana, have been drastically affected, leading to enhanced land erosion and flooding, and the natural loss of resilience [37,38].Mangroves in neighbouring Brazil are also threatened by deforestation and pollution [39].These regional threats are exacerbated by climate change and predicted sea level rises [40].
Although FG mangroves are Avicennia germinans dominated, benthic biodiversity is expected to respond rapidly to mangrove maturation at both the α and β diversity levels.α diversity refers to the number of species/taxa in communities, whereas β diversity relates to changes in community composition and structure among localities or along a gradient [41].We hypothesized that the natural mangrove development associated with changes in the physical and biogeochemical properties of the sediments would increase both benthic α and β diversities.Specifically, this study aimed to: (1) characterize the vegetation and sediment variables, from bare sediment though young mangrove stages, (2) identify which environmental variables likely drive benthic community structure along the mangrove growth gradient, and (3) describe how benthic communities are modified along this gradient on the FG coast, under the influence of mud supply from the Amazon and periodic erosion.

Study Area and Sampling Sites
The present study was conducted during the dry season months of October and November 2013 on a mud bank located northwest of the Sinnamary estuary (5 • 28 41 N, 53 • 02 05 W; Figure 1a).We used a time series of PLEIADE satellite images (0.5 m.pixel −1 ; https://www.intelligence-airbusds.com/geostore) (accessed on 19 September 2013) to observe the accretion and elevation processes of this mudbank from a bare mudflat stage that started in 2009.As the mud consolidated, the substrate elevation and the aerial exposure time increased.Mangroves colonized the intertidal and consolidated mud in early 2012 and the oldest trees were just under two years when the study commenced.The saplings spread southwestward onto the adjacent consolidated bare sediments.
Tides are semidiurnal in the area, with high spring and neap tide levels of up to 3.2 m and 2.5 m, respectively, and a mean tidal range of 1.68 m (SHOM, 2013).The entire study area was inundated during spring tides whereas the higher substrates remained aerially exposed for up to 4-5 days during neap tides.
The subtidal and mobile parts of the mud bank were sampled at four stations (AA, J, K, L) (Figure 1b).The mangrove area was investigated during low spring tides.A ca. 450 m × 50 m grid was set up, and consisted of 9 equidistant (ca.50 m) shore parallel transects (from A to I; Figure 1c).Three stations, about 25 m apart, were investigated per transect (stations 1, 2 and 3; Figure 1c).As the three stations of each transect were parallel to the shoreline, they were subjected to the same tidal level simultaneously.

Topographic Measurements
In order to evaluate the topographic modifications that match mangrove development, we built a digital elevation model (DEM) using a Real-Time Kinematic Differential GPS (RTK-DGPS) Trimble R8.Measurements were obtained in the Universal Transverse Mercator (UTM) 22 North zone combined with the WGS 84 datum, using the Earth Gravitational Model 1996 (EGM96) geoid as the vertical datum.We carried out 51 measurements over the grid (from transects A to G) within the same day in early October, using a mobile GNSS antenna mounted on a kayak at high tide to prevent disturbance of the mud.High tree density precluded measurement at transects H and I, where the elevations were estimated.
The topographic values were meshed using the classical Delaunay triangulation method in GIS software ESRI ArcGIS desktop 10.1, and interpolated in the DEM on a grid

Topographic Measurements
In order to evaluate the topographic modifications that match mangrove development, we built a digital elevation model (DEM) using a Real-Time Kinematic Differential GPS (RTK-DGPS) Trimble R8.Measurements were obtained in the Universal Transverse Mercator (UTM) 22 North zone combined with the WGS 84 datum, using the Earth Gravitational Model 1996 (EGM96) geoid as the vertical datum.We carried out 51 measurements over the grid (from transects A to G) within the same day in early October, using a mobile GNSS antenna mounted on a kayak at high tide to prevent disturbance of the mud.High tree density precluded measurement at transects H and I, where the elevations were estimated.
The topographic values were meshed using the classical Delaunay triangulation method in GIS software ESRI ArcGIS desktop 10.1, and interpolated in the DEM on a grid (resolution of 2 m).The DEM portion at transects H and I was estimated from linear regression interpolations yielded by stations 1, 2 and 3, along the gradient from D to G (to avoid overestimation) with Z = aD + b, where Z is the elevation (m), D the distance in the slope direction (m), and a and b the regression coefficients.

Hydrodynamic Measurements
Since small differences (e.g., 10 cm) in substrate elevation with mangrove development induce a difference in emersion/immersion time between the bare sediments and the vegetated sediments, we calculated flood durations for each transect.The mean flood durations during spring and neap tidal cycles were evaluated from a combination of in situ water level measurements at Sinnamary and offshore tide gauge records at the Salvation Islands, 50 km offshore southeast of the Sinnamary estuary.The latter were downloaded from the free data portal of the French Marine Hydrographic Office (SHOM) (https://data.shom.fr,19 September 2013), whereas the in situ water levels were recorded with a NKE Ltd.ALTUS submersible altimeter deployed on the consolidated mudflat in the pioneer mangrove stage, 100 m south of the study grid (5 • 28 33 N, 53 • 02 03 W).Data were collected from 4th October to 11th November 2013 at a frequency of 1 min.
A mean offset of −0.425 ± 0.012 m was detected between the water levels measured in situ (Sinnamary) and those routinely recorded by the Salvation Islands tide gauge.
Local water elevations at any tidal phase were retrieved by subtracting this offset from the water level tables (15 min step) calculated at the Salvation Islands using the following: , where φ (second) is the tide phase; T 1 , and T 2 , the tidal cycle initiation and ending times, respectively; t the tidal time (from the cycle initiation) corresponding to the sought water level h(t); and H 1 and H 2 the high and low tide water levels, respectively.A half tide cycle of 6 h 10 min was considered from mean flooding/immersion duration recorded by a pressure sensor and in accordance with the classical semidiurnal tide duration.Finally, the flooding duration at each station was evaluated from the number of 15-minute periods, during which the water level was above the elevation of the studied stations.

Vegetation Structure
At the end of the benthic sampling phase (3rd-5th November), we described the vegetation structure of each stand age during low tides.We defined the vegetation plots in terms of number and surface area as a function of stem density, as previously recommended and applied [35,42]: a large plot area of 27,500 m 2 covering the three stations of transects A to D with bare sediment, one plot of 400 m 2 per station for transects E and F, and one plot of 100 m 2 per station at transect G. Three 1-m 2 plots per station were inventoried for transects H and I because of the very high tree density.Tree species and abundances, stem diameters (cm), tree heights (m) and the number of pneumatophores were measured in each plot and then normalized to a square meter scale.Tree height and diameter at breadth height were measured and used for above ground biomass (AGB) calculation following allometric equations (expressed in ton dw .ha−1 ; [27,42]).Stem diameter was measured at breast height (~1.3 m) for trees > 2 m, and at half height for trees < 2 m.Mangrove litter (i.e., fresh and decaying leaves, flowers, twigs) was hand collected at the surface of each plot and then stored in plastic bags to evaluate plot biomass standardized to one square meter.

Sediment Bulk Variables
Abiotic (i.e., grain size, sediment bulk density, water content, porosity and salinity) and biotic variables related to the fauna's food (i.e., total carbon (TC), total nitrogen (TN), photosynthetic pigments and benthic microalgae composition) were sampled within the intertidal sediments.Variables were measured in the surface sediments of each station on the consolidated mudbank (stations A to I).We collected three surface sediment samples (3 cm Ø cut off syringes, 1 cm depth) within 10 m of each GPS marked station (Figure 2) for all variables except for grain size (1 sample) and microalgae composition (2 cores, 1.6 cm Ø cut off syringes, 1 cm depth).The 1-centimeter depth maximizes the quantification of benthic microalgal biomass which does not migrate deeper in sediment during high tides [43].The samples dedicated to carbon, nitrogen and pigment analyses were frozen immediately in the field (dry ice).Those dedicated to the composition of microphytobenthos were fixed with 1 mL of Lugol and stored in the dark at 4 • C in a cooler.Due to accessibility and in order to prevent sediment disturbance, sediment samples were collected at high tide from a kayak at all sites (13th-17th October) except for transect H (19th October) which was only accessible on foot at low tide.Pore-water salinity was determined from three additional surface sediment samples taken at each station at low tide (LT), between LT −3 h and LT + 3 h.
(3 cm Ø cut off syringes, 1 cm depth) within 10 m of each GPS marked station (Figure 2 for all variables except for grain size (1 sample) and microalgae composition (2 cores, 1. cm Ø cut off syringes, 1 cm depth).The 1-centimeter depth maximizes the quantificatio of benthic microalgal biomass which does not migrate deeper in sediment during hig tides [43].The samples dedicated to carbon, nitrogen and pigment analyses were froze immediately in the field (dry ice).Those dedicated to the composition of microphytoben thos were fixed with 1 mL of Lugol and stored in the dark at 4°C in a cooler.Due to acces sibility and in order to prevent sediment disturbance, sediment samples were collected a high tide from a kayak at all sites (13th-17th October) except for transect H (19th October which was only accessible on foot at low tide.Pore-water salinity was determined from three additional surface sediment samples taken at each station at low tide (LT), betwee LT −3 h and LT + 3 h.Surface sediment has the advantage of being quickly sampled and replicated ove large areas.However, some surface variables such as salinity can be tide dependent.I order to check how the abovementioned variables collected at the sediment surface rep resented benthic habitat over depth, we compared surface sediment variables with dept averaged sediment variables that integrate the tidal processes.Three additional sedimen cores (10.4 cm inner diameter, 18 cm height, 84.9 cm 2 area) were collected at transects A E and I (i.e., one core per station).These transects were assumed to be representative o the three main mangrove ages ("0" month, 6 months, 2 years old).Core sediments wer sliced into 0.5-centimeter layers from 0 to 5 cm depth, 1-centimeter layers from 5 to 10 cm depth, and 2-centimeter layers from 10 to 18 cm depth.One additional core per selecte transect was collected for deep pore-water salinity and divided into 0-10 cm and 10-2 cm horizons.For each core, sedimentary variables were averaged over depth.
The wet and dry (60 °C, 24 h) sediment densities (m/V) were obtained by weighin the wet and dried sediment (m), respectively, of a known volume (V) of the initial we sample.Water content, WC (%), was calculated from the weights of sediment using: WC = [(wet weight-dry weight)/wet weight] × 100.Then, the sediment porosity was est mated as (wet density × water content)/100 (the small salinity correction was ignored Pore-water salinity was measured using a refractometer (±1 ppt; practical salinity scale after extracting a drop of filtered interstitial water from the sediment [44].Grain size frac tions in the range 0.05 µm to 2000 µm were analyzed with a Malvern Mastersizer-S200 (Malvern Instruments Ltd., Malvern, UK) after deflocculating the aggregated particle [45].The median grain sizes (diameter) were separated as clay, mud (very fine, fine an medium silts), coarse silt, very fine sand and sand.
The biomass of mangrove litter collected in each vegetation plot was estimated b weighing collected dried matter to constant weight.Total carbon (TC) and nitrogen (TN Surface sediment has the advantage of being quickly sampled and replicated over large areas.However, some surface variables such as salinity can be tide dependent.In order to check how the abovementioned variables collected at the sediment surface represented benthic habitat over depth, we compared surface sediment variables with depth averaged sediment variables that integrate the tidal processes.Three additional sediment cores (10.4 cm inner diameter, 18 cm height, 84.9 cm 2 area) were collected at transects A, E and I (i.e., one core per station).These transects were assumed to be representative of the three main mangrove ages ("0" month, 6 months, 2 years old).Core sediments were sliced into 0.5-centimeter layers from 0 to 5 cm depth, 1-centimeter layers from 5 to 10 cm depth, and 2-centimeter layers from 10 to 18 cm depth.One additional core per selected transect was collected for deep pore-water salinity and divided into 0-10 cm and 10-20 cm horizons.For each core, sedimentary variables were averaged over depth.
The wet and dry (60 • C, 24 h) sediment densities (m/V) were obtained by weighing the wet and dried sediment (m), respectively, of a known volume (V) of the initial wet sample.Water content, WC (%), was calculated from the weights of sediment using: WC = [(wet weight-dry weight)/wet weight] × 100.Then, the sediment porosity was estimated as (wet density × water content)/100 (the small salinity correction was ignored).Pore-water salinity was measured using a refractometer (±1 ppt; practical salinity scale) after extracting a drop of filtered interstitial water from the sediment [44].Grain size fractions in the range 0.05 µm to 2000 µm were analyzed with a Malvern Mastersizer-S2000 (Malvern Instruments Ltd., Malvern, UK) after deflocculating the aggregated particles [45].The median grain sizes (diameter) were separated as clay, mud (very fine, fine and medium silts), coarse silt, very fine sand and sand.
The biomass of mangrove litter collected in each vegetation plot was estimated by weighing collected dried matter to constant weight.Total carbon (TC) and nitrogen (TN) were analyzed by combustion at 930 • C using a CHN carbon analyzer (FLASH-2000, Thermo Fisher Scientific Inc., Milan, Italy).Due to the lack of carbonates on the FG coast, the organic fraction largely dominates total carbon [29].TC and TN were thus used as proxies of the organic matter, and the molar TC:TN ratio (mol:mol) was calculated to infer a proxy of the refractory versus labile nature of the organic matter.Sedimentary pigment concentrations were determined by fluorimetry [46].The chlorophyll-a to phaeopigment ratio was used to assess the degradation state of phytodetritus.
Samples for identification and counting of microphytobenthos were collected using a 10 mL syringe (diameter 1.6 cm) cut at its end.The surface sample (0-1 cm of sediment) was then fixed immediately with 1 mL of Lugol's and kept cold and in the dark.The composition of the microphytobenthos was determined later in the lab, after staining the sample with a solution of rose bengal diluted in filtered seawater (30 µL/3 mL).Counts were performed using an inverted microscope (Zeiss Observer A1), on one or more sedimentation dish diameters for the most abundant species, or on the entire surface for the rarest [47,48].The number of cells counted was ratioed to the total volume of sediment, and to the area sampled, to determine the total concentrations of each species per m 2 of sediment.The identification of diatoms was carried out using SEM after cleaning the frustules (with H 2 O 2 at 60 • C).

Benthic Infauna
For the study of the benthic fauna, three additional sediment cores (10.4 cm diameter, 18 cm depth, 84.9 cm 2 area) were collected at each of the 3 stations in each transect along the mangrove grid.Cores were pooled within stations, giving a single sampling unit with an area of 254.7 cm 2 per station (Figure 2).Sampling was performed depending on site accessibility during high tides (transects A, B, C, F, G and I) and low tides (transects D, E and H) on 9th-19th October.In order to distinguish the infaunal size structure (mesovs.macrofauna), sediments were sieved through 0.25 mm and 1 mm meshes, and the organisms were fixed in a 4% buffered formaldehyde solution until identification.
Meso-and macro-organisms were sorted and identified to the lowest possible taxonomic level (hereafter referred to as the operational taxonomic unit, OUT, based on the morphological criteria) and counted.Image based measurements of organism length and width (Visilog ® software, Noesis, France, 1 px equals to 1.0-12.4µm depending on the magnification used) taken under a binocular microscope were converted into biovolume (mL) by relating the smallest organisms to geometrical forms and by immersing the largest organisms in graduated filled containers [35].For the rarest taxa, the biovolume of each individual was recorded, whereas for the dominant taxa (oligochaetes, ostracods, copepods, tanaids and nematodes), average individual biovolume was calculated from the total biovolume measurement of 100 individuals per size range (0.25 < x < 1 mm and x > 1 mm).For small sized organisms that cannot be weighed, the measurement of biovolume is useful as a proxy for biomass (e.g., meiobenthos in [49]).For unified comparisons with data from the literature, abundance and biovolume values were standardized per square meter for the entire benthic communities, and specifically for the benthic meso and macrofauna.
Additional small sized sediment cores (3 cm diameter, 5 cm depth) were collected for a preliminary analysis of the meiofauna structure (>62 µm) at a few consolidated mud stations, A (6 cores), B (3 cores), E (1 core), F (9 cores), G (9 cores), H (1 core), I (3 cores), and at mobile mud stations AA (1 core), J (3 cores), K (3 cores) and L (3 cores).Sediments for the meiofauna analysis were fixed in a 10% buffered formaldehyde solution until identification.Before sorting, meiofauna samples were stained with rose bengal, washed onto a 0.062 µm mesh sieve, and sorted under a dissecting microscope to the class level.Meiofaunal abundances were normalized to the weight of dry sediment after noting porosity ranges [31].

Statistical Analyses
A principal components analysis (PCA) on the correlation matrix was conducted to determine which environmental variables (among tree biomass and density, pneumatophores and litter biomasses, substrate elevation, porewater content, salinity, porosity, sediment grain size, pigments, Chla:Phaeo, TC and TN contents, TOC:TN and C-Chla:TC ratios) explained the most variability between the transects.Then, statistical differences in the environmental variables were tested between transects with PERMANOVA analysis with 999 permutations and threshold of probability at 0.05 [50].The differences in microphytobenthic composition (diatoms versus cyanophycae) between transects were tested using a nonparametric one way analysis of variance (Kruskal-Wallis) followed by a post hoc Wilcoxon test for pairwise comparisons.
Regarding the environmental variables sampled at the surface (i.e., C:N and Phaeo:Chla ratios, salinity, and median grain size), we tested their reliability in representing the infaunal habitats during high versus low tides.We examined the hypothesis of "no effect of sampling time on surface sediment parameters" in these mangrove stations.The variables that fulfilled the normality and homoscedasticity criteria were tested with an ANOVA and log transformed when necessary (C:N, log transformed Phaeo:Chl-a), while the others were analyzed with a nonparametric Kruskal-Wallis test setting the tide effect as a fixed factor (two conditions: low tide for transect H and high tide for transects G and I).The patterns of surface sediment variables on the one hand, and sediment depth averaged variables on the other, were compared both as a function of transect using graphical analysis and analyses of covariance (ANCOVA).Nonparametric ANCOVAs [51] were chosen to verify the hypothesis of curve parallelism.ANCOVA was, however, not applied to salinity because there was no core replication per transect.
The effects of transects were tested on the α diversity for the entire benthic infauna (i.e., total OTUs richness, Shannon's entropy and Pielou evenness (J) indexes, density and biovolume) by using parametric one way variance analysis (ANOVA).We also compared the α diversity per size class (mesofauna versus macrofauna) depending on the transect, by using two way variance analysis (ANOVA).All of these ANOVA were followed by post hoc Tukey tests.
A permutational multivariate analysis of variance, using a Bray-Curtis dissimilarity index (PERMANOVA; [50]), based on density and biovolume per matrix of OTUs, was computed to test the multivariate response of OTU assemblages to the transect location after verifying the multivariate homogeneity of group dispersions.
In order to visualize the overall β diversity and the contribution of OTUs in the differentiation of stations among transects, two principal components analysis (PCA) based on total density and biovolume were performed on the Hellinger transformed data to avoid giving excessive weight to rare OTUs [52].
Redundancy analysis (RDA) was used to quantify the effect of spatial structure and abiotic and biotic environmental variables on infaunal β diversity.Distance based Moran's eigenvector maps (dbMEMs; [53]) and variation partitioning [54], were combined in order to isolate the specific and shared spatial and environmental variance fractions of abundance and biovolume infaunal communities in terms of beta diversity.Only the dbMEMs corresponding to positive spatial autocorrelation were used as spatial variables.Any significant linear spatial trend was first identified and removed with multiple linear regressions.Redundancy analysis was used to model beta diversity in terms of spatial and environmental variables.For each subset (dbmems, abiotic and biotic) of explanatory variables, a global model was first computed and subsequently tested.If significant, and in order to optimize the predictive power of the model [55], forward selection of explanatory variables (dbMEMs, biotic and abiotic environmental variables) was carried out in each subset according to [56] (see Appendix A for more details).All analyses were conducted within the R environment (R Development Core Team 2021; [57]) and relied on the FactoMineR [58], vegan [59], adespatial [60] packages.

Environmental Variables
The environmental PCA analysis explained 68.7% of the variability among transects along the mangrove maturation gradient: the first axis explained 43.6% and the second axis 25.1% (Figure 3).The most important variables defining the first axis were TN, Cchla:TOC, TC, chla, chla:phaeo, elevation, phaeo, salinity, and grain size and TC:TN, while pneumatophore density, tree density, tree biomass, litter biomass and elevation defined the second axis.Transects were separated into three groups based on PCA results: group 1 (transects A, B, C), group 2 (transects D, E, F, G), and group 3 (transects H and I).The PERMANOVA analysis showed that environmental variables were significantly different between the groups of transects (p < 0.001).

Environmental Variables
The environmental PCA analysis explained 68.7% of the variability among transects along the mangrove maturation gradient: the first axis explained 43.6% and the second axis 25.1% (Figure 3).The most important variables defining the first axis were TN, Cchla:TOC, TC, chla, chla:phaeo, elevation, phaeo, salinity, and grain size and TC:TN, while pneumatophore density, tree density, tree biomass, litter biomass and elevation defined the second axis.Transects were separated into three groups based on PCA results: group 1 (transects A, B, C), group 2 (transects D, E, F, G), and group 3 (transects H and I).The PERMANOVA analysis showed that environmental variables were significantly different between the groups of transects (p < 0.001).Group 1 included transects with the lowest values of salinity, elevation, C-chla:TOC and chla:phaeo ratios and the highest values of sediment water content, median grain size and TC:TN ratio (Figures 3-5).Group 2 exhibited the highest chla and TC contents, and C-chla:TOC and chla:phaeo ratios (Figures 3-5).Group 3 was associated with mangrove biological variables (pneumatophore density, litter biomass, tree biomass and density) and increasing substrate elevation, salinity, and pheophytin content (Figures 3-5).Group 1 included transects with the lowest values of salinity, elevation, C-chla:TOC and chla:phaeo ratios and the highest values of sediment water content, median grain size and TC:TN ratio (Figures 3-5).Group 2 exhibited the highest chla and TC contents, and C-chla:TOC and chla:phaeo ratios (Figures 3-5).Group 3 was associated with mangrove biological variables (pneumatophore density, litter biomass, tree biomass and density) and increasing substrate elevation, salinity, and pheophytin content (Figures 3-5).

Validity of the Selected Surface Sediment Variables
Surface sediment biogeochemical variables (Phaeo:Chl-a, C:N; ANOVA, p > 0.3), and median grain size (ANOVA, p > 0.9) did not significantly differ between the low and high tides within transects A, E, I. Graphical and ANCOVA analyses showed that these variables, measured at the sediment surface and depth averaged over 18 cm, exhibited similar trends between the three transects.The slopes of the regressions of the surface and depth averaged variables as a function of mangrove stages did not violate the assumption of parallelism for median grain size and C:N ratio (p > 0.05), i.e., both curves had the same shape (Figure S2).The curves representing the surface and depth averaged Phaeo:Chl-a were not parallel (p < 0.05), though they showed relatively comparable patterns along the mangrove development gradient (Figure S2).ANCOVA was not applied to salinity, but the graph (Figure S2) clearly showed the same pattern at the sediment surface and at depth between the three transects.Substrate elevation decreased slightly (from 1 to 5 cm) along a north-south axis (i.e., from stations 1 to 3) in most of the shore parallel transects (A, B, C, F and G) and showed a significant increase in the order D, E, F, G (60 cm) and H, I (83 cm).Transects of groups 2 (D to G) and 3 (H, I) were always exposed to air during neap tides (nil flood duration) whereas group 1 (A to C) was flooded for 2.5 to 3.5 h per tidal cycle.During spring tides (i.e., sampling time), the flooding duration decreased from group 1 (5-5.5 h) to group 2 (4-4.5 h), and group 3 (3-4 h).
The total abundance of microphytobenthos (MPB) significantly increased in transects of groups 1 and 2 (p < 0.001, KW).This trend was explained by higher cyanobacteria abundances in the vegetated sediments (p < 0.0001, KW; Figure S1A, Supplementary Materials), except for transect I, whereas diatom abundances remained similar along the mangrove maturation gradient (p > 0.12, KW; Figure S1A).Relative abundances of diatoms were higher than those of cyanobacteria in group 1 (>90%), whereas cyanobacteria dominated in groups 2 and 3 (>80%; Figure S1B).The same trends were observed for the MPB biovolumes (data not shown).

Validity of the Selected Surface Sediment Variables
Surface sediment biogeochemical variables (Phaeo:Chl-a, C:N; ANOVA, p > 0.3), and median grain size (ANOVA, p > 0.9) did not significantly differ between the low and high tides within transects A, E, I. Graphical and ANCOVA analyses showed that these variables, measured at the sediment surface and depth averaged over 18 cm, exhibited similar trends between the three transects.The slopes of the regressions of the surface and depth averaged variables as a function of mangrove stages did not violate the assumption of parallelism for median grain size and C:N ratio (p > 0.05), i.e., both curves had the same shape (Figure S2).The curves representing the surface and depth averaged Phaeo:Chl-a were not parallel (p < 0.05), though they showed relatively comparable patterns along the mangrove development gradient (Figure S2).ANCOVA was not applied to salinity, but the graph (Figure S2) clearly showed the same pattern at the sediment surface and at depth between the three transects.

Infaunal Alpha Diversity
A total of 56 infaunal taxa belonging to 16 classes were identified across all stations (Table S1).Total OTU richness and the Shannon index significantly increased along the mangrove maturation gradient, reaching the highest values at the vegetated stations H and I (p < 0.02, Figure 6a,b).OTU richness was higher for mesofauna at all stations.At the low elevation stations (A to D), the Shannon index was higher for mesofauna (1.5-2).Intermediate stations (E-G) were characterized by a similar Shannon index between meso and macrofauna, while values were significantly higher (>2) for macrofauna at vegetated stations.The Pielou index (J) of the total infauna remained similar along the mangrove transects, and it only differed in terms of the infaunal size category (p < 0.05, Figure 6c).The Pielou index for mesofauna was significantly higher than that of macrofauna in low elevation bare mud transects (A-D), whereas the opposite trend was found in the upper vegetated transects (H-I).Differences in the Pielou index were not significant between the two infaunal size categories in intermediate transects (E-G).Only the macrofaunal Pielou index significantly increased along the mangrove maturation gradient, from stations A to F (similar values from stations F to I).The total density and biovolume of the infaunal communities (meso + macro) ranged from 25 to 300.10 3 ind.m−2 and from 5 to 130 mL.m −2 , respectively, across transects (Figure The total density and biovolume of the infaunal communities (meso + macro) ranged from 25 to 300.10 3 ind.m−2 and from 5 to 130 mL.m −2 , respectively, across transects (Figure 6 d,e).Total densities significantly decreased from the bare sediments to the most vegetated sediments (F = 12.63; p < 0.0005; Figure 6d).Mesofauna significantly dominated macrofauna in density in each transect (Anova, F = 101; p < 0.0005).Total biovolumes differed significantly among transects (ANOVA, F= 25.72; p < 0.0005; Figure 6e), and were lowest at the transects of group 2 (D, E, F; Tukey, p < 0.05).Biovolumes remained similar between the two size categories within a given transect.
The total abundance of the overall meiofauna ranged from 0.5 to 130 ind.g −1 dry weight of sediment, with the lowest values at the mobile mud stations (AA, J, K, L; KW test: 16.769, p < 0.005).Nematodes were the exclusive taxa in density in AA and L stations, dominant in B and H stations and equal to polychaete taxa in the station I (Figure S4A).The relative proportions of tanaids and ostracods increased at stations A, B, J and K (Figure S4A).Vegetated sediments were characterized by copepods mainly in transect E, and juvenile polychaetes and insects.Juvenile decapods were only observed at station K (Figure S4A).The tanaid population of the intertidal bare sediments (A-B) was dominated by individuals with body length comprised between 1 and 5 mm, whereas those of the mobile mud stations (J-K) were smaller (<2 mm in J; 2 to 3.5 mm in K) (Figure S4B).
PCA performed on total biovolumes explained 49.2% of the variability among sampling stations: the first axis explained 33.6% and the second axis 15.6% (Figure 7b).The PCA first axis still opposed the biovolumes of H. spaansi and Ostracoda spp. of group 1 transects (A to C), to those of Ceratopogonidae and Namalycastis sp.1 in the vegetated transects (groups 2 and 3; D to I).The second axis mainly explained the variability among the vegetated stations due to high biovolumes of the malacostracean decapod Uca cumulanta and of adult insects (Coleoptera sp.4) at I transect.

Disentangling the Effects of Environment Variables on Infaunal Multivariate Responses
Appendix A specifies how the variables were selected prior to the variance partitioning analysis.The variable selection process ended up with a dataset of 10 explanatory variables (spatial: dbMEMs 1, 5, 6 and 7, Figure S4; environmental: elevation, pneumatophore density, C:N, Phaeo:Chl-a, salinity and median grain size) used in the variation partitioning analyses.

Disentangling the Effects of Environment Variables on Infaunal Multivariate Responses
Appendix A specifies how the variables were selected prior to the variance partitioning analysis.The variable selection process ended up with a dataset of 10 explanatory variables (spatial: dbMEMs 1, 5, 6 and 7, Figure S4; environmental: elevation, pneumatophore density, C:N, Phaeo:Chl-a, salinity and median grain size) used in the variation partitioning analyses.
Among the abiotic and biological variables, substrate elevation, grain size, salinity, TC:TN and Phaeo:Chl-a ratios, and pneumatophore density had little to no effect (<5%) on the benthic community structure, as they act synchronously with the spatial variables.The linear spatial structure, together with the abiotic and biological variables, jointly explained most of the variations in the infaunal structure (30% for abundances and 20% for biovolumes).Likewise, the nonlinear spatial structures and the abiotic and biological variables jointly explained a non-negligible part of the variations in the infaunal structure (12% for abundances and 7% for biovolumes).Part of the variations in infaunal structure (9% for abundances and 7% for biovolumes) is only attributable to nonlinear spatial structures (dbMEMs), without any relation to the measured variables (Figure 8).Linear spatial predictors (XY coordinates) significantly explained the total abundance variations (R 2 Adj = 0.34, F = 7.811; p < 0.001) and total biovolume variations (R 2 Adj = 0.24, F = 5.135, p < 0.001; 9999 perm.)along the south-north and east-west directions (Appendix A2).The south-north effect was stronger (abundances: R 2 Adj = 0.33, F = 13.887,p < 0.001; biovolumes: R 2 Adj = 0.24, F = 9.030, p < 0.001; 9999 perm.)than the east-west one (abundances: R 2 Adj = 0.13, F = 4.811, p = 0.0029; biovolumes: R 2 Adj = 0.12, F = 4.659, p < 0.001; 9999 perm.).We partitioned the beta diversity of benthic faunal communities according to four sets of explanatory environmental variables: linear spatial structures, nonlinear spatial structures (dbMEMs), abiotic and biological variables, (Figure 8).Together, these variables explained 67% (R 2 adj total = 0.67) and 45% (R 2 adj total = 0.45) of the variations in the total abundance and biovolume, respectively.
Among the abiotic and biological variables, substrate elevation, grain size, salinity, TC:TN and Phaeo:Chl-a ratios, and pneumatophore density had little to no effect (<5%) on the benthic community structure, as they act synchronously with the spatial variables.The linear spatial structure, together with the abiotic and biological variables, jointly explained most of the variations in the infaunal structure (30% for abundances and 20% for biovolumes).Likewise, the nonlinear spatial structures and the abiotic and biological variables jointly explained a non-negligible part of the variations in the infaunal structure (12% for abundances and 7% for biovolumes).Part of the variations in infaunal structure (9% for abundances and 7% for biovolumes) is only attributable to nonlinear spatial structures (dbMEMs), without any relation to the measured variables (Figure 8).The two correlation triplots of RDA (based on OTU abundance and biovolume data; Figure S5) shows how the selected biotic and abiotic variables constrained the infaunal structure.

Environmental Changes versus Young Mangrove Development
Mangroves grow rapidly along the French Guiana coast as an adaptative response to the natural and recurrent cycles of accretion/erosion phases [27,42].In our study area, mangrove growth began when the mud was consolidated at 600 g.L −1 , physically stable, and at a substrate elevation of 30 cm (transects D-F, group 2).These values are in agreement with those found previously [28].The oldest trees, composed solely of A. germinans, were 2 years old with a height of 8 m.This growth rate averaged 4 meters a year and was superior to the previous studies in FG (~2 meters a year) and is a worldwide record for mangrove forests.A previous study has shown that Avicennia growth and biomass are widely enhanced in low latitude areas with high precipitation amounts [61].Such conditions result in diminished water salinity (<30), a condition suitable for Avicennia to thrive.This high growth rate yields an above ground biomass of 145 t DM.ha −1 and density of 190,000 ind.ha −1 , for a substrate elevation of 45 cm (H-I, group 3).This elevation decreased the flooding time in the west-east direction (A-I, groups 1 to 3), thus enhancing porewater evaporation and salinity in the eastern part (H-I, group 3) of the study area.Seedling settlement and development started within the D-F transects of group 2. A previous study from the same area found that remineralization rates in initially deposited fluid muds had demonstrated that the initially deposited sediments were extremely nutrient rich [62].Thus, mangroves can take advantage of the high reactivity of the substrate when the deposits are physically stabilized (group 2).We surmise that mangrove growth is enhanced by sediment stability, dominant suboxic conditions (Fe rich The two correlation triplots of RDA (based on OTU abundance and biovolume data; Figure S5) shows how the selected biotic and abiotic variables constrained the infaunal structure.

Environmental Changes versus Young Mangrove Development
Mangroves grow rapidly along the French Guiana coast as an adaptative response to the natural and recurrent cycles of accretion/erosion phases [27,42].In our study area, mangrove growth began when the mud was consolidated at 600 g.L −1 , physically stable, and at a substrate elevation of 30 cm (transects D-F, group 2).These values are in agreement with those found previously [28].The oldest trees, composed solely of A. germinans, were 2 years old with a height of 8 m.This growth rate averaged 4 meters a year and was superior to the previous studies in FG (~2 meters a year) and is a worldwide record for mangrove forests.A previous study has shown that Avicennia growth and biomass are widely enhanced in low latitude areas with high precipitation amounts [61].Such conditions result in diminished water salinity (<30), a condition suitable for Avicennia to thrive.This high growth rate yields an above ground biomass of 145 t DM.ha −1 and density of 190,000 ind.ha −1 , for a substrate elevation of 45 cm (H-I, group 3).This elevation decreased the flooding time in the west-east direction (A-I, groups 1 to 3), thus enhancing porewater evaporation and salinity in the eastern part (H-I, group 3) of the study area.Seedling settlement and development started within the D-F transects of group 2. A previous study from the same area found that remineralization rates in initially deposited fluid muds had demonstrated that the initially deposited sediments were extremely nutrient rich [62].Thus, mangroves can take advantage of the high reactivity of the substrate when the deposits are physically stabilized (group 2).We surmise that mangrove growth is enhanced by sediment stability, dominant suboxic conditions (Fe rich sediments) associated to high sediment reactivity (i.e., remineralization rates), and by nutrients uptake (i.e., phosphorus), the concentrations of which are important in resuspended sediments and in the overlying water column according to previous studies in the same area [29,44,62,63].Such uptake may be more important during spring tides, when the vegetated sediments were completely flooded for hours, but should also be possible during neap tides, when the creeks are continuously flooded [64].Flooding processes favor exchanges of dissolved nutrients between sediments and water.This relationship between nutrients, mangrove substrate elevation and creek dynamics was hypothesized from the morphology and stand structure of west African mangroves but not directly demonstrated from measurements [65].Mangroves are typically plastic in their ability to use nutrients opportunistically as they become available [66].For example, phosphorus, a limiting nutrient for mangroves, was measured in low concentrations in the pioneer mangrove stage, probably due to repeated P desorption from the resuspended sediments to the coastal seawater, thus rendering this element available for primary production (i.e., microphytobenthos and mangroves; [63]).The same authors measured higher percentages of sedimentary org-P phosphorus in the mixed and mature mangrove stages, which correspond to the most important mangrove biomass.In the same way, the benthic microalgae, making up the photosynthetic component of the biofilm, developed to very great abundances (>100.10 3 cells mL −1 ) in near surface sediments, rapidly increasing the Chl-a content (~53 µg Chl-a g −1 DW).This surface biofilm was composed of diatoms within the bare low elevation sediments (A-C, group 1), as also identified in other FG mudbanks [67], but it was dominated by cyanobacteria in the mangrove swamp stabilized by vegetation (D-I, groups 2 and 3).Decreasing oxic conditions within mangrove sediments [29] could explain the cyanobacteria dominance, the latter often taking advantage under nutrient stress conditions [68,69].However, the biomass of benthic microalgae was found to decrease (~26 µg Chl-a g −1 DW) where the trees densities were highest (H-I, group 3), probably due to reduced sunlight at the sediment surface (shading).We hypothesize that the higher Chl-a:Phaeo and Chl-a:TC ratios at intermediate vegetated stages (D-G, group 2) indicate fresh benthic phytopigments and ungrazed litter.Conversely, the decrease in these ratios in the denser mangrove reflects a greater fraction of relatively degraded organic matter.Vegetation development was also associated with an increase in litter deposits and pneumatophore density.We surmise that enhanced root transpiration raised sediment porewater salinity (13 instead of 4 in the bare sediments).

Benthic Biodiversity Spatial Structure versus Young Mangrove Growth
Benthic biodiversity correlated with mangrove development.The linear spatial variables along the west-east axis, associated with environmental variables, explained most of the variation in infaunal structure (30% for the abundances and 20% for the biovolumes, Figure 8).More specifically, the gradual increases in substrate elevation, sediment aerial exposure time and porewater salinity, and vegetation and pneumatophore densities with mangrove development (Figure 3, Table S1), were the main variables that explained benthic community variations along this west-east axis (Figure S5A).These variables are usually important drivers in structuring mangrove benthic communities [5,13].Our study showed that they changed linearly along the west-east mangrove gradient.The linear modification of infaunal community structure along the north-south axis, however, is more difficult to explain with our data.It could be influenced by the study grid design that shifted northward the two most vegetated stations (within H-I transects, group 3), both of which have a distinct infaunal composition, compared to the other stations (Figure 1, Figure 5 and Figure S5B), and by the slight northward increase in substrate elevation (from 1 to 5 cm).
Next, our measurements showed that 7% to 12% of the beta diversity was determined by nonlinear spatial variables jointly with biotic and abiotic sediment variables (Figure 8).Here, the infaunal structure correlated with an oscillation model along the west-east axis, differentiating the three groups of transects (dbMEM 1, Figure S5C).This is consistent with the biological PCA analyses showing that the benthic abundances and biovolumes were clustered according to these same groups of transects (Figure 7).In these conditions, the benthic community structure was related to the quality of sediment organic matter (i.e., TC:TN, Chl-a:Phaeo and Chl-a:TC ratios), which changed between the groups of transects.
A remaining part of the beta diversity (7% to 9%) was only defined by nonlinear spatial structures (dbMEMs) and, thus, not explainable by the variables used in our analysis.Other biotic or abiotic variables related to microscale heterogeneity, biological interactions or neutral processes, not described here, could explain this tendency.
It is likely that the soft and brackish sediments of the western consolidated part of the mudbank (group 1, A to C) were suitable substrates for the burrowing activities of Halmyrapseudes spansii (tanaid, family Parapseudidae) and Ostracod spp.[35,[70][71][72].The latter would benefit from low wave induced shear stress resulting from the attenuation of wave energy by the vegetated mudbank in the eastern sector.The presence of tanaids and ostracods only in the meiofaunal samples from stations J and K led us to conclude that the sediments at these stations were more stabilized than those at stations AA and L. These last stations were located in the western part of the mud bank which was still likely to be highly mobile.This highlighted the strong sensitivity of tanaids and ostracods to the sediment consolidation state.Another study also noticed the absence of tanaid species in highly unstable and dynamic sediments along the FG coast [33].Being herbivores on attached algae or detritivores, tanaids and ostracods preferentially feed on diatoms [73][74][75], which is consistent with their distribution in the consolidated bare sediments.Studies in Egypt and Turkey [76,77] found high abundances of ostracods usually associated with well oxygenated sediments due to diatom photosynthetic activities, and to areas of coarse silt and fine sand, consistent with our measurement in the bare sediment areas.
Although we did not take any porewater oxygen measurements, it is likely that sediments from group 2 (D to G) were suboxic-anoxic, given the greater proportion of cyanobacteria relative to diatoms and the dominance of Oligochaeta spp.Previous studies have shown the ecological preferences and physiological adaptations of Oligochaeta spp. to oxygen poor environments [78,79].In this study, the dominance of oligochaetes may result from organic matter enrichment by microalgal mats [80,81] and probably by pioneer mangrove seedlings.Oligochaeta spp., first considered as deposit feeders, would also be able to feed on dissolved organic substances, particularly sugars and amino acids [79].Dissolved organic matter can be released as exopolysaccharides from the benthic biofilm [82,83] following high frequent deposition at immersed tidal periods during spring tides [84,85].Although algal mats are known to inhibit copepod larval settlement in other ecosystems [86], we argue that the copepods we observed in greatest abundances at the sediment surface of the early vegetated transects have adapted their grazing activity to concentrate on cyanobacteria [87,88].
Desiccation is known to be an important factor limiting polychaete populations in consolidated bare tidal flats bordering mangroves [18], likely explaining their absence in group 2 transects (D-G).Polychaetes were mainly represented in this study by Capitellidae and Nereididae within the most vegetated transects (group 3, H-I), where degraded organic matter enrichment associated with higher silt clay fractions and substrate elevation became favorable for the growth of polychaete larvae and juveniles and for their burial [15,18,89].The benthic assemblages in group 3 (H-I) differed from the others by having higher densities of insect larvae (Cerotopoginidae and Chironomidae) the development of which was favored by a low wave energy habitat, substrate humidity and litter detritus as food resources [90,91].
Finally, 7-9% of infaunal beta diversity was correlated with nonlinear spatial variables unrelated to the mangrove gradient (Figure 8).This variation is related to small scale spatial heterogeneity, i.e., changes in the morphology and elevation of the order of meters between creeks, flat areas, and depressions [64], but also changes in root density, canopy cover and food resources [12,14,92].Likewise, such spatial heterogeneity is known to influence crab distributions in the study area [34,64].In this study, this may explain the patchy distribution of some taxa of the mesofauna (i.e., ostracods and tanaids) and of the macrofauna (i.e., Uca cumulanta and Coleoptera sp.4), particularly in some vegetated sediments.For instance, roots provide hard substrates for specific arthropod assemblages such as insect larvae or crabs [12].The strong association found in this study between pneumatophores and annelid (Capitellidae and Nereididae) occurrences (Figure 7) is likely to be due to deep sediment oxidation by roots, thus favoring worm burial.Roots can also offer refuges from predators [93].Although we did not study biotic relationships, our findings are consistent with the conclusions of [94], that brachyurans, fishes and birds can also directly (by predation) or indirectly (by competition or habitat physical disturbance) impact the small sized infauna distribution.

Benthic Biodiversity Patterns in Highly Dynamic Amazon Influenced Mangroves
Low alpha diversity values (0.5 < ENS < 3; 0.2 < J < 0.7) were observed overall, as classically and naturally found in benthic communities of transitional waters in estuaries, deltas or lagoons [95].In tropical areas, these communities are usually homogeneous with high infaunal densities and the dominance of a few OTUs referring to pioneer like species [20].In our study, this was more marked within the intertidal bare mud for the macrofaunal community (ENS < 1; J < 0.4) than for the mesofaunal ones (ENS > 1.5; J > 0.5).We explain this difference by the omnipresence of tanaids (H.spansii) within the macrofauna of bare sediments, while the mesofauna was a mix, in variable abundances, of "real mesofauna" (mature in this size class; e.g., H. spansii), permanent meiofauna (e.g., nematodes, ostracods, copepods) and temporary meiofauna (juveniles of a few OTUs of macrofauna, e.g., bivalves or polychaetes).The tanaid species (>0.5 mm) of the Sinnamary estuary would be more resistant and opportunistically developed according to their smaller size and earlier adulthood [33].Previous studies [30,32] also showed low macrofaunal diversity composed by small sized (mostly < 10 mm in length) polychaetes, bivalves, gastropods and a large fraction of tanaids in bare intertidal mudflats in Surinam and FG.They attributed the low diversity and small size of the macrofauna to unstable bare sediments along the Guianas coast.However, these authors analyzed macrofauna > 0.5 mm without taking into account the 0.25-0.5 mm size fraction, thus underestimating the mesofaunal component.Yet, another study in the Amazon influenced mobile muds found more than 50% of all infauna in the size range 0.3-0.5 mm while a large part of the remaining benthic biomass was attributed to the meiofauna and bacteria [31].The present study confirmed the dominance of mesofauna (>74%) and emphasized the fact that the mesofaunal communities (>0.25 mm) were more heterogeneous and diverse than the macrofaunal ones (>1 mm), and likely more adapted to such unstable mudbanks.It has been shown that intertidal bare muds preceding mangrove development in Australia are characterized by a more diverse mesofauna in the mudflat than at mangrove sites [21,22].Since we did not identify the most abundant mesofauna OTUs (nematodes, ostracods, copepods) at the level of families, genus or species, it is likely that we underestimated their species richness within FG bare muds.We know, for instance, that several nematode genera are specific to adult mangrove stands in another FG area (Cayenne estuary), thus reflecting the influence of differences in physical and biogeochemical composition of mangrove sediments [96].
The increases in the alpha diversity (OTU richness, diversity indices) and beta diversity (variance of the community, i.e., its heterogeneity) with mangrove growth emphasized the fact that macrofaunal communities become more heterogeneous at vegetated transects.This trend was explained mainly by the appearance of new OTUs with lower abundances but large sizes (Chironomidae and Ceratopogonidae larvae, Coleoptera sp., Polychaeta, Namalycastis sp., Uca sp.).The same OTUs occurred within the mesofauna, indicating the presence of juvenile macrofauna (low biovolumes).Mangrove roots and mudflat geomorphological variations associated with creeks and depressions created a variety of ecological niches including crab burrows and mounds [5,97], enhancing benthic diversity [98].Habitat complexity also favored the multiplicity of food resources for the invertebrates (litter, cyanobacteria, diatoms), their trophic modes and the associated food webs.We suppose this was suitable for the establishment of specialist species (longer life cycle species generally characterized by larger body size and low densities) as observed in subtropical mangroves of the Persian Gulf [15], but also for maintaining species with a wide range of life strategies [99].Large sized macrofaunal taxa were also observed during the development of restored young Indo-Pacific mangroves [97, [100][101][102].
As the FG coast is part of a unique sedimentary system influenced by the tremendous volumes of mud discharged from the Amazon River, the comparisons with mangroves from other biogeographical regions facing less dynamic environmental conditions remains difficult.Moreover, very few studies of benthic infauna have been carried out on the Guianas coast, and those studies concerned the bare intertidal mud flats but not the mangrove covered upper part of the mudbanks.This underscores the difficulty Alongi [20] encountered when he attempted a comparison between Amazon influenced mangroves and Australian bare tidal flats and mangroves facing the Great Barrier Reef.The main differences in infaunal abundances between these biogeographical regions come from the very great abundances that we measured in FG.At similar mesh sizes and sampling area units, the total densities of mesofauna (0.25 mm < x < 1 mm) were 18-fold higher compared to those of tropical tidal flats in Australia, Thailand and Chile [21].Considering the similar ages of mangroves and only macrofaunal densities (>1 mm), the latter were at least 11-fold higher than in areas of restored pioneer and young mangroves [97, 100,102] and in naturally developing mangroves in Malaysia and China [16].Although we know that the vegetation structure impacts benthic composition, we contend that the persistently high and reoccurring sediment inputs and the exceptionally rapid environmental changes associated with mangrove development promote the unique characteristics of the infaunal benthic succession in FG.The exceptionally high abundances have not been observed in the mangroves east of the Amazon delta, an area of the coast of Brazil not influenced by the Amazon system, and depleted of mud [39].
The diversity indices and the OTU richness that we measured at the vegetated stages are similar to those of other mangrove benthic assemblages not influenced by the Amazon muds [16,39,100].Clearly, benthic species diversity in FG mangroves is not impoverished compared with the other regions, despite the regional sediment dynamics, rapid mangrove growth, and monospecific mangrove forests.However, further studies are needed on the benthic diversity within muds associated with older stages of mangroves in FG in order to verify this hypothesis.Indeed, sediments associated with mature mangroves are generally characterized by harsh conditions (long substrate air exposure times, sediment compaction and drought, dense root networks and low redox conditions) potentially reducing the infaunal density, and stabilizing or decreasing the OTUs richness [16,103].
The unique "mangrove belt" ecosystem along the coasts of the Guianas in South America is the result of sediment accretion and erosion phases and model ecological resilience [38].The benthic diversity of FG mangroves follows the succession of repetitive erosion/accretion phases.In his renowned study on the resilience and stability of ecological systems, Holling [104] highlighted the importance of knowing the amplitude and frequency of oscillations to understand the behavior of particular systems.Whether benthic infaunal assemblages tend toward a persistent structure or permanently adjust to this constantly evolving environment [105] remains to be determined, given the long term oscillation of erosion/accretion phases along the FG coast [25].Long term temporal monitoring of growing mangroves (including the oldest stands, and seasonal variability with potential effects of rain and freshwater inputs) is thus required to validate the hypothesis of permanent benthic adaptability to constantly varying environmental conditions in these coastal mangrove systems.

Conclusions
This study is the first to quantify variance partitioning (i.e., beta diversity) in benthic community structure in response to mangrove development.We showed that, on the Amazon River influenced coast of the Guianas, during the dry season, four sets of variables can explain the spatial structure of the benthic infaunal communities.The single most important one was a linear spatial structure associated with changes in substrate elevation, porewater salinity, TC content and pneumatophore density, which explained 20-30% of the infaunal structure variations along a west to east gradient of mangrove maturation.Up to 7-12% of the variance partitioning was attributable to nonlinear spatial patterns and sediment organic matter quality changes between the three main mangrove stages (bare sediments, pioneer and young mangroves).However, 7% to 9% of the beta diversity was not explainable by the variables used in our analysis, and, rather, could be the result of other biotic or abiotic variables, biological interactions or neutral processes, likely related to small scale spatial heterogeneity (i.e., substrate geomorphology, pneumatophores and crab burrows), but not described here.This study has highlighted initial hypotheses that are relevant to initiatives on biodiversity evaluation or the success of mangrove restoration through monitoring of benthic biodiversity.Variables do not change only linearly with mangrove growth, but they also depend strongly on microhabitat (along with sediment parameters).Scientists and managers involved in such initiatives are urged to identify the spatial heterogeneity of their marine area and should identify the kind of microhabitats they could include in their monitoring designs.
By examining two infaunal size structures, we demonstrated that the mesofaunal communities (0.25 mm) were more heterogeneous and diverse than the macrofaunal ones (>1 mm) in the bare sediments.On the contrary, the macrofaunal communities became heterogeneous and diversified in the vegetated sediments.Overall, mangrove early development increased the alpha and beta diversities of the infaunal communities for the two size classes.
This study raises questions regarding the meso-and macro infaunal benthic community composition and structure within older stages of mangroves.Whether this community structure is repeated at each new mangrove cycle or is permanently modified along the FG coast, depending on local environmental conditions and seasonal variables, remains to be determined.In a global context of threats to mangroves worldwide, future studies in FG should consider additional sampling of young and older mangrove stands over both the wet and dry seasons, in order to be able to generalize patterns in infaunal succession associated with the entire mangrove community life cycle.

Supplementary Materials:
The following are available online at https://www.mdpi.com/article/10.3390/su14052829/s1. Figure S1: Microphytobenthos (diatoms and cyanobacteria) (a) mean abundances (Cells mL −1 ) and (b) mean relative density (%) along the mangrove maturation gradient (transects B-I; mean ± SD, 2 < n < 6). Figure S2: Relationships between mangrove stages and surface sediment (black symbols) and 18 centimeter-depth averaged (grey dots) parameters at transects A, E and I. Empty and filled dots denote replicates and means, respectively.The significance values of the ANCOVAs are shown on the graphs.ANCOVA was not applied to salinity (no core replicates per transect).Figure S3: (A) mean relative density (n ≤ 9; %) of the meiofauna OTUs at mobile mud stations (AA, J, K and L) and at some consolidated stations of the mangrove maturation gradient (A, B, E, F, G, H and I).(B) mean abundances of tanaids (No.g −1 DW) per size class (1 to 6 mm) at the stations (A, B, J, K) containing tanaids (mean ± SD, n ≤ 6). Figure S4: Bubble plot of predicted position (linear constraints) on the only canonical axis for the linear models of Hellinger transformed abundances as a function of (A) east-west and (B) south-north coordinates.(C) Bubble plots depicting the selected positive dbMEM spatial Eigen functions after forward selection.The dbMEM1, dbMEM5 and dbMEM6 were included in the models for detrended Hellinger transformed abundances and biovolumes.The dbMEM7 was only included in the model of abundances.Figure S5: RDA correlation triplot of infaunal community data and selected environmental variables based on Hellinger transformed (A) abundance and (B) biovolume infaunal data.Stations of the three groups of transects are represented by specific colored circles (white: A-C, group 1; grey: D-G, group 2; black: H-I, group 3)."Grain size" refers to median grain size, "Pneumatophore" refers to roots density, "Elevation" refers to substrate height.In the top graph "Gr" denotes three taxa clustered too close to be distinguished Gastropoda sp.1, Bivalvia sp.1 and Melitidae sp.1.Table S1: Taxonomic list of the OTUs (nd: not defined).a little over 30% of variation in biovolume (9999 perm., R 2 Adj = 0.33, F = 3.590, p < 0.001).Forward selection of variables resulted in a model with C:N ratio, pneumatophore density and Phaeo:Chl-a ratio for abundances (9999 perm., R 2 Adj = 0.49, F = 9.230, p < 0.001), whereas only C:N ratio and pneumatophore density were selected for biovolumes (9999 perm., R 2 Adj = 0.29, F = 6.401, p < 0.001).The selected biotic and abiotic environmental variables were retained for variation partitioning analysis.

Figure 2 .
Figure 2. Study design employed for the meso and macro benthic community structure and surfac sediment characteristics and their relationships along a mangrove growth gradient (0 to 2 years old referring to the transects from A to I).

Figure 2 .
Figure 2. Study design employed for the meso and macro benthic community structure and surface sediment characteristics and their relationships along a mangrove growth gradient (0 to 2 years old, referring to the transects from A to I).

Figure 3 .
Figure 3. Principal component analysis performed on the correlation matrix of environmental variables (tree biomass, density and height, pneumatophore density, litter biomass, grain size, water content (WC), porosity, salinity, total carbon (TC), total nitrogen (TN), Chla and phaeopigments, TC:TN, Chla:phaeo and C-Chla:TOC ratios, along the mangrove maturation gradient (0 to 2 years old, referring to the transects from A to I).Together, axes 1 and 2 explain 68.7% of total variability in station data.The larger dots illustrate the centroids of the group of individuals for each transect.

Figure 3 .
Figure 3. Principal component analysis performed on the correlation matrix of environmental variables (tree biomass, density and height, pneumatophore density, litter biomass, grain size, water content (WC), porosity, salinity, total carbon (TC), total nitrogen (TN), Chla and phaeopigments, TC:TN, Chla:phaeo and C-Chla:TOC ratios, along the mangrove maturation gradient (0 to 2 years old, referring to the transects from A to I).Together, axes 1 and 2 explain 68.7% of total variability in station data.The larger dots illustrate the centroids of the group of individuals for each transect.

Figure 4 .
Figure 4. (a) Bed elevation (m) along the study grid area.Pink dots indicate the locations of the topographic measurements; interpolations between dots were carried out following the natural neighbor method.Letters and numbers refer to the sampling stations of the grid.(b-d) Vegetation structure along the study grid: (b) tree biomass (tons of dry matter per ha −1 ), (c) tree density (individuals per ha −1 ), (d) tree height (m); values averaged at each station and linearly extrapolated between stations along the mangrove gradient.

Figure 4 .
Figure 4. (a) Bed elevation (m) along the study grid area.Pink dots indicate the locations of the topographic measurements; interpolations between dots were carried out following the natural neighbor method.Letters and numbers refer to the sampling stations of the grid.(b-d) Vegetation structure along the study grid: (b) tree biomass (tons of dry matter per ha −1 ), (c) tree density (individuals per ha −1 ), (d) tree height (m); values averaged at each station and linearly extrapolated between stations along the mangrove gradient.

Figure 7 .
Figure 7. Principal component analysis based on Hellingher transformed data for OTU abundances (a) and biovolumes (b) along the mangrove maturation gradient (transects A to I).For the abundance (a) and biovolume (b), the two axes explained 64.41% and 49.15% of total variance, respectively.Plotting characters for sites correspond to K-means clustering into 3 groups (circles, triangles, squares) with shades of grey corresponding to the established groups of transects based on the environmental PCA (white circles: A-C (group 1); grey boxes: D-G (group 2) ; black triangles: H-I (group 3).

Figure 7 .
Figure 7. Principal component analysis based on Hellingher transformed data for OTU abundances (a) and biovolumes (b) along the mangrove maturation gradient (transects A to I).For the abundance (a) and biovolume (b), the two axes explained 64.41% and 49.15% of total variance, respectively.Plotting characters for sites correspond to K-means clustering into 3 groups (circles, triangles, squares) with shades of grey corresponding to the established groups of transects based on the environmental PCA (white circles: A-C (group 1); grey boxes: D-G (group 2); black triangles: H-I (group 3).

Figure 8 .
Figure 8. Variation partitioning of the Hellinger transformed abundance (a) and biovolume (b) data among the four sets of explanatory variables: selected abiotic environmental variables (i.e., elevation, grain size, salinity), cartesian spatial coordinates (linear spatial structure); selected dbMEM spatial eigenfunctions (non linear spatial) and selected biotic environmental variables (i.e., phae:chla, TC:TN, pneumatophores).Values are R 2 Adj; values < 0.005% are not shown.The underlined values denote the variable or convergence of variables that mostly explain the infaunal community structure.

Figure 8 .
Figure 8. Variation partitioning of the Hellinger transformed abundance (a) and biovolume (b) data among the four sets of explanatory variables: selected abiotic environmental variables (i.e., elevation, grain size, salinity), cartesian spatial coordinates (linear spatial structure); selected dbMEM spatial eigenfunctions (non linear spatial) and selected biotic environmental variables (i.e., phae:chl-a, TC:TN, pneumatophores).Values are R 2 Adj ; values < 0.005% are not shown.The underlined values denote the variable or convergence of variables that mostly explain the infaunal community structure.