The Brazilian coast extends over ca. 7,490 km, from 4° N to 34° S (Figure 1), with an enormous diversity of climatic, geomorphologic, oceanographic, and ecologic features. As the ecological dynamics of classic rocky shores and flat intertidal reefs are distinct (Mahon and Mahon, 1994; Macieira and Joyeux, 2011; Andrades et al., 2018), the present study surveyed sites between 0°S and 21°S, which encompasses the transition of intertidal flat reefs to rocky shore (slope is present).

The northern coast has a humid tropical climate (Köppen Am type), while the northeast region is divided in two main climates, a tropical coast recording dry summer (semi-arid) (Köppen As type) and another section recording humid tropical type (Köppen Af and Aw types), and subtropical climate dominating southern coast (Köppen Cfa and Cfb types) (Bernardes et al., 2012; Alvares et al., 2013). The Brazilian coast can be divided into four major geological sectors (Bernardes et al., 2012), the Northern Quaternary or Amazon plain (4° N to 3° S), Eastern Tertiary (3° to 20° S), Southeastern Granitic (20° to 29° S), and Southern Quaternary (29° to 34° S). The Brazilian coast is dominated by semi-diurnal tides with the greatest tidal range in the North region (>4 m), declining southward to less than 1 m (Dominguez, 2009). Three principal surface currents flow across the Brazilian Continental Shelf: the oligotrophic South Equatorial Current off the northeast coast (12–15° S), the North Brazil Current, which meets the Amazon plume off the northern coast, and the southward meandering Brazil Current off the eastern, southeastern, and southern coasts as far as its confluence with the Malvinas current (Peterson and Stramma, 1991; Campos et al., 1996). Based on this set of features, we divided the tropical intertidal Brazilian coast (0–21° S) into four major subprovinces (Table 1 and Figure 1):

A total of 19 sites (Figure 1), intertidal flat reefs, were sampled along ca. 4,900 km of the Brazilian coast. The distance between each site ranged from 10 to 200 km. Field sampling was conducted in a single season (March–June/2012) in order to minimize the effect of seasonal shifts in environmental data among sites. Eight tidepools (N = 152) were randomly chosen at each site to describe its fish assemblage. Sampling was always carried out during diurnal low tides in isolated tidepools (i.e., without connection to either the sea or other tidepools).

Multivariate analyses were performed in PRIMER-E (following Clarke and Gorley, 2006). Firstly, the full dataset containing all 25 environmental variables (Supplementary Table S1) was tested for collinearity using a Draftsman plot and Spearman correlation matrix. Where pairs of variables had a correlation coefficient of 0.9 or more, one of the correlated variables was excluded from the analysis (Clarke and Ainsworth, 1993). Two pairs of variables (tidepool area vs. planar area, and perimeter of the rock formation vs. Frac) were highly correlated (r² > 0.9), so one of each pair (tidepool area and perimeter of the rock formation) was omitted from the analysis, to avoid collinearity. Eight variables required transformation (e.g., log- or fourth root-transformation) to compensate for skew.

All the selected variables were standardized and a Principal Components Analysis (PCA) was used to evaluate their variation among sites and subprovinces. A permutational multivariate analysis was also used to test for differences in tidepool fish composition and environmental characteristics in the different “subprovinces” (4 fixed levels: Amazon estuary, Northern mangrove, Northeastern semiarid and Tropical warm), and “site,” nested within “subprovinces” factor, 19 random levels (see PERMANOVA; Anderson et al., 2008). Statistical significance was tested using 9,999 permutations of the residuals in a reduced model (according to Freedman and Lane, 1983) and the Type III (partial) sums of squares (see Anderson et al., 2008). In these analyses, the PERMANOVA was run on a Bray-Curtis similarity matrix, calculated from the fourth-root transformed biological data and a Euclidean similarity matrix for the environmental data (see Anderson, 2001; McArdle and Anderson, 2001). When the results were significant, post-hoc pair-wise contrasts were performed to establish where the differences in assemblage structure and environmental variables were found. In addition to the PERMANOVA, we conducted a Permutational Analysis of Multivariate Dispersions (PERMDISP) on the same Bray–Curtis matrix to test for differences in the multivariate dispersions of the “subprovinces” and “sites.”

