Reef Larval Recruitment in Response to Seascape Dynamics in the SW Atlantic

Advances in satellite observation have improved our capacity to track changes in the ocean and seascapes with numerous ecological and conservation applications, but yet under explored for coastal ecology. In this study, we assessed dynamics in the Seascape Pelagic Habitat Classification, a satellite remote-sensing product developed by NOAA to monitor biodiversity globally, and invertebrate larval recruitment in order to identify and predict changes in coastal benthic assemblages at tropical reefs in the SW Atlantic. Our results revealed that pelagic Seascapes correlated with monthly and seasonal variations in recruitment rates and assemblage composition. Recruitment was strongly influenced by subtropical Seascapes and was reduced during warm, blooms, and high-nutrient waters, likely to affect reef communities in the long term. Modeling indicate that Seascapes may be more efficient than temperature in predicting benthic larval dynamics. Based on historical Seascape patterns, we identified seven events that may have impacted benthic recruitment in this region in the last decades, which not surprisingly, coincided with consistent global heatwaves. These findings provide new insights into the application of novel satellite remote-sensing Seascape categorizations in benthic ecology and evidenced how reef larval supply in the SW Atlantic could be impacted by recent and future ocean changes. 17 18 21 23


Introduction
Climate change has caused severe impacts to marine ecosystems during the last decades. Identifying the drivers and ecological responses from warming oceans and associated effects is imperative, and long-term monitoring of marine ecosystems will be a key step to this process. Global warming is changing the pelagic ocean Seascapes through altered temperature and nutrients, which is resulting in changes of biological and ecological processes in many marine species from microorganisms to top predators [1,2,3,4]. Ocean warming is among the major causes of altered ecological dynamics in marine ecosystems, shifting patterns of species distribution and increasing species loss at unprecedented rates [5,6,7]. The tropicalization of temperate and subtropical oceans is already ongoing while tropical sites are becoming inhabitable for several species, a scenario that may become worse in the near future [8,9,10,11]. In this changing ocean, there is an urgent need to anticipate biodiversity patterns and understand how ecosystem processes will behave across climatic conditions and altered seascapes. This approach can improve our capacity to predict biodiversity loss and recovery, and thus, provide valuable information to maximize management and mitigate climate impacts with more specificity [12,13].
As many marine organisms have a life cycle with an intermediate larval stage, local biodiversity of marine ecosystems is strongly linked to oceanographic and ecological processes that influence the larval supply and recruitment (e.g., [14,15]). Processes influencing benthic recruitment may be acting at distinct spatial and temporal scales related to pelagic conditions that affect planktonic larval development and supply rates [16,17,18]. Thus, lower numbers of benthic recruits over an area may be a response to changes in seascapes, which will have later consequences to the overall biodiversity of that region [19,20]. For these benthic species, recruitment also requires surviving strong local post-settlement pressures, such as competition, predation, and unfavorable environmental conditions that impact juvenile abundances [21,22]. The reduction of local populations is the first step towards species extirpation and an alarm to potential biodiversity loss at larger scales [23,24]. Although vulnerable populations are commonly identified through adult mortality, reduction of offspring supply can be used as an alert for potential species disappearance and change in the health status of populations and communities (e.g. [25,26]).
Monitoring larval dynamics can potentially anticipate ecological responses to changes in marine seascapes, but it is a challenging task and often neglected by long-term observing programs. For rocky shore and coral reef populations, which are widely studied marine ecosystems with respect to larval dynamics, variations in larval recruitment may be estimated using pelagic variables such as sea surface temperature, chlorophyll-a concentrations, current velocities, wave height, and the occurrence of oceanographic features such as upwellings, cold fronts, and storms [27,28,29,30,31]. Several studies indicate that the early life stages of benthic invertebrates are highly vulnerable to a changing ocean, where even small changes in pelagic primary productivity or temperature fluctuations may have strong impacts on benthic recruitment (e.g. [32]). Modeling larval dynamics from pelagic variables could be a useful approach to develop ecological metrics of benthic community change and, in turn, to monitor the coastal ecosystems.
Currently, oceanographic variables are readily available online in a wide range of spatial and temporal resolutions that are derived from in situ monitoring networks or satellite products [33,34].
These ocean observing databases are rich, diverse, and can be easily retrieved to develop panoramas of the pelagic systems at multiple levels and for different purposes. Among them, the Seascapes Pelagic Habitat Classification developed by NOAA is a novel satellite remote sensing approach that translates spatial and temporal variations in seascapes from multiple change in pelagic conditions (CoastWatch/OceanWatch) [35,36]. The NOAA's Seascape classification is obtained from dynamic fields of satellite and modeled data designed to support biogeographical assessments and marine observation networks, providing an environmental tool to track changes in marine habitats with relatively high resolution. By integrating key pelagic conditions and state, this Seascape classification goes beyond single-variable categorizations (temperature anomaly; e.g., [37]) by reflecting the quality and extent of different oceanographic domains or features. Therefore, global ocean classification is a promising tool to assess and predict pelagic-driven biological processes (phytoplankton, e.g., [38]).
Considering that larval dynamics is key to the ecology of benthic communities and strongly dependent on oceanographic conditions, we examined whether changes in integrated Seascape Pelagic Habitat Classification (hereafter Seascapes) can be associated with larval recruitment at tropical intertidal reefs and could be applied to forecast recruitment variability. The study was carried at a Brazilian LTER site (HCES) on the coastal SW Atlantic in a marine protected area located in the Eastern Brazil Marine Ecoregion. This region has experienced significant warming trends during the last four decades (1970-2010; [39]) and previous work indicated that larval recruitment and adult intertidal reef assemblages in this region are vulnerable to periods of high sea temperatures [40,41].
There has been a significant advance in long-term monitoring of the coastal ecosystems in the SW Atlantic in the last decades, particularly in for the tropical reefs [42,43,44]. However, the scarcity of basic knowledge about reef dynamics in particular regions and punctual ecological modeling efforts imposes several challenges to reef management, such as preventing the anticipation of potentially harmful conditions and forecasting ecosystem changes.
Using the LTER HCES dataset from that represents regional recruitment trends [40,41], we quantified recruitment across multiple temporal scales (years, seasons, and months) for several reef invertebrate species. We then tested the association between recruitment patterns (abundance and composition of species) and Seascapes, identifying the Seascape classes that were key to recruitment variability. We then evaluated the effectiveness of Seascapes in predicting recruitment success (i.e. higher or lower abundance and diversity), when compared to a single-variable approach using sea surface temperature as a driver of larval dynamics in this region [40]. We also tested if past Seascape dynamics could indicate periods of decreased recruitment potential by,revisiting the previous 16 years (2003 to 2019) of monthly Seascape conditions in the area of study. Our hypothesis was that recruitment rates (i.e. abundance) and taxonomic composition would be highly correlated to changes in Seascapes through time, reflecting both local and regional variability in pelagic habitats. This study provides new insights into the application of pelagic Seascape classifications to the supply-side ecology at tropical reefs in the tropical South Atlantic and reveals how these benthic communities may be highly vulnerable to a warmer ocean. We believe that identifying environmental and connectivity drivers for marine assemblages in these reef ecosystems is the first step towards improving ecological predictions and could be an important subside for planning and mitigation of environmental impacts at SW Atlantic reefs and coastal MPAs [45,46,47].

