Meiofauna at a tropical sandy beach in the SW Atlantic: the influence of seasonality on diversity

Study area and sampling

The study was carried out at Gramuté, a sandy beach located within a marine protected area in the Eastern Brazilian Marine Ecoregion (Fig. 1A). The region is geomorphically characterized by abrasion terraces of Barreiras Formation from the coast to the inner continental shelf (Martin et al., 1996). Gramuté beach is marked by scattered intertidal lateritic reefs (Mazzuco, Stelzer & Bernardino, 2020), with the presence of carbonate secreting organism which contributes to the deposition of bioclastic sediment (Albino & Suguio, 2011). Additionally, Gramuté beach is marked by strong internal tidal currents, and E-SE wave swells with upwelling events occurring mostly during spring and summer (Pereira et al., 2005). This tropical region is marked by dry winters and rainy summers (Bernardino et al., 2015), with sea surface temperatures ranging between 21 °C and 27 °C, and salinity ranging from 34.6 to 36 ppt (Quintana et al., 2015; Mazzuco et al., 2019; Mazzuco, Stelzer & Bernardino, 2020). This region has also experienced significant warming in the last 40 years (Bernardino et al., 2015; Mazzuco, Stelzer & Bernardino, 2020).

The study region is marked by frequent exposure to waves generated mainly by the South Atlantic Subtropical Anticyclone (ASAS), with northeast (NE) swells mainly. Although there is dominance of NE waves throughout the year, in the autumn and winter period the wind regime changes to E-SE, strengthening the presence of waves from these directions (E-SE), with average heights of 1.5 m. During winter, the region is also affected by the passage of frontal systems, making it susceptible to wave action coming from the south-southwest (S-SW) (Silva, Hansen & Cavalheiro, 1984).