The Similarity Percentages – species contributions (SIMPER; Clarke and Warwick, 2001) procedure was applied to reveal which fish species contributed most to the multivariate community patterns within the “subprovinces.” A constrained ordination technique, the Canonical Analysis of Principal Coordinates (CAP; Anderson and Robinson, 2003; Anderson and Willis, 2012), was performed to test the null hypothesis, given that no differences were found in the species composition among the “subprovinces” established a priori. A Spearman correlation of over 0.7 was used as an arbitrary limit for the display of the potential correlations between fish densities by feeding category and the canonical axes.

A Permutational Distance-based Linear Modeling (PERMDISP; Anderson, 2001; McArdle and Anderson, 2001) was used to test the relationship between fish assemblages and environmental variables. This procedure consists of a multivariate multiple regression method in which the ordination axes from a resemblance matrix of the response dataset (species composition) are regressed on a matrix of explanatory variables (the 23 environmental variables). A stepwise DISTLM routine was also run with 9,999 permutations selected using the Bayesian Information Criterion, BIC (Kuha, 2004). A distance-based redundancy analysis (dbRDA) was used as a constrained ordination technique to visualize the results of the DISTLM (Anderson et al., 2008).

The PCA revealed a clear segregation of the four proposed intertidal subprovinces (Figure 2). Five of the components of the environmental data modeled approximately 61% of the total variation. The first principal component (PC1) explained 21.8% of the total variation, and clearly separated the equatorial (AE and NM) from the tropical (NS and TW) subprovinces, based on salinity, rainfall, and the environmental features of the tidepools. Except for rainfall and rock coverage, the environmental data were negatively correlated with PC2, explaining 13.8% of the total variation. Similarly, the PERMANOVA indicated significant differences among sites (p = 0.0001) and subprovinces (p = 0.0001; Supplementary Table S2), where “site” was the main effect, explaining 37.2% of the variability in the data, followed by “subprovince,” which explained 30.6%.

A total of 5,113 fish, representing 67 taxa belonging to 27 families was sampled from the 152 tidepools studied (Supplementary Table S3). No fish were recorded in three tidepools. Scartella cristata was the most abundant species, accounting for 33.5% of all the individuals collected (1,309 ind.), followed by Bathygobius soporator, with 21.4% (899 ind.), and Bathygobius geminatus with 9.3% (388 ind.). The most representative families were Blenniidae (38.8%) and Gobiidae (14.7%), and only three families (Blenniidae, Gobiidae, and Batrachoididae) were recorded in all subprovinces (Supplementary Table S4). Species richness increased gradually from the AE to the TW subprovince, although a larger number of local endemics and exclusive species were recorded in the AE in comparison with the other subprovinces (Supplementary Table S5). Four species (B. soporator, B. geminatus, Labrisomus nuchipinnis, and S. cristata) were broadly distributed in the intertidal zone of the Brazilian coast, occurring in 17, 11 and 10 (both L. nuchipinnis and S. cristata) sites, respectively. By contrast, 28 taxa were restricted to a single site. Two species were endemic to subprovince AE (Gobiosoma sp. and Eleotris sp.), and one to NS (Paraclinus sp.).

The PERMANOVA explained 32.5% of the variation at the local scale (“site”) and 22.5% at the regional scale (“subprovince”), demonstrating the strong influence of local-scale factors shaping fish assemblages along the tropical Brazilian coast (p-values < 0.001; Supplementary Table S6). The PERMDISP analysis also revealed that the dispersion of the fish assemblages in tidepools varied significantly (F = 4.3169, p = 0.0003) across sites. The results of the SIMPER analysis indicated that B. soporator is typical of all four subprovinces, contributing most to similarity with the AE (94.79%). The most representative fish species in NM subprovince were Gobiesox barbatulus, Omobranchus punctatus, Barbulifer enigmaticus, Bathygobius soporator, Mugil curema, Genyatremus luteus, and Thalassophryne nattereri, whereas Scartella cristata, B. soporator, Abudefduf saxatilis, Lutjanus alexandrei, and Bathygobius geminatus were typical of the NS province. The TW subprovince was characterized by B. geminatus, S. cristata, Labrisomus nuchipinnis, B. soporator, Stegastes fuscus, and Ctenogobius boleosoma (Figure 3). Mean dissimilarity between subprovinces was 88.1, ranging from 81.8 for AE vs. NS to 96.2 for NM vs. TW.

