Non-deterministic reef fish community assembly in an upwelling-influenced transitional subprovince of the southwestern Atlantic

1 Introduction

Community structure theory has long acknowledged the importance of historical and regional processes to explain local species diversity (Ricklefs and Schluter, 1993). Embedded in these spatio-temporal scales are deterministic and stochastic mechanisms involving coexistence, migration and speciation that potentially influence local communities. It follows that communities are structured around different spatio-temporal scales according to the availability of species in the regional pool and the processes underlying their establishment and persistence. These, in turn, determine their local distribution, abundance and functional relationships (Ebeling and Hixon, 1991). According to this spatio-temporal hierarchy, higher scale properties constrain lower-scale processes and large-scale patterns may emerge from local relationships (MacNeil and Connolly, 2015). Eventually, scales and processes in community assembly will be structured by a combination of processes including selection, ecological drift, speciation, and dispersal, highlighting the importance of clearly defining principles before looking for patterns (Vellend, 2010). It is well accepted that these processes act concertedly to structure communities, but their relative roles are still a matter of contemporary research.

1.1 Community assembly: local habitat filtering or random process?

Under niche theory, the regional pool provides species that will immigrate and undergo environmental or habitat filtering to form local assemblages (Cornell and Lawton, 1992). On the other hand, neutral theory predicts that species assembly (and richness) forms as a random process, assuming that the per capita relative fitness within a trophic level are equivalent, and the occurrence of a strong local contribution through ecological drift (Hubbell, 2001). In a nutshell, the main drivers of community assembly equate into species composition responding to speciation and dispersal, whereas relative abundances are driven by selection, ecological drift, and ongoing dispersal (Figure 1). Disturbance can reduce abundances in a community, enhancing ecological drift and impacting speciation. So, in diverse environments, biogeographical processes, environment filtering and local adaptation are expected to exert a detectable influence on community assembly. It is necessary, however, to factor out the different drivers of species abundance by looking at the relative roles played by regional pool of species, abiotic disturbance and biotic interactions to the assembly of local communities.

Figure 1 Schematic representation of the main processes underlying community assembly (Hubbell, 2001; Vellend, 2010). Insert map depicts spatial patterns of surface chlorophyll-a concentration induced by coastal upwelling in Arraial do Cabo, southwestern Atlantic. Red indicates high surface concentration of chlorophyll-a and blue indicates low concentration.

There is evidence that neutrality is important in both tropical and temperate marine environments, with dispersal limiting species assemblage in tropical sites and environmental filtering becoming more important for temperate communities (Ford and Roberts, 2018). In fact, at local scales (<100 m), the physical structure of the substrate may act as a protection to environmental variability or predation and plays a role in organizing reef fish communities (Tuya et al., 2011). The limited spatial scale also makes species more susceptible to unpredictable environmental conditions and natural disturbances. In vulnerable groups, the loss of functionally similar species induced by disturbances may cause functional redundancy to decline (Brandl et al., 2016). Therefore, clarification about the role of historical contingency (phylogeny), dispersal, traits and niche in shaping local communities along environmental gradients is necessary to allow the scaling-up of assembly processes.

1.2 Reef fish community: the regional pool

Within the Brazilian Biogeographic Province (BBP), reef fish are distributed from the Amazon River mouth in the north to the Santos Basin in the south (Santa Catarina state) encompassing a large latitudinal gradient (2°N-29°S). The Brazilian Province can be subdivided into six biogeographic sub-provinces (Pinheiro et al., 2018). Differences in species richness among sub-provinces have been attributed to factors ranging from the southern branch of the South Equatorial Current, coastline length, habitat diversity and differences in temperature driven by upwelling influence in the southern regions (Floeter et al., 2001; Floeter et al., 2007; Floeter et al., 2008; Mazzei et al., 2017; Pinheiro et al., 2018). The coastal central east-southeast sub-province (between 12–25°S) has the highest richness with 330 species of reef fishes, where 49 species are exclusive. This sub-province has been further segmented in two subgroups according to species occurrence, one to the north centered at Abrolhos (state of Bahia) and another to the south, centered at Arraial do Cabo (state of Rio de Janeiro). In fact, this is a tropical-subtropical transition zone (Pinheiro et al., 2018), where 27% of species are not found in the north-northeast sub-province close to the Caribbean. A high diversity of habitats, ranging from tropical coral and coralline algal reefs to the north and subtropical rocky reefs towards the south, makes this region to harbor the largest number of endemics (16 species) in the BBP. Besides the largest coral reef buildups in the South Atlantic (the Abrolhos reefs, 19°S), there are two important upwelling regions centered at Cabo Frio (23°S) and Cabo de Santa Marta (28°S). The former is restricted to the central portion of the sub-province close to Arraial do Cabo and the latter occurs at the southern limit. The ecological speciation and genetic structure of these sub provinces are likely influenced by habitat diversity, temperature, primary productivity and larval transport (da Cunha et al., 2014), also driven by western boundary currents and persistent mesoscale eddies (Endo et al., 2019; Gouveia et al., 2021). It has been advocated that this transition subprovince should be considered a secondary center of biodiversity of reef fishes on the basis of a relatively high species richness and endemism (Pinheiro et al., 2018).

At the scale of sub-provinces, the composition of reef fish communities tends to be driven by habitat use, body size, depth range and diet (Pinheiro et al., 2018). However, it is not yet clear what is the shortest scale threshold that it holds true for communities submitted to natural environmental changes (Cordeiro et al., 2016). In some regions, like in Western Australia, the assembly of temperate marine fish communities can be mostly driven by local-scale processes (species interaction with biotic and abiotic environments), whereas dispersal limitation is more important for tropical communities (Ford and Roberts, 2018). The concerted influence of neutral and niche theory processes has been demonstrated to structure cryptobenthic fish assemblages, suggesting that they should not be viewed as mutually exclusive (Ahmadia et al., 2018). Therefore, areas submitted to natural environmental disruptions, such as coastal upwelling with its typical alternation of cooler (upwelling) and warmer (non-upwelling) conditions are suitable test sites for studies designed to investigate community assembly processes.