Seascapes
We identified six Seascape classes within the study region, including the Tropical-Subtropical  The Seascapes Warm, Blooms, High Nutrients, and Tropical conditions were dominant along the studied period with over 70% of coverage during the study period (Fig. 1). As expected, the predominance of Seascape classes differed seasonally and monthly (seasons F = 7.64, p = 0.01, month F = 1.82, p = 0.04, Table 1; p < 0.05, See Supplementary Table S3). Subtropical conditions (Seascape classes 3, 5, and 13) were more frequent during the fall (Mar-May) when the pelagic conditions were similar and more influenced by southern Atlantic water masses (Fig. 1). On the other hand, Hypersaline Eutrophic conditions occurred frequently during summer months (Dec-Fev) at both regional and mesoscale spatial scales, with the predominant influence of tropical warm currents ( Fig. 1). The seasonal variability on Seascape conditions did not change at the annual scale (F = 0.62, p = 0.57, Table 1).
Monthly changes in regional Seascape classes were correlated with mesoscale temporal patterns

Trends in recruitment
Over 60 taxa recruited on the Sargassum beds during the 24-month period (See Supplementary Table   S4), with a predominance of polychaete (32%), gastropods (30%), and holothurians (14%) ( Fig. 2; see Supplementary Table S4). Recruitment was continuous over the year with an average of 2 to15 taxa recruiting simultaneously every month (Fig. 2). Recruit abundance, composition, diversity, and richness varied significantly between seasonal and monthly time scales (p < 0.01, Table 2). We observed a lower abundance of recruits during winter (May or Jun), lower diversity during summer (Dec-Jan), and higher richness during fall and spring (Tukey HSD p < 0.05, see Supplementary Table   S5). Recruits peaked in abundance on two observed events that occurred during summer 2017 (Nov-Dec) and fall 2018 (May-Jun), when abundance was almost twice above average (Tukey HSD p < 0.005, see Supplementary Table S6; Fig. 2). The lowest diversity of recruits was registered during the peak abundance on summer 2017 (Dec) and spring (Mar) (Tukey HSD p < 0.001, see Supplementary   Table S5 Table 2 and Supplementary   Table S6; Fig. 2), but temporal variability was taxa specific ( Fig. 2 and Supplementary Fig. S3).
Monthly variations were significant for snails and limpets, and polychaetes had a higher number of recruits during summer and fall (Nov-Apr) (p < 0.05, see Supplementary Table S8 and S9; Fig. 2 and S3). Crabs and shrimps had a stronger recruitment during the fall (Mar), whereas limpets, holothurians, and bivalves did not show a consistent seasonal recruitment during the time sampled (p < 0.05, see Supplementary Table S8 and S9; Fig. 2 and Supplementary Fig. S3). Despite monthly and year differences in the number of recruits sampled, the number of samples represented well the regional species diversity after the first months of study (see Supplementary Fig. S4).