The Canonical Analysis of Principal coordinates (CAP; Figure 4) revealed a clear spatial trophic structure, with the fish taxa grouped in the four subprovinces (δ²₁ = 0.80, p = 0.0001). Overall, the leave-one-out allocation success was 82.5% (i.e., misclassification error of only 17.5%) for the “subprovince” factor. The CAP axis scores indicated higher herbivore densities associated with the TW and NS. Omnivorous fish were correlated with the NM and carnivores with the AE (Table 2).

Distance based Linear Model (DistLM) identified water salinity explaining the largest proportion (11.3%) of the variation in the fish assemblages. Similarly, 20 of the other 23 predictive environmental variables tested were also significant (p < 0.05) in the marginal tests, with the only non-significant results being recorded for rugosity, temperature, and percentage sand cover (Supplementary Table S7). The optimal model solution consisted of the combination of seven environmental variables: salinity, tidal range, the shape of the rocky shore formation (SRF), percentage algal cover, distance to the subtidal area (DSA), latitude, and rainfall (Supplementary Table S7 and Figure 5). The results plotted in the RDA ordination indicated the segregation of AE and NM from the TW and NS samples (Figure 5). The most important variables contributing to the first axis of the dbRDA plot, explaining 44.4% of the fitted variation, were related to the DSA, salinity, rainfall, algal cover, and SRF. On the second axis, which explained 19.4% of the fitted variation, tidal range and latitude were the principal environmental variables.

Subprovinces evidenced here clearly shape intertidal fish communities. The main drivers acting on the spatial distribution of the fish were low salinity and high tidal ranges in the northernmost subprovinces (AE and NM), and high algal cover (NS and TW) and latitude (TW) toward the south. While the substrate has been identified as an important variable structuring intertidal fish communities (Arakaki et al., 2014), local-scale factors (the mouth of the Amazon, mangroves, and low-salinity estuarine waters) were the main factors defining the coastal subprovinces (e.g., AE and NM). However, regional-scale latitudinal trends were also key factors determining the configuration of the Tropical West subprovince. Local and regional drivers also act synergistically to shape subprovinces in intertidal and subtidal reefs worldwide (Rivadeneira et al., 2002; Floeter et al., 2008; Travers et al., 2010; Arakaki et al., 2014).

The spatial segregation of the geomorphological and oceanographic characteristics observed in the present study, together with the positive correlations between physical tidepool features (e.g., salinity, DSA, algal cover) and community variables further support the subprovinces proposed here. In this case, the established Brazilian Province (Briggs and Bowen, 2012; Pinheiro et al., 2018) can be divided into at least four distinct coastal subprovinces, based on the tropical intertidal seascapes. No strong correlations were found between physical tidepool features (e.g., pool volume and rugosity) and community variables, however, although a number of other studies have identified the area and volume of tidepools among the principal factors that determine fish density and diversity (Prochazka and Griffiths, 1992; Mahon and Mahon, 1994; Prochazka et al., 1999; Meager et al., 2005). The analytical approach we adopted allowed us to incorporate local and regional variables in a step-wise procedure, which indicated that species distributions and its composition may be the result of responses to the predominant environmental gradient. This contrasts with local-scale studies, which point to factors such as competition, foraging, and predator avoidance strategies, as principal drivers of fish assemblage structuring, which is the typical perspective in other marine environments (Martino and Able, 2003). In particular, intertidal reefs in the low latitude subprovinces (AE and NM) represent relatively harsh environments for marine fish, given their high rainfall and reduced salinity, emphasizing the influence of estuarine features (see Crabtree and Dean, 1982). In northern Brazil, the influence of the Amazon and other estuaries on coastal ecosystems is fundamental to the distribution and organization of marine species (Barletta et al., 2003; Rocha, 2003; Giarrizzo et al., 2010). For example, many of the intertidal ecotones of these two subprovinces are located within estuaries (Giarrizzo and Krumme, 2008; Oliveira et al., 2016), which not only determines the presence of a specific intertidal fish community, but also provides habitats for relatively large populations in comparison with the subprovinces located at higher latitudes (see Supplementary Table S3). Crabtree and Dean (1982) also found that estuarine tidepools harbor high fish densities. However, the stability of the environmental variables in the NS and TW subprovinces is demonstrated by the closer grouping of their samples in PCA ordination in comparison with the AE and NM, highlighting the heterogeneity of the intertidal habitats in these subprovinces, which also likely affects population density. In fact, Oliveira et al. (2016) observed three distinct intertidal habitats (rocky shore, mangrove, and salt marsh) within one coastal bay in the AE subprovince.