Here, we present a process-oriented study aimed at testing the hypothesis that niche limitation induced by coastal upwelling is the main process influencing reef fish community assembly of Arraial do Cabo. These communities are located in the coastal central east-southeast Brazilian sub-province under the influence of a local upwelling and are naturally distributed between rocky reefs. This geographic setting can be used as a natural laboratory to test the hypothesis that changes in environmental conditions triggered by coastal upwelling control the assembly of reef fish communities. Reef fish communities from Arraial do Cabo were partially described by Ferreira et al. (2001) and are now more thoroughly assessed, focusing on community assembly drivers in this natural sheltered vs. exposed environment. We used a hierarchical modeling scheme that allowed disentangling the influence of phylogenetic history, traits and biotic interactions embedded in a spatio-temporal context (Tikhonov et al., 2019; Ovaskainen and Abrego, 2020). Our results highlight the role played by environmental disturbance and biogeographical history on the dynamic of these communities in a tropical-subtropical transition subprovince.

2 Data and methods

Arraial do Cabo (-22.966, -42.0288) is a cape located in the southwestern Atlantic (Brazil) formed by a granitic rocky embayment in the lee of Cabo Frio Island (Figure 2) that is under the influence of a local upwelling. This is a seasonal phenomenon that occurs mainly during the austral spring and summer (Valentin, 1984) resulting from the offshore wind-driven Ekman transport of coastal waters forced by the persistence during a few days of the northeast surface winds. The upwelling circulation pattern is characterized by a coastal jet, which interacts with the topography resulting in gradients of relative vorticity. These gradients create vertical velocities that bring the deeper, colder, nutrient-rich South Atlantic Central Water (SACW) closer to the surface (Rodrigues and Lorenzzetti, 2001; Castelão and Barth, 2006; Mazzini and Barth, 2013). Direct environmental changes caused by the coastal upwelling are negative sea surface temperature (SST) anomalies and an increase in primary productivity (Matsuura, 1996). The cape divides the coast into north-south and east-west sections, providing a natural set of sheltered and exposed substrates. Sheltered sites are shallower, with depths of up to 15 m, and the benthic cover includes a diverse epilithic algal community and scleractinian corals typical of tropical systems (Ferreira et al., 2001; Cordeiro et al., 2016). On the other hand, exposed sites are usually deeper (up to 30 m depth) and the benthos is dominated by macroalgae, a thin epilithic algal community and lacks scleractinian corals (Cordeiro et al., 2016). The whole region is a Marine Extractive Reserve of Arraial do Cabo, created by Federal Decree in 1997, with an area of 51,601.46 hectares (ICMBio, 2020), which limits only local fishermen to explore the marine resources, while no take areas are inexistent. We selected six sampling sites with three representing sheltered sites: Anequim (-22.980571, -41.984114), Pedra Vermelha (-22.986476, -41.992504), Porcos (-22.963476, -41.993369), and three exposed sites: Franceses (-22.981626, -42.037293), Sometudo (-22.99 2232, -42.016005) and Ingleses (-23.006032, -42.008244). Fish abundance was obtained from underwater visual censuses (UVCs) consisting of strip transects measuring 20 × 2 m (40 m²) (Ferreira et al., 2004; Floeter et al., 2007), surveyed between 2017 and 2021, from which a total of 153 species were recorded.

Figure 2 Map showing the location of reef fish sampling sites in Arraial do Cabo. Red dots denote sheltered and blue dots denote exposed sites to upwelling influence.

2.1 Environmental data