Association between recruitment and Seascapes
The abundance and richness of benthic recruits were temporally correlated with Seascapes (Tables 3   and Supplementary Table S10). Peaks in recruit abundances occurred one month after intrusions of Subtropical Gyre Transitions Seascape conditions (r = 0.51, p < 0.01, Table 3), and when Tropical waters were frequent in the region (r = 0.49, p < 0.01, Table 3). The peaks in abundance described above were particularly associated to a higher dominance of snails (r = 0.64, r = 0.41, p < 0.01, Table   3). Other taxa peaked following the dominance of Hypersaline Eutrophic waters, including crab and shrimp recruits (r = 0.44, p < 0.01, Table 3). Peaks in bivalve recruitment occurred after the onset of Subtropical Gyre Mesoscale Influence (r = -0.45, p < 0.01, Table 3). Therefore, recruit richness increased when coastal waters were more heterogeneous and with a weakening of the Warm, Blooms, High Nutrients Seascape (r = -0.45, p < 0.01, Table 3).
Seascape multivariate patterns explained 51 to 21% of the differences in taxa abundance and taxonomic composition of recruits (F = 2.64, p = 0.004, see Supplementary Table S11 and Fig. S5).

Modeling and validation
Overall, the generalized additive models described well the associations between Seascape classes or SST and benthic recruit abundance and richness ( Fig. 3; Table 3). According to these models, Seascape classes showed stronger correlations with recruitment than SST, with less deviance from predicting model (p < 0.05, Table 3). A higher model fitness was observed for snails (> 70%) and bivalves (> 90%) in association with the predominance of Tropical Seas and Subtropical Gyre Mesoscale Influenced seascapes, respectively. Validation with data from additional monitoring months showed similar trends to model predictions for richness and abundances of snails, crabs and shrimps, and bivalves (Fig 3). This data was within the spectrum of variability from previous months, with the exception of total abundances which were much higher.

Historical regional Seascape patterns
During the years of 2003 to 2019, nine Seascape classes were registered in the study region, characterizing both subtropical and tropical pelagic conditions (Fig. 4). Two Seascape classes were registered more frequently but differed in their spatial influence within the study region. Warm, Blooms, High Nutrients had a mean spatial coverage of 39% and predominated on the northern section (North of 20.03 o S), whereas Tropical Seas was predominant in the southern section (South of 20.03 o S) with a mean spatial cover of 46%. Both the average class and spatial coverage varied significantly between years, months, and seasons (p < 0.05, Table 1). Seascapes were highly seasonal with a higher cover of subtropical seascapes during summer and fall, whereas tropical nutrient-rich conditions dominated during late winter and spring (p = 0.01, see Supplementary Table S3; Fig. 4).
Seascapes patterns were historically distinct in 2015 and 2016 (p < 0.05, see Supplementary Table   S3), with increased spatial coverage of tropical waters along the northern and southern sections of the study area (Fig. 4).
The historical seascape database showed seven periods of potentially reduced benthic recruitment richness that were potentially associated with the dominance of Warm, Blooms, High Nutrients