In general, the paucity of consistent significant effects of the local and structural environmental variables on the abundance and distribution of intertidal fishes is surprising given the apparent prevalence and importance micro-habitat arrangements in other habitats, terrestrial or marine (Pickett and White, 1985; Hector and Bagchi, 2007). However, the approach adopted in the present study provides alternative insights through a broader analysis of the factors regulating the biogeographical processes that mold intertidal reef fish communities, that is, the large-scale environmental factors that influence the distribution of fish most intensely. Further research incorporating data on the life cycle of the species that inhabit tidepools (i.e., residents and non-residents) will contribute to a better understanding of the micro- and meso-scale regulatory factors, which will be necessary to determine the exact nature of the processes that determine fish distribution patterns.

The plume of the Amazon River forms a natural biogeographic barrier to many marine fish species, limiting the passage between the Caribbean and Brazilian provinces (Gilbert, 1973; Collette and Rützier, 1977; Rocha, 2003; Floeter et al., 2008). However, this barrier is considered to be semi-porous, allowing many marine species to cross when permitted by specific environmental conditions (e.g., Amazon run-off, sedimentation rates) and the ecological traits of the species, e.g., reproductive mode and body size (Joyeux et al., 2001; Rocha, 2003; Floeter et al., 2008; Pinheiro et al., 2018). Highly mobile species of large body size, with pelagic spawning (e.g., the Carangidae and Serranidae), and with capacity to tolerate a diversity of conditions will generally have a much greater chance of crossing the Amazon plume than small-sized species with restrict ranges, demersal spawning, and parental care, e.g., the Gobiidae (Luiz et al., 2012; Pinheiro et al., 2018). In this context, the dispersal capacity of species resident in tidepools is limited by their ecological traits, which include small body size, as well as their reduced swimming mobility, demersal spawning, and cryptic behavior (Barlow, 1999; DeMartini, 1999; Macieira et al., 2015).

In the present study, the AE and NM subprovinces limit the presence of many resident species of both the Caribbean and the NS–TW subprovinces. This spatial pattern is expressed through (i) number of families/species that are permanent residents in the AE–NM vs. the NS–TW subprovinces (Supplementary Tables S3, S4); (ii) the absence of permanent marine resident species in the AE subprovince (Supplementary Table S4), and (iii) the exclusive and endemic species of the AE subprovince (Supplementary Table S5). Considering the two adjacent pairs of subprovinces, the AE–NM and NS–TW, only six species – Barbulifer enigmaticus, Bathygobius soporator, Ctenogobius boleosoma, Gobiesox barbatulus, Microdesmus bahianus, and Omobranchus punctatus – occur in both pairs of subprovinces, and as O. punctuatus is an invasive species, it does not represent a true biogeographic pattern. Low species richness to the north (i.e., the Amazon Estuary; Supplementary Table S5) may be related to the same environmental filters (i.e., the Amazon plume and sedimentation rates, with low salinity and turbid waters) that limit the presence of typical marine species.

The environmental filters of the subprovinces are best exemplified by the distribution of the resident tidepool frillfin gobies, B. soporator and B. geminatus, which each present distinct habitat requirements and niche. Bathygobius soporator was distributed in all subprovinces, but at particularly high mean densities at the AE sites (225.9 individuals per m³), which is likely related to its preference for tidepools with estuarine features (Macieira and Joyeux, 2011; Oliveira et al., 2016). By contrast, B. geminatus, which is distributed between Florida and the TW subprovince, typically inhabits classical marine tidepools (Macieira and Joyeux, 2011; Andrades et al., 2018), and is absent from the AE and NM sites, while increasing in density toward the south. The subprovincial endemism of tidepool fishes in the AE also highlights the local features that act as barriers. Van Tassell et al. (2015) reported two resident species of Gobiosoma (denominated here as Gobiosoma spp.) endemic to intertidal areas of the Amazon estuary, as we observed in Eleotris (Guimarães-Costa et al., 2016), which has two endemic intertidal resident species (currently undescribed) in the AE and NM subprovinces, one of which (from the NM) was not record in the present study. This high degree of endemism may be accounted for by sympatric speciation events, through which the local drivers (e.g., environmental filters) of a subprovince may promote scenarios favorable to speciation at a local-scale (Doebeli and Dieckmann, 2003; Rocha et al., 2005; Bowen et al., 2013).