Environmental predictors of niche assembly driven by coastal upwelling were obtained from daily remote sensing data on sea surface temperature (SST, °C), surface chlorophyll-a concentration (as a proxy for primary productivity; SCC, log10 mg/m−3) reanalysis data for meridional surface (at 9 m) and the meridional wind speed component (MWC, m/s-1), covering the period between January 2017 and December 2021. The zonal surface wind velocity component was also considered, but was later dropped from the analysis because of the weak response of species niche. The foundation SST level 4 product was acquired from the Multi-scale Ultra-high Resolution (MUR) analysis that combines data from infrared and microwave satellite sensors and in situ observations with a spatial resolution of 0.01° x 0.01° (approx. 1.1 km) with daily resolution (https://doi.org/10.5067/GHGMR-4FJ04) (Chin et al., 2017). The SCC data were obtained from different ocean color satellite sensors merged in the interpolated Global Ocean Color level 4 product from the Copernicus program, with a spatial resolution of 4 km and daily temporal resolution (https://doi.org/10.48670/moi-00281). Finally, MWC data from ERA5 reanalysis (a combination of model data and observations) were downloaded from the European Center for Medium-Range Weather Forecasts (ECMWF), with spatial resolution of 0.25° x 0.25° and hourly temporal resolution (https://doi.org/10.24381/cds.adbb2d47) (Hersbach et al., 2023).

Environmental covariates were pre-processed for a proper characterization of the upwelling environment, taking into consideration its spatio-temporal variability. First, we looked at a ten-years long series of satellite SST and SCC data to assess the upwelling variability at the surface. We observed the existence of warmer/cooler SST anomalies upstream/downstream of the cape and a separation of the upwelling coastal jet at Arraial do Cabo, as previously described by Mazzini and Barth (2013). This feature roughly separates sheltered and exposed sampling sites to the north and south of the cape, respectively. Based on that, we selected two polygons with approximately 5.655 km² each, representing the environmental condition of sheltered and exposed sites (Figure S1), from which we extracted environmental data used in the analysis. Second, we considered a time window of five days as a minimum necessary persistence time of the northeast wind to develop a mature upwelling condition. Third, we adopted the following strategy to calculate the environmental variables: (1) identification of sampling occasions at each sampling site, (2) download of daily satellite products for a five-day period back in time counting from (and including) the sampling date, (3) calculation of the five-day average SST prior to each sampling date, (4) calculation of the average SCC above the 75% percentile of the 5 days prior to each sampling date, and (5) calculation of the circular mean of MWC for the 5 days before the sampling date. This integration time is long enough to allow the development of the baroclinic coastal jet, as shown by Rodrigues and Lorenzzetti (2001) in their numerical upwelling experiments. These factors are not completely independent, but they are complementary predictors representing important thermal and energetic constraints to reef fishes (Quimbayo et al., 2019). The impact of multicollinearity in our models is expected to be low, with |r| ≤ 0.6 for all covariates. Upwelling events also differ at the surface due to natural spatio-temporal variability, with negative SST anomalies generally observed during high primary productivity events (Rodrigues and Lorenzzetti, 2001). Other environmental covariates used to describe coastal upwelling such as nonlinear physical interactions, sea breeze (Franchito et al., 1998) and wind stress curl (Castelão and Barth, 2006) were not considered in the present study as they are of secondary importance to the upwelling process.

2.2 Species niches and community assembly modeling

Species niches were modeled using a multivariate hierarchical generalized linear model fitted with Bayesian inference from the Hierarchical Modeling of Species Communities (HMSC) package in R (Ovaskainen et al., 2017; Tikhonov et al., 2019; Ovaskainen and Abrego, 2020). This is a correlative model designed to allow different kinds of ecological data to be used for making inferences about community assembly processes. As a process-oriented model, it demands from the user a previous knowledge of the ecological system that will be modeled to orient the necessary statistical decisions and support interpretation (Ovaskainen and Abrego, 2020). At the core of the HMSC model there are six data matrices: species abundances, environmental data, species traits, phylogenetic data, the spatio-temporal structure (study design) and spatial coordinates (detailed below). The pre-processing of abundance data from the UVC involved the trimming of the 153 species recorded to retain only those that: (1) were detected in all years, (2) were present in at least 5 censuses for at least 50% of the total years, or (3) reached an average abundance greater than 0.2 individuals per square meter per year (sensu Garcia, 2021; Gomes et al., 2023). This resulted in 72 species that were retained, of which six species are exclusive to the exposed area, 20 species are exclusive to the sheltered area, and the remaining 46 species occur in both areas. As expected, the resulting abundance (response) matrix Y has many zero entries, so we applied the hurdle approach, which tackles the problem of too many zeroes in the original species abundance data, separating it into two components (Barry and Welsh, 2002). The first one is a log-transformed abundance matrix conditioned to presence for each species modeled with a linear function, and the second one is a presence-absence matrix modeled with a probit link function.

A preliminary assessment of the species matrix Y showed that some species have been consistently recorded exclusively in the sheltered or the exposed habitats. Based on that, we decided to investigate if this could be related to an environmental filter by repeating the hurdle procedure on a separate subset of the same abundance data composed only of a list of all species that were found exclusively on the sheltered and the exposed sites. This resulted in a total of four models with an explicit spatio-temporal structure that accounts for site (Franceses, Ingleses, Sometudo, Pedra Vermelha, Porcos and Anequim), location (exposed and sheltered), month and year. This hierarchy has to be defined early in the data structure, called study design, in which the explanatory variables are the environmental covariates (SST, SCC and MWC) organized in matrix X and treated as fixed effects, while the spatio-temporal structure (sampling site lat/lon, exposed vs. sheltered location, month and year) was treated as random effects. Random effects are used to account for spatial and temporal dependencies and to indicate pairs of species that co-occur at the specified spatial and temporal scales (Ovaskainen and Abrego, 2020). Other input data are life history traits for all species defined in matrix T compiled from Quimbayo et al. (2021) and the phylogenetic correlation matrix C. Different model runs were used to test a total of eight traits (trophic level, minimum temperature of occurrence, home range, level of water, diet, maximum body size, depth range and geographic distribution index) from which maximum body size, diet, and geographic range index were kept for the present analysis (Table S1). The phylogenetic history of the community was computed as the fraction of evolutionary time shared by species pairs obtained from the molecular phylogeny of teleosts of Rabosky et al. (2018), available at https://fishtreeoflife.org/. The complete phylogeny was pruned down to 72 species of the Y matrix and the mean branch length is computed for each taxon from 100 phylogenetic trees using the phytools package running in RStudio (R Studio Team, 2020). The species Serranus aliceae (Carvalho-Filho and Ferreira, 2013) passed the cut criteria applied, however, because it is not included in the phylogenetic tree of Rabosky et al. (2018), was replaced by the phylogenetically closest species, Serranus tortugarum (Carvalho-Filho and Ferreira, 2013).

All four models (two for the complete set of species and two for those exclusive to each location) were fitted from the default prior distribution adjusted to the scale of variation in the data. The posterior distribution was sampled with four MCMC (Markov chain Monte Carlo) chains, each with 150,000 iterations, where the first 50,000 were discarded as burn-in and the remaining interactions were thinned by sampling the posterior distribution at an interval of 100 samples. This generated a total of 4,000 chain samples, from which convergence was analyzed (Ovaskainen and Abrego, 2020). The explanatory power of the linear model for abundances was evaluated using the squared Pearson correlation coefficient between fish abundance and model prediction (R²) and the square root of the mean error (RMSE).

The discrimination power of the presence-absence model was assessed using the coefficient of discrimination Tjur R² (predicted probability of occurrence where species occur minus where it does not; Tjur, 2009), the area under the curve (AUC) and the square root of the mean error (RMSE). The HMSC allows the estimation of realized niche of species denoted by the β parameter, the phylogenetic signal between species denoted by the ρ parameter, the influence of traits on niches denoted by γ parameter and species associations denoted by the Ω parameter (see Table 4.4 in Ovaskainen and Abrego, 2020). Species associations (denoted by Ω) are computed at all random levels from residual covariances after controlling for shared responses to the environment (Tikhonov et al., 2019). All results have a level of statistical support based on a posterior probability > 0.9. The predictive power was computed from cross-validation using the K-fold cross-validation method (Hastie et al., 2009), based on 4 folds. Sampling sites were randomly grouped into four folds and abundance predictions were computed for each separately, avoiding the use of data from the focal fold to be used for fitting and prediction (Ovaskainen and Abrego, 2020). A short description of model functionalities has been reported by Gomes et al. (2023) that applied a similar setup to communities of reef fishes from remote oceanic islands of the southwestern Atlantic.

The hurdle model applied to the entire community and to the exclusive species needed to be fitted separately for abundance (conditioned to presence) and presence-absence data. The convergence parameters for the four resulting models returned effective sample sizes (ess, the minimum size of posterior samples for an adequate sampling of a given parameter) close to 4,000 and the potential scale reduction factors (psrf) close to 1 for the linear and probit models. The ess close to one indicates all chains give consistent results, while the psrf showed little autocorrelation among consecutive samples and approached the actual MCMC sample size (Ovaskainen and Abrego, 2020) (Figure S2). In addition to the convergence analysis, the explanatory and predictive power was calculated for each species (Tables S2, S3). For the presence and absence model, the AUC, Tjur R² and RMSE coefficients were used, for the linear model, the R² and RMSE were used, where each coefficient varies from 0 to 1, with 1 being defined as the best fit and 0 as bad.

3 Results

The average explanatory power for the abundance model (R2) of all species was 0.42 (RMSE = 0.36) and 0.41 (RMSE = 0.31) for the exclusive species (Figure S3). The explanatory power of the presence-absence model (TjurR2) was 0.39 (AUC = 0.92, RMSE = 0.28) and 0.46 for the exclusive species (AUC = 0.96, RMSE = 0.23) (Figure S4). The mean predictive power of the abundance model was low for both all species (R2 = 0.1, RMSE = 0.46) and the exclusive species (R2 = -0.01, RMSE = 0.46) (Figure S5) models. The predicting power of the presence-absence model of all species was higher, scoring TjurR2 = 0.18 (AUC = 0.68, RMSE = 0.36) and for the exclusive species a TjurR2 = 0.21 (AUC = 0.72, RMSE = 0.32) (Figure S6). Differences in species abundance response to environmental covariates and random effects across all species were measured by the average variance partitioning (Figure 3). We found a stronger contribution of random effects, specifically sampling sites with 24% and location with 20.7% of the explained variance, and weaker fixed effects from the environmental covariates, with 7.1% for SCC and 5.4% for SST. This contrasts with the variance partitioning of the presence-absence model where environmental effects had a stronger contribution, with the highest mean of 20.1% for SCC. This symmetric behavior is evident in the columns, where effects due to sampling years and site are bigger for the abundance model. Species like Dactylopterus volitans, Holocentrus adscensionis (Figure 3, columns 1 and 55, top panel) responded strongly for sampling year, whereas Pronotogrammus martinicensis and Scarus zelindae (Figure 3, columns 20 and 46, top panel) responded to the site effect. The predominant environmental effect identified in the presence-absence model is generally more constrained, as seen in columns 5 and 51 (Figure 3, bottom panel) representing Gnatholepis thompsoni and Synodus intermedius, respectively. It should be noted that the random effect assigned to sheltered and exposed sites (location) have comparable variance contribution among species in the abundance and the presence-absence models (20.7% and 18.1%, respectively) exemplified by Stegastes fuscus (Figure 3, column 69, bottom panel).

Figure 3 Species-specific and global variance partitioning (VP) among environmental (fixed) and spatio-temporal (random) effects for the abundance model of all species (top panel) and presence-absence (probit) model (bottom panel). MWC, meridional wind speed component; SST, sea surface temperature; SCC, surface chlorophyll-a concentration. Correspondence between numbers and species name is shown in the Supporting information (Table S4).

Looking at the influence of upwelling conditions on all species (β parameter) it is clear that SCC and MWC exert a negative influence on both their abundance and occurrence (Figure 4, left and right panels, respectively). Some closely related species like Gymnothorax moringa, Synodus intermedius, S. synodus and Chilomycterus spinosus spinosus, Canthigaster figueiredoi showed negative responses to SCC in the abundance and presence-absence models. Although restricted to fewer species, the same negative relation was found for SST for a group of related species exemplified by Parablennius pilicornis, P. marmoreus, Emblemariopsis signifer and Malacoctenus delalandii. Positive relation between SCC and abundance was observed for Caranx latus, Decapterus macarellus and Heteropriacanthus cruentatus. These are highly mobile, pelagic species that will likely approach coastal habitats during upwelling events when SCC and food availability increase. Not all species living exclusively in sheltered and exposed sites responded to upwelling effects, but those who did, effects are mostly negative, as will be detailed below.

Figure 4 Heat map indicating species-specific response of abundance to environmental covariates (species niche, estimated β parameter) for the abundance model of all species (left panel) and presence-absence (probit) model (right panel), with posterior probability > 0.9. Blue indicates a significant negative association and red a significant positive association. SCC, surface chlorophyll-a concentration; SST, sea surface temperature; MWC, meridional wind speed component.

After extensive tests with different sets of life history traits to determine their role in explaining species abundance response to upwelling (γ parameter), we found a limited set of traits (Table S1) that helped explain species niche, limited to the abundance model (Figure 5). Interestingly, the geographic range index showed a positive relation with abundance response to SCC and SST. The positive relation found between the intercept and macroalgal feeder can be interpreted as the relation of this specific diet to the mean abundance of species in the community when SCC and SST increase. The presence-absence model failed to detect any influence of traits on niches above the 0.9 threshold for the statistical support (posterior probability). Similarly, the phylogenetic parameter ρ explained around 95% of species niches in the abundance model, but only 15% in the presence-absence model.

Figure 5 Heat map depicting the links between species traits to species niches (estimated γ parameters) for the abundance model of all species with posterior probability > 0.9. There were no detectable links between species traits and niches for the presence-absence (probit) model. Blue indicates a significant negative association and red a significant positive association. SCC, surface chlorophyll-a concentration; SST, sea surface temperature; MWC, meridional wind speed component.

The residual species associations (or co-occurrence, denoted by the Ω parameter), when upwelling-related covariates fail to explain their abundances, have been computed for different spatio-temporal scales (random effects). This can be interpreted as a proxy for biotic filtering, with the caveat that such interpretation applies only to the environmental covariates used in the models. Both abundance and presence-absence models showed a stronger response to location, and less so for sampling sites (Figure 6). This indicates that the seasonal upwelling recurrence in AC exerts a control in communities between sheltered and exposed sites, defining how species become more or less abundant among sites (Figure 6A) and among locations (Figure 6B) or move between locations (Figure 6F). Temporal effects had less influence on species associations, indicating local interactions are consistent within months (Figures 6C, G) and years (Figures 6D, H).

Figure 6 Heat map indicating associations between species (estimated Ω parameters) for the abundance model of all species (A–D) and presence-absence (probit) model (E–H), with posterior probability > 0.9. Blue indicates a significant negative association and red a significant positive association. Correspondence between numbers and species name is shown in the Supporting information, Table S5.

We further detailed the above analysis by running the model using only those species that were exclusively found in the sheltered and exposed habitats (locations). The variance partitioning among fixed and random effects (Figure 7) are in general agreement with the pattern observed when modeling the complete set of species (see Figure 3). Specifically, the symmetric pattern, i.e. high random/low fixed factors for abundance model – low random/high fixed factors for presence-absence model. The influence of random effects on abundance predominates with the mean for site and location effects around 19%, followed by year with 17% and month with 14.8% (Figure 7, top panel). However, a relatively strong influence of SCC was also found, explaining up to 13.1% of the across-species variance. There is a relatively homogeneous per species response to all factors in the abundance model with only a few species having a strong response to site: Scarus zelindae (Figure 7, column 15, upper panel) and Dules auriga (Figure 7, column 7, upper panel) that occur only in the sheltered and exposed habitats respectively, and location: Stegastes fuscus (Figure 7, column 25, upper panel), a species that occurs only in the sheltered region. Note that the percentage of global variance partitioning for site and location, in the abundance model, have very close values (19.8% and 19%, respectively). The presence-absence model showed a stronger response of species to the environment, with 21.1% of total variance explained by SCC and 18.2% for MWC. Adding to the environmental niche variances there is also a strong random effect related to sheltered/exposed habitats denoted by 23.2% of total variance explained by the location factor. The MWC covariate (18.2% of the total variance) becomes quite expressive for Decapterus macarellus (column 20, lower panel) and Caranx crysos (Figure 7, column 21, lower panel), both species frequently observed in the sheltered area, despite not being limited to it.

Figure 7 Species-specific and global variance partitioning among environmental (fixed) and spatio-temporal (random) effects for the abundance model of exclusive species (top panel) and presence-absence (probit) model (bottom panel). MWC, meridional wind speed component; SST, sea surface temperature; SCC, surface chlorophyll-a concentration. Correspondence between numbers and species name is shown in the Supporting information, Table S4.

The environmental niche (β parameter) estimated for the exclusive species showed that SCC, SST and MWC have mostly negative effects on the abundance of species. In fact, most species for which environmental niche has a statistical support > 0.9 are those found exclusively inside the upwelling sheltered area. Positive SST relations with abundances are restricted to Stegastes pictus and S. fuscus, that exclusive of the sheltered area, and Serranus tortugarum and Heteropriacanthus cruentatus exclusive from the upwelling exposed areas (Figure 8, left panel). Likewise, the presence-absence model showed mostly the negative influence of SCC and MWC on the occurrence of species, meaning that the increase in SCC and northeast wind velocity reduce the occurrence of these species (Figure 8, right panel). The contribution of the phylogenetic signal to species niche, denoted by the ρ parameter, explained 88% of the assembly of exclusive species in the abundance model, and 34% in the presence-absence model.

Figure 8 Heat map indicating species-specific response of abundance to environmental covariates (species niche, estimated β parameter) for the abundance model of exclusive species (left panel) and presence-absence (probit) model (right panel), with posterior probability > 0.9. Blue indicating a significant negative association and red a significant positive association. MWC, meridional wind speed component; SST, sea surface temperature; SCC, surface chlorophyll-a concentration.

A very limited influence of life history traits on species niche (γ parameter) was also found for those living exclusively in sheltered or exposed locations (Figure 9). Geographic range index and diet (planktivore) were found to influence niches only in the presence-absence model, differently from what was observed when modeling the complete set of species (see Figure 5). Only SCC appears to influence niches with a positive relation, indicating that more planktivores and wide geographic range species are more present during upwelling events, when SCC tend to increase.

Figure 9 Heat map depicting the links between species functional traits to species niches (estimated γ parameters) for the presence-absence model of exclusive species, with posterior probability > 0.9. There were no detectable links between species traits and niches for the abundance model. Blue indicates a significant negative association and red a significant positive association. SCC, surface chlorophyll-a concentration; SST, sea surface temperature; MWC, meridional wind speed component.

The residual correlation not explained by the environmental niche was also computed for the location-specific species, modeled for abundance and presence-absence, to check if random effects could help explain why some species were observed only in a specific location (Figure 10). The result of both models reinforces the importance of this habitat effect over other random effects (site, month and year), and agrees with the influence of upwelling exposure estimated for the complete set of species (see Figure 6). Again, random effects capturing spatial dependency (location and sites) achieved a noticeable influence on both species abundance and occurrence (Figures 10A, F). Sampling location showed greater influence, making it possible to perceive the relationship of species exclusiveness to upwelling exposure (Figure 10F – 2 to 6) responding negatively to the presence of species exclusive to the sheltered region (Figure 10F – 8 to 23 and 25).

Figure 10 Heat map indicating associations between species (estimated Ω parameters) for the abundance model of exclusive species (A–D) and presence-absence (probit) model (E–H), with posterior probability > 0.9. Blue indicates a significant negative association and red a significant positive association. Correspondence between numbers and species name is shown in the Supporting information, Table S6.

A preliminary 10-year analysis highlighted some minor ocean warming over the studied period (2017–2021), shown in Figure S7. This provided the background on which we analyzed the covariates between 2017 and 2021, where we found the expected inverse relationship between SST and SCC (Figure 11, top panel). Note that the exposed location is generally cooler than the sheltered one but peaks of SCC may occur in both locations with higher values found in the exposed habitat. The negative values of MWC, denoting upwelling-favorable northeast wind, also agree with the expected increase in SCC (Figure 11, bottom panel), which confirms the general upwelling pattern where strong and persistent northeast wind induces the increase of SCC and a cooling of SST (Matsuura, 1996; Rodrigues and Lorenzzetti, 2001).

Figure 11 Comparative environmental analysis performed using daily SCC from the copernicus-GlobColour data, daily MUR for SST data, and ERA 5 reanalysis data for MWC. Top panel shows the relation between SCC and SST values, where periods with higher SCC are generally associated with lower SST. Bottom panel shows that high SCC values tend to occur simultaneously with negative MWC values, indicative of the predominant winds coming from the northeast.

4 Discussion

Ecologists often debate the predominance of niche limitation, ecological drift and migration-driven control on community assembly (Vellend, 2010; Gilbert and Levine, 2017). We show evidence of the joint action of migration-niche mechanisms on reef fish species living in sheltered and exposed habitats to coastal upwelling. This is inferred from the dominance of spatio-temporal structure and limited influence of life history traits in abundance models. The spatio-temporal dependency structure observed from the residual variation of abundance was dominant in both the entire community and the exclusive (sheltered vs. exposed) species models. The separate modeling of location-exclusive species abundance (designed to detect differences in assembly rules between exposed or sheltered habitats) failed to show a clear niche-determined structure. However, environmental filtering became important to define reef fish occurrence, with a strong influence of upwelling-induced chlorophyll-a concentration. Even so, the influence of SST in niche assembly was lower than expected in an upwelling scenario, indicating the bulk of species are generalist to local changes in oceanographic conditions (e.g. temperature plus wave surge). It has been suggested that roving herbivorous fishes in Arraial do Cabo may be able to avoid exposed areas with a high overlap of their distribution between sheltered and exposed habitats (Cordeiro et al., 2016). In fact, Arraial do Cabo communities are representative of a tropical to temperate gradient encapsulated in a short spatial scale (e.g. 10 km), mimicking the Brazilian coast latitudinal gradient pattern (Ferreira et al., 2004; Pinheiro et al., 2018).

The predominance of spatio-temporal structure in the abundance models and environmental niche in the presence-absence models suggests that coastal upwelling is not the major driver of community assembly in Arraial do Cabo. Despite upwelling initiated and having stronger influence on exposed habitats, short time-scale dynamics (days to weeks) of associated oceanographic (currents and tides) and atmosphere (local wind regimes) conditions may enhance water exchange from the exposed to sheltered rocky shores (Ferreira et al., 2001). It is also possible that upwelling-induced disturbances might favor stochastic processes, including the dispersal of individuals between exposed and sheltered habitats. We, therefore, hypothesize that specific assembly opportunities are initiated during each upwelling event by allowing some species to move in and out of exposed and sheltered habitats according to environmental cues such as changes in SCC and SST. The fact that some of the species that appear as exclusive from the sheltered (e.g., Decapterus macarellus and Aluterus monoceros) and exposed (e.g., Heteropriacanthus cruentatus) locations have high geographic range index (between 133 and 374) also points to the likelihood that dispersal assembly from the regional pool of species repeatedly contributes to the community via short distance migration. All those species possess high dispersal capacity which enable them to occupy either sheltered or exposed habitats. Evidence detailed below suggests that, in this upwelling-dominated coastal habitat, dispersal assembly and phylogenetic history are the main drivers of reef fish abundance. Niche assembly would mostly result from upwelling acting as a disturbance event that affects the occurrence of some species, specially in the sheltered area. The Cabo Frio upwelling can, therefore, be viewed as a historical contingency that affects community assembly with past events varying in time of occurrence, intensity and spatial extent (Fukami, 2015). The compounded effects of dispersal and niche controls possibly contribute to the overall spatio-temporal stability of reef fish abundance on a local spatial scale.

The spatio-temporal structure (measured as random effects) is two to three times stronger in regulating species abundances compared to occurrences depicted from the presence-absence models (Figures 3, 7 bottom panels). Therefore, from the standpoint of fish abundance alone, the community assembly process is only weakly controlled by the environment, despite the putative upwelling exposure gradient (Ferreira et al., 2001). Therefore, we interpret the high importance of the spatial structure as resulting from stochastic processes related to the historical colonization of the coast and migration from the regional pool (Hubbell, 2001; Gilbert and Lechowicz, 2004). This assumption is further supported by a strong correlation between phylogeny and the abundance response to the environment (ρ parameter). The mean ρ parameter for the entire community abundance model is 0.94 (SD = 0.05) and for the exclusive species model is 0.89 (SD = 0.15). Interestingly, the percentage contribution of temporal structure to total abundance is of the same order of the spatial structure, as observed in the entire community model (18.1% for month and 20.3% for year) and the exclusive species model (14.8% for month and 17% for year). This temporal structure could result from the seasonal nature of the Cabo Frio upwelling that occurs mainly in the austral spring and summer months. During this period of the year, the bottom thermal front associated with the cooler South Atlantic Central Water (SACW) under the surface Ekman layer separates inner and mid-shelf waters 10-20 km from the coast (Castro Filho and Miranda, 1998). This interpretation should, however, be taken with some caution because upwelling occurrence and intensity may vary among different years. Also, underwater visual censuses were performed at different months in each site, and were not coincident among years (Table S7), eventually weakening the upwelling seasonality in the temporal structure. It is, nevertheless, an indication of temporal similarity in species abundance at seasonal and interannual scales. Another candidate process contributing to the temporal structure of the abundance of reef fish communities is the observed occurrence of annual recruitment pulses of species such as Diplodus argenteus, Haemulon aurolineatum and Abudefduf saxatilis. These pulses are mediated by settlement and post settlement processes (Ault and Johnson, 1998) and may lead to samples of the same year being slightly more similar than those from the same month, as the influence of recruitment on abundance may vanish in two to three months.

The environmental filter found in the presence-absence models, particularly when only site-exclusive species are modeled (Figures 3, 7 bottom panels), is characterized by maximum percentage variance explained by environmental covariates between 20.1% and 21.2% for SCC (entire community model and the exclusive species model, respectively). A possible explanation for the low thermal signal in species occurrence is the slight warming of the ocean observed between 2017 and 2021 (Figure S7). In addition, the SACW upwelling is not always able to leave a clear thermal signature at the surface (Castelão et al., 2004). Despite this warming trend, the presence of upwelling was detected in our study as shown in Figure 11, when periods of higher SCC converge coincide with periods of stronger northerly winds (negative MWC) and lower sea surface temperature values. Emphasis on SCC may indicate that at least some species are capable of tracking environmental changes for energetic reasons. The mostly negative response of a number of species to SCC, SST and MWC seen in the β parameter plots in Figures 4 and 8 also points to the possibility of species tracking environmental conditions for different reasons. For example, the abundance of the cleaner fish Elacatinus figaro responds negatively to SCC, which is found together with Abudefduf saxatilis, that responds positively to SST. These species occur at both sheltered and exposed sites and show a positive co-occurrence at the level of location (species no. 45 and 42, respectively in Figure 6B) suggesting that the cooling of coastal water or the increase in SCC could negatively affect their abundances. It is also possible that species exclusive of the sheltered location, with high geographic ranges (given by the number of 5° × 5° grid cells occupied by each species, corresponding to approximately 550 × 550 km at the equator) like Acanthurus coeruleus (53, occurring from the USA to Brazil and Ascension Island) (Table S1) and Decapterus macarellus (133, occurring from Canada to Brazil and oceanic islands) could venture in and out (offshore) of the sheltered location to take advantage of shelter or higher productivity induced by the increase of SCC (Quimbayo et al., 2019). The niche-like effect in the presence-absence models is mainly found in the site-exclusive community and may not correspond to an actual environmental filter (Kraft et al., 2015), but more like an environmental signaling that stimulates short-distance migration between coastal and offshore habitats. This is possibly why we still find a relatively strong spatial structure of 23.2% in the global variance partitioning of the presence-absence exclusive species model (Figure 7 bottom panel). Abiotic and biotic historical contingencies that cause differences in species arrival are called priority effects, and are capable of influencing the community structure via population growth and species interaction (Fukami, 2015). The residual species associations depicted from the Ω parameter are controlled for by SST, SCC and MWC, which means that it probably indicates biotic interactions. This comes with a caveat that it could also be responding to a missing environmental covariate (i.e., species loadings depend on the environmental covariates; Ovaskainen and Abrego, 2020). It is, however, plausible to suggest that the interactions identified at the location level (Figures 6B, F, 10F) may reflect species associations induced by their response to SCC at both exposed and sheltered locations.

The shared response of species to the environment (niches) depicted from the β parameter reveals a consistently negative relation of abundance and occurrence to SCC and MWC. This reinforces the role of coastal upwelling as evidence of the joint action of migration-niche assembly. This pattern is also similar across the complete set of species found in the sheltered and exposed locations (Figures 4, 8), with a slight predominance of exclusive species found in the sheltered habitat (Figure 8). The fact that more species occur in the sheltered area than in the exposed location may indicate that the Arraial do Cabo sheltered bay is a focal habitat for the interaction between environmental variability and species occurrence. Indeed, these species are called exclusive in the sense that they are predominantly found in one location, but it does not mean that they do not potentially occur elsewhere. Although the abundance and occurrence of most species decrease with increasing SCC, those with high geographic range tend to increase their abundances with increasing SCC (Figure 5). The influence of geographic range index on all species of the community shown by the negative γ parameter in Figure 5 is represented by the intercept and can be translated as species with high range value having low abundance. If mobile species come close to the coast during upwelling events, then their abundance would increase, while those with low range will have their abundances decreased. The co-occurrence matrix (estimated Ω parameter) of residual correlations between pairs of species structured at the level of location show that a group of five species (species no. 14, and 16 to 19 in Figure 6B) increase (decrease) in abundance more often (than by chance) when most of the species from Figure 4 (left panel) reduce their abundances (species no. 36 to 58 in Figure 6B), showing a negative response to SCC. A similar pattern is depicted from the presence-absence model (Figure 6F) but, in this case, SCC contributes only 7.1% to the global variance of abundance (Figure 3 top panel). For the presence-absence model with higher SCC influence (20.1%) traits did not contribute to species response to the environment. No clear causal relation between niche and traits was found, but it is possible that some relevant phylogenetic-structured trait is missing from the analysis, which could lead to the observed phylogenetic structuring of the community. It is true that phylogenetic closeness does not imply similarity of realized niches, but in terms of abundance, spatio-temporal structuring is more important than niche in these reef fish communities from Arraial do Cabo.

The recurrence of natural environmental changes provides new opportunities for more vagile species to take advantage of exploiting available food supply. For example, the abundance of pelagic and coastal carangids with wide depth ranges (up to 400 m depth) like Decapterus macarellus and Caranx latus increase with SCC (Figure 4) in the sheltered and exposed locations. This is not restricted to mobile pelagic species, as some reef dwellers like Chaetodon sedentarius and Heteropriacanthus cruentatus also correlate positively with SCC. When upwelling conditions fade and SST increases, so does the abundance of typical shallow water reef dwellers like Abudefduf saxatilis and Cryptotomus roseus. Other species such as holocentrids (Holocentrus adscencionis and Sargocentron bullisi) are common to both sheltered and exposed locations and show a positive response to SST. Such versatility of habitat choices and environmental response of species found in Arraial do Cabo suggest that niche filtering can be approached not as a deterministic process, but more like an on-off mechanism triggered by environmental cues that some species might explore to their advantage.

A non-deterministic assembly process in Arraial do Cabo could have evolved from the southwestern Atlantic (SWA) species pool where nearly half of reef fishes are generalists (48% of species) with roving mobility, being pelagic or demersal spawners that can be found in up to five different habitats (Pinheiro et al., 2018). A flexible diet and enhanced dispersal capability would have made it easier for a relatively large number of species to move around and explore transient environmental states provided by coastal upwelling. In fact, thermal affinity alone may not accurately explain changes in abundance in transitional regions with water stratification, range extension and depth distribution of species also being important (Gomes et al., 2023; Silva et al., 2023). The schematic representation shown in Figure 12 summarizes the assembly process discussed above. We argue that the history of colonization attempts by species may have led to different rounds of community assembly, with initial niches being filled according to the arrival order. Non-upwelling conditions (Figure 12A) would favor warm water species to be the first to establish in sheltered and exposed sites. Then, the onset of upwelling conditions (Figure 12B) would facilitate the establishment of cooler water species, with the sheltered location possibly offering some level of physical protection against wave action to warm-affinity species. This upwelling condition can be considered as a transient state that regularly modifies the initial niches and allows new species to immigrate, but without excluding all species that arrived during non-upwelling conditions (Figure 12, species represented by red and blue symbols). The community we analyzed probably lives through changing transient states because upwelling conditions occur every summer lasting for up to three months, hence regularly opening the opportunity for new rounds of community assembly. A possible long-term outcome of priority effects on the reef community assembly of Arraial do Cabo is the functional divergence resulting from differences in arrival order and limiting similarity (Fukami, 2015; Ford and Roberts, 2018). Early arriving species would limit the abundance of those arriving late if resource use by both overlaps. We found a low level of functional redundancy (not shown) in both the exposed and sheltered locations, with nearly one species per functional group. It has been long acknowledged that the dynamic nature of environmental change across the species pool may cause range shifts and facilitate the coexistence of different species as a result of niche divergence or reduced niche overlap (Mittelbach and Schemske, 2015). Many marine assemblages show a positive relation between species and functional diversity, and reef fishes submitted to natural disturbances can attain low levels of functional redundancy (Micheli and Halpern, 2005).

Figure 12 Schematic representation of the assembly process of local communities assuming a priority effect acting on the Arraial do Cabo rocky embayment without upwelling (A) and with the influence of coastal upwelling (B), as a transient state that regularly modifies the initial niches and allows new species to immigrate. Symbols represent species belonging to different functional groups, red for warm water, blue for cold water and yellow represents typical marine pelagic species associated with upwelling conditions. Black arrows indicate the direction of one or more immigration attempts from the regional species pool. Green arrows depict cycles of community assembly, with occurrence and abundance varying after each cycle. The Cabo Frio island is represented by the black rectangle.

5 Conclusions

We devised a process-oriented study aimed at determining the role of coastal upwelling in the assembly of reef fish communities using a hierarchical analysis with Bayesian inference. The overarching finding is that communities of reef fish species living in sheltered and exposed habitats are assembled by the joint action of migration-niche mechanisms. The spatio-temporal structure controls abundance models, while environmental filtering better explains reef fish occurrence, with a predominant influence of upwelling-induced chlorophyll-a concentration. These communities are located in the coastal central east-southeast Atlantic where the regional pool of reef fish species is characterized as transitional between tropical and subtropical, centered at Arraial do Cabo. We found that environmental filtering and dispersal combine to drive community assembly. This is a hotspot capable of supporting novel communities across different environmental cycles that result from evolutionary habitat partitioning among similar species (Fairclough, 2021). Changes in connectivity between communities and the regional pool mediated by seasonal upwelling conditions favor the priority effect and probably contribute to increase species richness by enhancing coexistence (Palamara et al., 2023). We contend that the strong spatio-temporal structure of community assembly seen in Arraial do Cabo may have evolved from a stepping-stone driver of community structure along the SWA (Pinheiro et al., 2018), promoted by the large diversity of habitats and environmental conditions.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding author.

Author contributions

EC: conceptualization (equal), formal analysis (lead), investigation (lead), methodology (equal), writing – original draft (lead). LG: investigation (supporting), methodology (equal), formal analysis (supporting), review and editing (supporting). GG: investigation (equal), methodology (equal), formal analysis (supporting), review and editing (supporting). CC: data curation (lead), investigation (supporting), review and editing (supporting). TM: data curation (lead), investigation (supporting), review and editing (supporting). CF: conceptualization (equal), formal analysis (supporting), investigation (supporting), data curation (lead), review and editing (supporting). DG: Conceptualization (lead), investigation (equal), formal analysis (equal), methodology (supporting);, supervision (lead), writing – original draft (lead), review and editing (equal). All authors contributed to the article and approved the submitted version.

Funding

The author(s) declare financial support was received for the research, authorship, and/or publication of this article. EAAC was supported by a National Council for Scientific and Technological Development (CNPq) fellowship (n° 131353/2021-7) and received additional support from the Coordination for the Improvement of Level Personnel Superior – Brazil (CAPES) – Finance Code 001. TCM thanks CNPq for a postdoctoral fellowship (#102450/2022-6). CELF is gratefully acknowledged for continued support from FAPERJ and CNPq. Fieldwork was supported by the Projeto Costão Rochoso and funded by FUNBIO. The publication of the article was supported by the Brazilian Space Agency.

Acknowledgments

We thank NASA Making Earth System Data Records for Use in Research Environments (MEaSUREs) Program for providing SST data, the E.U. Copernicus Marine Service Information for providing SCC data and the Copernicus Climate Change Service (C3S) Climate Data Store (CDS) for the distribution of wind reanalysis data.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmars.2023.1272323/full#supplementary-material

References