Discussion
This is the first description and analysis of the NOAA Coast Watch Pelagic Seascape Habitat Classification for the Southern Atlantic on the Brazilian coast (20.3º to 19.8º S), revealing that warm and tropical Seascape classes are the dominant pelagic conditions in the region. As expected, we found a marked seasonal and regional variability in Seascape classes, with Hypersaline eutrophic We observed that the dominant Seascapes on the Eastern Brazil Marine Ecoregion correlated with patterns of larval recruitment of reef benthic species at multiple temporal scales, corroborating our initial hypothesis. This finding supports the interdependence of pelagic larval development and pelagic conditions such as water temperatures, food supply, and current transport in marine ecosystems [48,49]. Although these relationships have long been demonstrated for a range of coastal marine ecosystems, this study suggests that the Seascape classification may capture the dynamic nature of coastal pelagic oceanography and partially explain variations in larval recruitment across reef habitats in the South Atlantic at multiple spatial and temporal scales.
Seascape patterns provided several insights about potential larval mechanisms in the studied tropical rocky reef communities. The Seascapes revealed a lower recruitment rate in rocky reefs during warmer and higher nutrient conditions, which suggests that expected warmer temperatures would decrease recruitment success in this region [40,41]. Seascape analysis also revealed that intrusions of meteorological cold fronts change the coastal oceanography to subtropical conditions and likely promote rapid increases in recruitment through increased larval supply and dissolved nutrients. We found that the relationship between Seascape dynamics and recruitment was taxon-specific, revealing a range of reproductive strategies and larval dynamics on this area. Snails recruited preferentially with the onset of subtropical conditions, whereas Annelid worms recruited more during typical warm conditions of summer months. The spatial-temporal relationship of Seascape classification and recruitment also suggests that larval sources in the region can originate from large areas as larval recruitment was strongly associated was mesoscale seascape dynamics. This has important implications for the management of coastal reef ecosystems as it reveals that local assemblages are dependent on reproductive populations that are outside the managed boundaries of this MPA, and suggest that these external larval inputs are a strong driver of benthic biodiversity [50,51,52].
The capacity of identifying changes in past, current and future climatic conditions is crucial to prevent further declines on marine biodiversity [1, 2, 3, 4]. We examined the past Seascape dynamics within our study area and identified years with significant climate anomalies that significantly differed from the 19-yr time-series. During these anomalies, the predominant pelagic conditions Niño, and resembled current Seascapes with typical low recruitment events [53,54]. During the past years of climate anomalies, Warm, Blooms, High Nutrients Seascape covered most of the study region (60-100%) for 2 to 6 continuous months. The prolonged duration of warmer Seascapes suggests that observed changes in diversity during years of climate anomalies may be partially related to lower recruitment events. These observations have also strong implications for the longterm health of marine ecosystems if these warming events become more frequent and with prolonged duration [55,56].
Considering that this pelagic Seascape classification method is primarily derived from remote sensing measurements, it is important to be aware of the limitations in analyzing nearshore areas. For instance, river plumes in this region may add confounding effects to satellite images and overestimate the discrimination of warm blooms [57,58]. Likewise, monthly Seascapes may not detect short-term key larval transport features (e.g. upwellings, [30]), which could account for additional variability in recruitment patterns and associations with Seascapes. We shall also consider 2-yr duration of the study and the limited spatial coverage, which are fundamental aspects to be improved as Seascapes are applied to spatial planning and regional forecasting.
Our study revealed that Seascape ecology has a promising potential to improve the understanding of the long-term dynamics of coastal marine ecosystems. Seascapes provide a more holistic picture of the pelagic system than single-variable classification (i.e. based only in SST) and modeling of larval dynamics, which is a common approach in larval ecology (e.g., [28,59]). Consequently, Seascape offers many advantages in identifying past and current changes to marine ecosystems and has key applications in marine spatial planning in the search for climate-resilient protected areas. The historical assessment and cross-correlation of Seascape classification with current biological dynamics also offers a tool to model past changes in marine ecosystems and evaluate their ecological dynamics along broader spatial-temporal scales. In addition, Seascapes may prove valuable to evaluate variability in the essential ocean variables EOVs (e.g. invertebrate richness and abundance, [60]), and may aid the process of identifying temporal responses of marine ecosystems on a changing climate [61]. With the challenges ahead, Seascapes classification requires reasonable technical training with biological data acquisition, which could facilitate its application within marine observation networks.