During a year (December 2019 to November 2020), we monthly collected sediment samples (approximately 200 g each replicate) on three different stations (n = 9 sediment samples per month) at Gramuté beach on the low-tide shoreline, spaced 20 m apart (Fig. 1B). Sediment samples were collected manually using sterile, DNA-free corers, over all seasons during the sampling period (Table S1). Additionally, we collected samples for sediment analysis (grain size, total organic matter, carbonate content, and sedimentary organic biopolymers). All samples were transported in thermic bags with ice, and stored at −20 °C until analysis. Field sampling was authorized by the Biodiversity Authorization and Information System of the Brazilian Institute for the Environment and Renewable Natural Resources (SISBIO-IBAMA, sampling license number 24700-1). Total monthly rainfall data for the sampling period (December 2019—November 2020) was obtained from the National Water Resources Information System (SNIRH) portal, made available by the National Water and Sanitation Agency (Agência Nacional de Águas (ANA) (2023); ANA—https://www.snirh.gov.br/hidroweb/), considering the station of Santa Cruz -Litoral (code: 1940002; Lat:−19.9578, Lon:−40.1544), which is approximately 4 Km from Gramuté beach.

Sediment analysis

Sediment samples were dried at 60 °C for 48 h before all granulometric analysis. Dried sediment was macerated and sieved in mesh openings of 63 µm to 2 mm in a sieve shaker to determinate the carbonate content by muffle combustion at 550 °C for 4 h with an additional hour at 800 °C. Additionally, we quantified total organic matter (TOM) by weight loss after combustion (500 °C for 3 h) (Suguio, 1973).

Sedimentary organic biopolymers (proteins, carbohydrates, and lipids) were analyzed following Danovaro (2010). After extraction with NaOH 0.5 M we determinated total protein (PRT) content according to Hartree (1972) as modified by Rice (1982) to compensate for phenol interference. For total carbohydrates (CHO) analysis, we followed the protocol from Gerchacov & Hatcher (1972). Total lipids (LIP) were extracted from 1 g of homogenized sediment lyophilized by ultrasonication in 10 ml of chloroform: methanol (2:0 1 v/v) and analyzed according to Marsh & Weinstein (1966). The concentrations of PRT, CHO, and LIP are presented respectively as bovine serum albumin, glucose, and tripalmitin equivalents. Concentrations of PRT, CHO, and LIP were converter to carbon equivalents following Fabiano & Danovaro (1994) using conversion factors of 0.49, 0.40, and 0.75, respectively. The sum of PRT, CHO, and LIP carbon equivalents are presented as biopolymeric carbon (BPC) (Fabiano, Danovaro & Fraschetti, 1995). Further, protein to carbohydrate (PRT: CHO) and carbohydrate to lipid (CHO: LIP) ratios were used to assess biochemical degradation processes (Galois et al., 2000). All analyzes were performed in triplicate and blanks were carried out for all analysis with pre-combusted sediments at 450 °C and 480 °C for 4 h.

Seascapes identification

The Marine Biodiversity Observation Network (MBON) Seascapes are a characterization of water masses with particular biogeochemical features obtained from satellite and modeled data that comprises different oceanic parameters (sea surface temperature—SST, sea surface salinity—SSS, absolute dynamic topography—ADT, chromophoric dissolved organic material—CDOM, surface chlorophyll-a–Chl-a, and normalized fluorescent line height—NFLH). These variables are used for a categorization system of 33 water masses that represents different marine scenarios/conditions (Montes et al., 2020).

Oceanographic conditions were characterized according to the variation in MBON Seascape Pelagic Habitats Classification (Kavanaugh et al., 2014, 2016; Mazzuco & Bernardino, 2022) using the database available in the NOAA Coast and Ocean Watch Programs, with monthly frequency on a 5 Km² grid (Kavanaugh et al., 2014, 2016), to characterize the seascapes for the Área de Proteção Ambiental Costa das Algas (~30 km coastline, 465 Km², Longitude–40.3° to–39.8°, Latitude 20.3° to 19.8°) for the study period (December 2019–November 2020). Additionally, to determine seasonal SST and SSS for the study area, we calculated a weighted average based on the monthly coverage area of each identified MBON marine seascape.

MBON Seascapes are presented as seascape coverage (%), which represents the extent of an area that is encompassed within any of the MBON Marine Seascapes categories. This percentage represents how much of the area of the Área de Proteção Ambiental Costa das Algas (which is a Marine Protected Area) is encompassed with one of the seascape categories. Each MBON Seascape category is defined by a fixed value for each oceanic variable, and the seascape product is generated as monthly and 8-day composites at 5 Km spatial resolution. In this manuscript we used the seasonal mean, calculated as the mean of the monthly seascape coverage for all 3 months per season.

DNA extraction and sequencing

Prior to DNA extraction, sediment samples were elutriated using sieves of 45 μm mash in an attempt to increase the meiofaunal abundance and enrich metazoan DNA recovery as suggested by Brannock & Halanych (2015).

A total of 1 L flasks were filled with 950 mL of filtered seawater and sediment samples were added to it, then homogenized and let to sit (for approximately 30 s) before decanting the liquid over the 45 μm sieve. This procedure was repeated 10 times for each sediment sample. The sediment retained on the sieve was transferred to 50 mL Falcon tubes and centrifuged at room temperature for 3 min at 1,342 × g in an Eppendorf Centrifuge 5430, then the sample volumes were standardized to 20 mL. We mixed the samples in the Falcon tubes, then aliquots of 1 mL (two tubes per sample) were stored at −20 °C (Brannock & Halanych, 2015). All glassware was cleaned and autoclaved between samples to avoid cross contamination between samples, and sieves were sterilized by soaking for 45 min in 10% sodium metabisulfite solution (Creer et al., 2010; Brannock & Halanych, 2015). DNA was extracted from 1 mL aliquots using the PowerSoil DNA® (Qiagen) kit following the manufacturer’s instructions. DNA integrity was verified in 1% agarose gel, and purity using NanoDrop One spectrophotometer (Thermo Fisher Scientific Inc., Waltham, MA, USA). We measured DNA concentration using the Qubit® 4 Fluorometer (Qubit™ 1X dsDNA HS Assay Kit—Life Technologies-Invitrogen, Carlsbad, CA, USA). Blank samples were carried in triplicate for each step before metabarcoding sequencing (sediment elutriation, DNA extraction, and integrity, purity, and concentration checking).

DNA samples extracted from the same sediment sampling station, for each month separately, were combined into a single pool, totaling three samples per month (nine samples per season, totaling 36 samples). PCR, library preparation, and sequencing were conducted by ©NGS Genomic Solutions (Piracicaba, SP, Brazil). Metabarcoding sequencing was performed using the MiSeq Illumina platform (2 × 250 bp, with a coverage of 100,000 paired-end reads per sample), sequencing the V9 hypervariable region from 18S SSU rRNA gene using the primers Euk_1391 forward (GTACACACCGCCCGTC) and EukBr reverse (TGATCCTTCTGCAGGTTCACCTAC) (Medlin et al., 1988; Amaral-Zettler et al., 2008; Stoeck et al., 2010), generating amplicons that could vary in size (mean 260 ± 50 bp), once the reverse primer doesn’t have its exact position conserved as the forward one.

Bioinformatic pipeline

Bioinformatic analysis were conducted using an AMD Ryzen 1950× Crucial 64 GB (16 × 4) DDR4 2,666 MHz computer. We used the QIIME2 2022.8 software to identify sequences with the demultiplexed raw paired-end reads (Bolyen et al., 2019). Firstly, we imported FastQ files as QIIME2 artifacts, then denoised them via DADA2 (Callahan et al., 2016) using the denoise-paired plugin, and removed low-quality bases and primer sequences.

The taxonomic assignment of amplicon sequence variants (ASV) generated by the DADA2 plugin (p-trim = 10, p-trunc = 160, and mean phred score = 39 ± 1; Table S2) was performed using the machine learning Python library scikit-learn (Pedregosa et al., 2011). A pre-trained Naïve Bayes classifier trained on the Silva 138 database (Quast et al., 2013) clustered at 99% similarity was used to identify taxonomically the DNA sequences. Due to the differences on the number of identified sequences, the dataset was normalized to allow analysis and comparisons under equal sampling depth. We used the spring dataset minimum sampling depth (1,384 reads) and resampled each sample to the same depth. This normalized dataset was used to calculate all diversity metrics. We performed rarefaction curves for all four sampled seasons (summer, winter, spring, and autumn) with the ASVs. We calculated the Faith’s Phylogenetic Diversity (PD) for each sample using the diversity core-metrics-phylogenetic pipeline from QIIME2. The PD was calculated based on phylogenetic trees generated using the align-to-tree-mafft-fasttree pipeline from the q2-phylogeny plugin from QIIME2. Shannon diversity was calculated using the qiime diversity alpha pipeline and setting the p-metric parameter to “Shannon”. Raw sequences are available online on NCBI (SRR24675047) and on Brazilian Biodiversity Information System (SiBBr; Bernardino & Coppo, 2024).

Statistical analysis

In our statistical analyses, only meiofaunal metazoan sequences used. The dataset was subset to contain only metazoan sequences based on taxonomic annotations using the qiime taxa filter-table function from the QIIME2 q2-taxa plugin. Here we considered all the exclusively meiofaunal phyla (Gnathostomulida, Kinorhyncha, Loricifera, Gastrotricha, and Tardigrada) and other metazoans that can be meiofaunal-sized during their life cycle (temporary meiofaunal taxa) (Higgins & Thiel, 1988; Giere, 2009). After this filtering, the final dataset contained 10 phyla (Annelida, Cnidaria, Crustacea, Echinodermata, Gastrotricha, Mollusca, Nematoda, Nemertea, Platyhleminthes, and Rotifera) as previously implemented in other studies (Brannock & Halanych, 2015; Bernardino et al., 2019; Fais et al., 2020; Bellisario et al., 2021; Castro et al., 2021; Coppo et al., 2023). Permutational analysis of variance (PERMANOVA; Anderson, Gorley & Clarke, 2008) was performed to compare environmental variables (rainfall, temperature, salinity, carbonate content, grain size, total organic matter, and biopolymeric composition), seascape coverage (the extent of an area that is encompassed within any of the MBON Marine Seascapes categories), and meiofaunal data (diversity metrics—Shannnon’s diversity index, phylogenetic diversity, and abundance of sequence reads) among sampled seasons (summer, autumn, winter, and spring) and sampled stations at Gramuté beach. Bonferroni corrections was used to adjust p-values for pairwise comparisons (Bonferroni, 1936). A canonical analysis of principal coordinates (CAP; Anderson & Willis, 2003) was performed with environmental variables (rainfall and sediment variables) and the meiofaunal assemblage composition at Phylum level (square-root transformed). Additionally, a similarity percentage routine (SIMPER; Clarke, 1993) was applied to define the taxa that most contributed to the dissimilarity among seasons, based on a Bray-Curtis dissimilarity matrix. A multiple linear regression was fit using Shannon’s diversity Index and phylogenetic diversity as response variables, and the assessed environmental variables as predictive variables. After testing for multicollinearity among variables, we removed carbohydrate content (CHO) which was highly correlated to carbohydrate-to-lipids ratio (CHO:LIP) and biopolymeric carbon (BPC). Normality tests were run on model’s residuals through QQ-plots and Shapiro-Wilk normality tests. After obtaining the multiple linear regression values, we used the Akaike Information Criterion (AIC), through a stepwise backward model configuration and the final model was chosen based on the lowest AIC value (Akaike, 1978). Significative differences were defined at p < 0.05. All graphical and analytical procedures were performed in the R environment (R Core Team, 2022).

Environmental conditions and seascape coverage

Significant seasonal variability was observed at Gramuté beach during the period studied (PERMANOVA, df = 3; Pseudo-F = 6.916; p = 0.001; Table S3), with higher LIP content on autumn (Fig. 2; Table S3) and lower on spring (Fig. 2; Table S3), meanwhile the CHO:LIP ratio was higher on winter than on other seasons (Fig. 2; Table S3). Total rainfall ranged from 80.2 ± 38.4 mm in summer to 193.0 ± 42.2 mm in autumn (Fig. 2). The sediment is completely composed of sand, mainly by medium and coarse sand, with carbonate content ranging from 26 ± 10% during spring to 55 ± 10% in winter (Fig. 2). Total organic matter (TOM) had its lower concentration in summer (8.6 ± 3.8%), and higher in spring (10.4 ± 7.9%) (Fig. 2; Table S3). The protein fraction of the organic matter content in the sediment ranged from 48.2 ± 21.6 mg/g in autumn to 96.9 ± 14.8 mg/g in summer (Fig. 2; Table S3), while the carbohydrate fraction ranged from 997.1 ± 193.5 mg/g in autumn to 2,102.0 ± 1,435.0 mg/g in winter (Fig. 2; Table S3). The labile fraction of the organic matter, which is represented by the BPC, was similar among seasons (Fig. 2; Table S3).

Overall, the Seascapes categories in this region were characterized by high sea surface temperature (SST > 20.9 °C), high sea surface salinity (SSS > 33.6 psu) and calm waters (absolute dynamic topography-ADT ranging from 0.51 to 0.83 m). The seascapes had wide ranges in dissolved organic matter (CDOM; 0.00 to 0.07 m⁻¹), chlorophyll-a concentration (CHLA; 0.07 to 2.09 mg.m⁻³), and fluorescence (NFLH; 0.02 to 0.24 W.m⁻².um⁻²sr⁻¹) (Fig. 3). We observed changes in the frequency of seascapes in the studied area along the year (PERMANOVA, df = 3; Pseudo-F = 8.014; p = 0.001; Table 1). Seascapes Tropical Seas (class 15–38.4% of area coverage during sampling period), Subtropical Gyre Transition (class 5–19.0% of area coverage during sampling period), Subtropical Gyre Mesoscale Influenced (class 13–18.3% of area coverage during sampling period), and Warm, Blooms, High Nutrients (class 21–12.4% of area coverage during sampling period) were the most frequent, with more than 80% of area coverage during the study period (Fig. 3).

Water masses at Gramuté beach during summer (Dec–Feb), autumn (Mar–May) and winter (Jun–Aug) were dominated by the Seascape Tropical Seas (class 15), which is characterized by high temperatures (25.4 °C) and salinity (35.4 psu), and covered with 40.9% (summer), 43.1% (autumn), and 45.1% (winter) of the study area (Fig. 3). During spring, the dominance of seascapes at Gramuté changed due to an intrusion of subtropical a water mass (class 13–42.7% of area coverage; Fig. 3), characterized by significant lower temperature (23.5 °C) and higher salinity (35.9 psu).

Meiofaunal assemblage

A total of 9,692 sequences from meiofaunal taxa were identified in the dataset. We did not observe significative differences in meiofauna composition among sampled stations (PERMANOVA, df = 2; Pseudo-F = 0.963; p = 0.491; Table 2). Nonetheless, there were significant seasonal variations, with winter differing from all other seasons (PERMANOVA, df = 3; Pseudo-F = 2.307; p = 0.002; Table 2; Fig. 4). During the summer (35% and 40% of reads), autumn (43% and 34% of reads), and spring (59% and 27% of reads), Crustacea and Annelida were the most prevalent taxa. The two taxa with the highest abundance throughout the winter were Crustacea (57% of reads) and Nematoda (17% of reads) (Fig. 4). Nemertea was not detected during autumn, Gastrotricha was not detected in spring, and Rotifera was not detected in neither. Only 11 taxa (e.g., Harpacticoida, Podocopida, and Chromadorea) were detected on all seasons (Table S4), meanwhile 14 taxa were detected only on one sampled season (Table S4).

Fewer meiofaunal taxa were detected in spring than all other seasons (Fig. 5), influencing on significant seasonal differences on diversity patterns in Gramuté beach. Lower phylogenetic diversity was registered in spring (9.23 ± 1.88) and in autumn (11.88 ± 1.82) than in summer (17.93 ± 3.11) and winter (19.37 ± 4.85) (PERMANOVA; Pseudo-F = 18.863; df = 3; p < 0.001; Table 3). Similarly, Shannon’s diversity was 1.7-fold, 2.0-fold, and 1.9-fold lower in spring than in autumn, summer, and winter, respectively (PERMANOVA; Pseudo-F = 13.129; df = 3; p < 0.001; Table 4).

Meiofaunal assemblages differed significantly among the seasons in Gramuté beach (PERMANOVA, df = 3; Pseudo-F = 2.353; p = 0.001; Table 2; Table S4). Dissimilarity levels ranged from 49.7% (between winter and summer) to 68.6% (between autumn and summer). SIMPER analysis revealed that Annelida (ranging from 16.5% to 28.3%; Table S5), Crustacea (ranging from 21.8% to 26.7%; Table S5) and Nematoda (ranging from 13.9% to 21.8%; Table S5) were the taxa that most contributed to the differences among all seasons. Platyhelminthes contributed 15.4% to the total dissimilarity of 49.5% between autumn and spring (Table S5). Annelids, crustaceans, and nematodes were more abundant in summer and winter. Furthermore, these taxa were highly associated with higher organic matter content and quality (total organic matter content, biopolymeric carbon, protein content, and protein-to-carbohydrate ratio; Fig. 6).

Protein content (PRT), lipid content (LIP), biopolymeric carbon (BPC), and carbohydrate-to-lipids ratio (CHO:LIP) composed a significant model of variables likely to drive meiofauna diversity (Shannon’s diversity Index and phylogenetic diversity) at Gramuté beach (Adjusted R² = 0.602; F = 7.319; p < 0.001; Table 5). We observe significant positive relationship between meiofaunal diversity (Shannon’s diversity and phylogenetic diversity) and lipid content (LIP) (t = 2.513; p = 0.018; Table 5). Biopolymeric carbon content—BPC showed a significant negative relationship to diversity (t = −2.584; p = 0.015; Table 5), as well as protein content (PRT), although it was not significative (t = −1.719; p = 0.096; Table 5).