One other process that should be considered is ephemeral ecological speciation (Cutter and Gray, 2016; Pinheiro et al., 2017), given the observation of one new Gobiosoma species closely related to Gobiosoma alfiei (Van Tassell et al., 2015), and the other two new, cryptic Eleotris species (Guimarães-Costa et al., 2016), which probably evolved in the context of dynamic changes in the characteristics of the Amazon barrier (run-off, sedimentation rates) over geological time (Joyeux et al., 2001; Rocha, 2003). This incipient speciation would be the result of intermittent dispersal and colonization events (Cutter and Gray, 2016; Pinheiro et al., 2017). Pinheiro et al. (2017) concluded that the fusion process between nascent species and source area populations explains why most marine species are not differentiated without strong genetic selection. This further reinforces the role of the Amazon River as a biogeographic barrier for several tidepool species, which may be related to the high degree of endemism recorded in the AE subprovince. A second barrier highlighted by Pinheiro et al. (2018) is the South Equatorial Current bifurcation, where the Brazil Current flows to the south and the North Brazil Current flows to the north-west, between 10 and 14°S (Rodrigues et al., 2007), which may also constrain population movements and dispersal. Although, the São Francisco River run-off cited as possible barrier in this same study, it’s not helps to define our intertidal coastal subprovinces.

In general, herbivores are expected to prevail at lower latitudes and carnivores at higher ones in the marine realm, due to the species’ physiological constraints and the relative availability of primary resources in the tropical region (Choat, 1991; Ferreira et al., 2004; Floeter et al., 2005). In the present study, however, the observed trophic patterns did not correspond exactly to latitudinal trends, although we did identify local and regional patterns of trophic group distribution, with carnivores (AE) and omnivores (NM) dominating the subprovinces at lower latitudes. This can be at partly accounted for by the turbid waters and high river discharge of this region, which restrict algal development, reducing the availability of resources for herbivores in the AE and NM subprovinces. The carnivore and omnivore trophic structures of assemblages were also influenced strongly by the high densities of B. soporator and O. punctatus in the AE and NM subprovinces, respectively. The presence of these species and other estuarine-related species, such as Mugil curema, Genyatremus luteus, and Thalassophryne nattereri, further reinforces the role of local-scale factors (e.g., river discharge) in the structuring of tidepool assemblages in these subprovinces. At the nine TW sites, herbivores were the predominant group structuring the intertidal communities, which is expected in shallow reef habitats of tropical regions (Steneck et al., 2017).

From a conservation perspective, the subprovince profiles observed in the present study represent an essential step toward the more effective management of Brazilian coastal zones. Overall, the Amazon Estuary subprovince appeared to be the most vulnerable, due to its endemism and environmental particularities. Conservation priority actions in the Amazon coastal region were also suggested by Vilar et al. (2017) for estuarine fishes. Currently, oil drilling and coastal urban development are the principal threats to the marine biota of the Brazilian Province (Buruaem et al., 2012; Moura et al., 2016). Urban development, in particular, not only causes habitat losses, but also favors biological invasions through the discharge of ballast water by international shipping (Carlton, 1996; Wonham et al., 2001), which invariably occurs close to intertidal areas. All the Brazilian subprovinces identified here have already undergone invasions from species such as the sun corals, Tubastrea spp. (Creed, 2006), and the muzzled blenny Omobranchus punctatus (Lasso-Alcalá et al., 2011), and ongoing invasions, such as that of the lionfish, Pterois spp., are even more preoccupying (Luiz et al., 2013). We recorded O. punctatus at seven sites. Potential niche overlap and competition with native intertidal species seems to be unavoidable, and will be harmful, given that the ecological traits and habitat requirements of this species are highly similar to those of endemic species of the Brazilian Province, such as Entomacrodus vomerinus, Ophioblennius trinitatis, and Gobiosoma alfiei.

Currently, while Marine Protected Areas contemplate 26% of Brazilian marine territory, less than 2% of these areas include coastal zones and intertidal seascapes. This leaves these coastal zones extremely vulnerable to countless impacts (oil and gas exploration, tourism, mariculture, fisheries, and port activities), which could be minimized through more participative management, firm intention and concerted actions by environmental agencies (Jablonski and Filet, 2008). In this context, we would recommend an integrated approach, involving environmental agencies (at federal, state, and municipal levels), research institutions (universities and research centers), and general society (community leaders and NGOs), for the establishment of priorities for the use and conservation of coastal zones. While it will be impossible to eliminate conflicts altogether, collaborative networks have led the way to effective management and conservation in large-scale territories under multiple jurisdictions (Wyborn and Bixler, 2013).