Study area and sampling
This study was carried out at a coastal reef located inside two marine protected area in the Eastern Brazil Marine Ecoregion (Refúgio da Vida Silvestre de Santa Cruz and Área de Proteção Costa das Algas; environmental permit by Instituto Chico Mendes #57819-1; see Supplementary Fig. S1). This is a tropical region with an average air temperature of 25ºC that has experienced significant warming trends during the last four decades [39]. The coastal zone is characterized by sparce intertidal lateritic reefs with abundant macroalgal and rhodolith beds [41]. Coastal oceanographic conditions are typically influenced by E-NE winds from the South Atlantic high-pressure system, strong internal tidal currents, and E-SE wave swells [62,63]. Meteorological cold fronts occur periodically and influence the vertical mixing of the water column and wave action on the coast [63]. Episodic upwelling events occur mostly during spring and summer [62].
Our study site was located at a reef area that has been monitored by the Brazilian LTER site Coastal Habitats of Espírito Santo (HCES). The reefs in this region are mainly composed by ferruginous sandstone with small portions of biogenic calcareous (laterites), forming a very heterogeneous substrate with various shapes covered with dense and rich macroalgal beds [40,41,64]. Laterites occupy large littoral areas (tens of kilometers alongshore) extending from the upper intertidal zone to approximately 15-20 meters deep and about 5km offshore (Fig.1). Benthic invertebrate larval recruitment was monitored for 24 months (May 2017 to November 2019; Supplementary Table S1; [40]). Invertebrate recruitment was measured in macroalgal beds (Sargassum sp.), natural substrates in the region commonly present in the upper subtidal zone along the reef fringe and inside tidal pools.
Benthic invertebrate recruits were sampled monthly at the same location, where five samples of Sargassum fronds with multiple branches (approximately 300g of algae each) were manually collected during spring low tides. Sampling protocol was developed based on previous studies that investigated Sargassum associated fauna (e.g. [65]), which allowed collecting very abundant and rich assemblages of invertebrate recruits from resident and non-resident species (average 80 and maxima >1000 individuals per sample). Samples were placed in plastic bags, taken to the laboratory, and frozen at -20ºC for at least 24 hours, then washed with freshwater on a 100 μm sieve to collect the fauna (e.g. [65,66]). All recruits (post larvae, settlers, and first juveniles) were sorted, counted, and identified to the lowest possible taxonomic level under a stereomicroscope. For this study we considered post-larvae as the planktonic or benthic larval form that are capable of settling, settlers as metamorphosed larval stages that hold specific morphological structures to live in the benthos, and first juveniles as small sized individuals that are similar to the adult forms and may or not be able to reproduce; [67,68]). Discrimination between recruits and adults was carried out based on the size (adult average size for the species or group) and/or morphological characteristics particular of invertebrate post-larval and settlement phases (e.g., [69]). Recruit morphological differences included the presence of eye-spot and shell morphology for bivalves and gastropods, number and size of body segments for polychaetes, body morphology for cnidarians, echinoderms, barnacles, sponges, and tunicates, abdomen attachment for crabs, and shell thickness for polyplacophorans (See Supplementary Table S2).

Data analysis
Seascape monthly patterns were described using (1) differences in the Seascape % Cover (relative and absolute), calculated by the number of pixels of each Seascape class occupying the study region per month, (2) variations in the mean Seascape class, the quantitative value from averaging all class values available for the study region (average n = 76 pixels) and mesoregion (average n= 9800 pixels), and (3) seasonal Seascape maps to illustrate differences at the mesoscale. Month, seasonal, and year differences in Seascapes within the study region were evaluated using 3-way permutational multivariate analysis of variance (PERMANOVA; [71]), using distance matrices of relative % cover of Seascape classes during the sampling window (2017-2019) and in the long-term (2003-2019).
PERMANOVAs designs tested year (factor 1) fixed with 2 levels (2017 to 2018, 2018 to 2019), seasons (factor 2) fixed with 4 levels (fall, winter, spring, summer), and months (factor 3) random with 3 levels (3 months in each season), and were complemented by post-hoc all pair-wise comparisons. Variations in regional and mesoscale Seascape averages were compared by Pearson correlation analysis [72].
Recruitment patterns were determined from the total and relative abundance of recruits (recruits per kg of algae), recruits diversity (Shanon-Wiener index), and richness (S) in each period sampled.
Representativeness of samplings and low sample size were evaluated through taxa accumulation and rarefaction curves for the monitored years (2017-2018, and 2018-2019; [73,74]). Variability in recruitment between year, seasons, and months was tested by 3-way ANOVAs (ANOVA, with unbalanced replicates; [72]), complemented by a PERMANOVA to assess differences in recruit assemblages (taxonomic composition). Both ANOVAs and PERMANOVA followed the same design described above for the Seascape evaluation, followed by Tukey HSD and PERMANOVA post-hoc pairwise tests. The long-term trend in recruitment was assessed using regression analyses adjusted to the diversity and richness time series [72].
The association between recruitment patterns (abundance, diversity, and richness) and Seascapes classes was assessed by cross-correlation analyses. Correlations were considered significant when p < 0.01 and lags were between 0 to 2 months, based on the average development time for tropical coastal invertebrate larvae (weeks to months, [75]). Additionally, a canonical analysis of principal coordinates (CAP; [76,77]) was carried out to identify the class or set of classes with a higher potential to drive changes in recruit assemblages. CAP was performed comparing variations in recruit abundances for each taxon and the frequency of occurrence of each Seascape class. Generalized linear models (GAM; [78,79]) were applied to the dataset (May 2017 to April 2019) to examine Seascapes and SST as predictors of temporal changes in recruitment. The correlation coefficient, adjusted R-squared, and % deviance were used as parameters to evaluate model fitness and the levels of association between variables. The observed patterns for the following months (July to November 2019) were graphically compared to GAM trends. Only Seascapes that were significantly correlated to recruitment and corresponding lags were included in this analysis.
Long-term analyses of Seascape variability included: (1) graphical assessment of the monthly variations of Seascape classes along the study region; (2) estimation of average Seascape classes for the northern and southern part of the study region, using numerical Seascape classification to perform zonal averages for the available grid points North to the sampling site n = 5, and South to the sampling site n = 6); (3) calculation of the Seascape climatology (Jan to Dec) by averaging numerical classes from the 16 years available in the dataset for each month; (4) estimations of the residuals, the differences between observed Seascape numerical classes and climatological means; (5) identification regionally increase in the cover of Warm, Blooms, High Nutrients Seascape, Seascape conditions for potentially reduced recruitment, according to our results.
Data were Log x + 1 or square-root transformed when needed to fit the assumptions of ANOVAs, correlation analysis, and CAP (normality and homogeneity of variances), verified by Kolmogorov-Smirnov and Cochran tests; and square-root x+1 prior to PERMANOVAs. We assumed α = 0.05 and determined the significant p values using the Benjamin-Hochberg false discovery rate method [80,81]. Graphical and analytical processing were performed in Panoply 4.8.1 [82] for Seascape visualization and extraction, Numbers (Apple Inc.) for charts, and R project [83] for statistics (R packages 'stats' for general calculations, 'GAD' [84] and 'outliers' [85] for ANOVAs and pos-hoc tests, 'vegan' [86] for PERMANOVA, 'rich' [87] for ecologic indexes, and 'gam' for modeling [88].  Tables: Table 1. Results of the 3-way permutational multivariate analysis of variance (PERMANOVA) to compare the differences in Seascape % Cover (number of pixels) per class between months, seasons. Analysis included the sampling years (2017-2019) and the previous 16 years (2003-2018). Note: F for statistic, significant results (p < 0.05) are in bold.  Table 2. Results of 3-way analyses of variance and permutational analysis of variance to compare the differences in recruit assemblage abundance (total number of recruits), composition (number of recruits per taxa), diversity (Shannon-Wiener Index S-W), and richness (number of species/taxa) between months, seasons, and years. Note: F for statistic, significant results (p < 0.05) are in bold.