The principal questions in community ecology are: how are communities structured; what mechanisms direct this structure; and do these mechanisms change over spatial, temporal, or abiotic gradients. The species diversity within ecological communities greatly influences the functioning of an ecosystem (Edwards et al., 2010; Belley and Snelgrove, 2016; Brisson et al., 2020). Understanding the direction and strength of the drivers of these observed ecological patterns is especially important as it pertains to predicted climatic changes. Ecological community patterns are enmeshed in the interplay between abiotic factors and interspecies interactions, and untangling this labyrinth of linkages is essential to unlocking the latitudinal gradient first described by Darwin and Wallace (Darwin, 1859; Darwin and Wallace, 1958; Wallace, 1905; Willig et al., 2003; Scheiner and Willig, 2005). With so many factors influencing ecological communities, assessing which of these factors have large effects on community structure is at the core of modeling and understanding ecosystems.

One general pattern of ecological biodiversity is that of higher biodiversity in tropical latitudes and declining biodiversity with increasing latitude (Hillebrand, 2004; Yasuhara et al., 2012; Parravicini et al., 2013; Tittensor et al., 2010). One hypothesis explaining the latitudinal gradient is the species richness–energy hypothesis, which states that there is more energy—in the form of solar radiation—near the equator. Moreover, this has been proven to be a good measure of biodiversity along the latitudinal gradient. Central to this hypothesis is the biological fact that biochemical activities, in general, increase with increased temperature (Aquino, 1968), which results in a faster rate of evolutionary change, resulting in higher biodiversity (Rohde, 1992; Taylor and Gaines, 1999; Willig et al., 2003; Mittelbach et al., 2007; Tittensor et al., 2010). This hypothesis has been supported by many studies (Macpherson, 2002), in which the mean annual temperature has been proven to be a good measure of energy in the system.

Predation has been put forth as another of the driving forces of trait divergence in the biotic interaction hypothesis (Pianka, 1966; Rohde, 1992). The biotic interaction hypothesis states that, in the tropics, species interactions are stronger and comprise a greater portion of evolutionary selection, which results in traits evolving faster (Willig et al., 2003; Mittelbach et al., 2007; Brown, 2014). Of the interactions that help explain community development, a number of studies have assessed the effects of predation across a latitudinal gradient (Reynolds et al., 2018). Top-down pressure from predation has a strong influence on the community structure, especially in marine environments (Edwards et al., 2010; Sheppard-Brennand et al., 2017; Reynolds et al., 2018).

Current ecological theory states that predation decreases with increasing latitude (Sheppard-Brennand et al., 2017). Predation pressure affecting community diversity differently across a latitudinal gradient has been recorded in shallow water systems worldwide (Pianka, 1966; Freestone et al., 2011; Sheppard-Brennand et al., 2017). The marine benthic community is structured by two important biological processes: colonization by larvae (Schwamborn and Bonecker, 1996; Leal et al., 2019; Rossi et al., 2018) and post-settlement survival (Thorson, 1966; Sutherland, 1974; Osman and Whitlatch, 1995; Freestone et al., 2011; Yakovis and Artemieva, 2017). Predation plays a prominent role in post-settlement survival (Gosselin and Qian, 1997; Kulp and Peterson, 2016). The predation hypothesis states that predation pressure lowers the level of competition, thus allowing for higher biodiversity through the coexistence of more prey species (Pianka, 1966).

Another strong influence of community structure is the bottom-up impact of foundation species (Cusseddu et al., 2016; Sala et al., 2008). As the old adage in architecture states, “structure dictates function.” The organization of an ecosystem is based on the structure of the foundation species as this dictates diversity, resilience, food web complexity, and productivity (Dayton, 1972; Bruno and Bertness, 2001; Yakovis and Artemieva, 2017). The interactions between co-occurring foundation species can have significant influences on predator–prey interactions (Ware et al., 2019), the trophic structure, and other ecological functions (Altieri et al., 2007; Angelini and Silliman, 2014; Yakovis and Artemieva, 2017). Much of the ecological effects of foundation species are due to their habitat-modifying capabilities (O'Brien et al., 2018; Vozzo and Bishop, 2019; Yakovis and Artemieva, 2019). The interstitial spatial structure of foundation species affects predator–prey interactions (Ware et al., 2019). In Florida’s mangrove, ecosystem shifts in the secondary foundation species—foundation species that reside on mangrove prop roots—had a strong influence on the epifaunal biodiversity on mangrove prop roots (Aquino-Thomas and Proffitt, 2025).

Red mangroves are a foundation species and ecological engineer that influences the microclimate, sedimentation, and subsidence rates and provides many ecological functions and services (Feller et al., 2010; Aquino-Thomas and Proffitt, 2014; Doughty et al., 2016; Aquino-Thomas, 2020). Collectively, mangrove species will be referred to here as the primary foundation species. Encrusting communities of invertebrates and algae form on the prop roots of red mangroves and can serve as secondary foundation species by influencing the diversity of many small motile and sessile species. By modifying the physical environment of the communities, these species create microhabitats with varying levels of light, water flow, and nutrient availability (Bishop et al., 2012; Angelini and Silliman, 2014; Thomsen et al., 2018; Ellison, 2019). This heterogeneity supports a greater variety of other organisms, including filter feeders, photosynthetic species, and detritivores (Sebens, 1991; Gallucci et al., 2020; Thomsen et al., 2022). This combination of foundation species, and the associated community, provides a unique opportunity to assess the effects of predation in a two-foundation species system along an environmental gradient. The space-for-time substitution (Pickett, 1989) along the latitudinal gradient will affect the dominance of the foundation species (Aquino-Thomas and Proffitt, 2025) and the influence of predators, which may vary with the differences in interstitial spaces that occur with different dominant taxa comprising the secondary foundation species. The supplemental refuge effects of having multiple foundation species is greatly influenced by the extent of the functional redundancy between foundation species (Vozzo and Bishop, 2019). This feature of foundation species creating refuge from predation has been shown to be important in varied ecosystems, such as mobile sea urchins (Altieri and Witman, 2014), Australian mangroves (Vozzo and Bishop, 2019), macroalgal ecosystems (O'Brien et al., 2018; Ware et al., 2019), red mangroves (Schutte and Byers, 2017), kelp forests (Efird and Konar, 2014), coral reefs (Dunn et al., 2017), oyster reefs (McAfee and Bishop, 2019), and even dead oysters (Tomatsuri and Kon, 2017).

Predation and secondary foundation species play crucial roles in structuring sessile mangrove prop root communities. However, their relative importance and their interactions across biogeographic gradients remain poorly understood. While previous studies have examined predation effects or foundation species’ impacts individually, only a few have investigated how these factors interact across large spatial scales to shape community assembly. This study aimed to fill this knowledge gap by assessing the combined effects of predation and secondary foundation species on mangrove epibiont communities along a latitudinal gradient in Florida.

The overall objective of this field experiment was to assess the relative influence of predation and secondary foundation species on the mangrove prop root species diversity along the latitudinal gradient. By deploying predator exclusion treatments at four sites spanning from tropical to temperate zones and by monitoring community development over 6 months, we aimed to elucidate how top-down and bottom-up forces interact to structure these communities across biogeographic regions. We hypothesize that: a) the effects of predation on community structure will shift from positive in tropical sites to negative in temperate sites, in line with the predation hypothesis; b) secondary foundation species will generally increase biodiversity when space is not limiting, but their effects will vary depending on the predation pressure and latitude; c) the interaction between predation and secondary foundation species will explain more variation in the community structure than either factor alone; and d) the relative importance of predation versus facilitation by secondary foundation species in shaping communities will shift across the latitudinal gradient.

Initial deployments of the treatments at each site were completed in June 2015. Throughout the field experiment, a total of 62 species were recorded. From south to north, 17 species were observed at Curry Hammock, 37 at Arch Creek, 42 at Rutherford, and 27 at HBOI. In the southern half of the experimental zone, sponge and ascidian species were more prevalent, regardless of treatment type. Tropical species such as the sponge (Dysidea etheria) and several Botrylloides ascidian species were absent at the northernmost site, while temperate species such as the eastern oyster (Crassostrea virginica) were not found at the southernmost site.

Species richness was assessed during the three sampling periods. The first data collection point (approximately 1 month after deployment) identified 33 species. The subsequent data collection, which was restricted to Rutherford and HBOI, yielded 30 species. The last data collection point (November and December) recorded the highest species richness, with 51 species. The total species richness over the length of the entire experiment was further examined by treatment within each site (Table 1). At all sites, an increase in species richness was observed over the course of the study: Curry Hammock, from 3 to 17; Arch Creek, from 25 to 32; James Rutherford, from 14 to 30; and HBOI, from 9 to 17. The total species richness at the end of the study was further examined by treatment within each site (Table 2).

Predation produced opposite effects at either end of the extremes of the latitudinal gradient: increasing species richness (Figure 3) in the north and less clearly reducing species richness in the south. Species richness was initially highest in Arch Creek, the transition zone, and remained the highest over the course of the study. Interestingly, at this site, initially, the exclosures that prevented access of at least some predators had the highest species richness; however, by the end of the study, they had the lowest species richness at the site. The southernmost site, Curry Hammock, had the lowest initial species richness and a weak current flow at the site, as it was situated snugly between two keys, Long Point Key and Little Crawl Key, likely contributing to the results more than any experimental factor.

Community diversity was also measured using Shannon–Wiener. A comparison of the treatment types over the latitudinal gradient found that the open treatment type had lower diversity values in the northern half of the study when compared with its protected treatment exclosure counter type, while in the south, the predation did not have as clear of an affect. In the northern half of the study, certain soft-bodied species were not observed, unless they were in macro- and/or micro-predation exclosures (Figure 4). The Shannon–Wiener for the initial data collection visit across all sites was 2.69; for the second data collection (which did not include two sites, Curry Hammock and Arch Creek), it was 2.199. The biodiversity for the end of the experiment across all sites had a Shannon–Weiner measure of 2.402: 2.244 for Curry Hammock, 2.716 for Arch Creek, 2.171 for Rutherford, and 0.900 for HBOI. The Shannon–Weiner at the end of the study was further parceled out by treatment within each site (Figure 4). Overall, Arch Creek had the highest diversity, and abundance was more evenly distributed (Pielou’s evenness: Curry Hammock, 0.544; Arch Creek, 0.658; Rutherford, 0.526; HBOI, 0.218). On the other hand, HBOI had one or a few dominant species. However, the HBOI micro exclosure treatment type broke with the trend and had an appreciably higher species evenness (open = 0.064, partial macro = 0.132, partial micro = 0.189, macro = 0.163, and micro = 0.446). Measures of pairwise dissimilarity can have limited usefulness when the local species richness is very low in comparison to the regional species richness.

Arch Creek had one treatment type, i.e., macro exclosure, which did not see a change in its mean Shannon–Weiner value. This appears to have been in relation to the species abundance/percent cover as the species richness for the macro exclosure treatment increased over the same time span. Over the course of the study, HBOI, not including the micro exclosure treatment type, had a decrease in its Shannon–Weiner value. The final HBOI micro exclosure treatment had values that were within the range of the next site to the south, Rutherford, demonstrating the important impact of micro-predation at HBOI.

The initial colonization, from deployment to 1 month later, at Curry Hammock was extremely low. Arch Creek had the highest initial overall biodiversity, with all the treatment types exhibiting high values. At James Rutherford Park, things were a little more difficult to interpret, with three of the treatment types having high initial biodiversity and two having some of the lowest. The three treatment types—open, partial micro, and micro exclosures—had high biodiversity indices and saw a moderate reduction in biodiversity in the next collection date. Conversely, the James Rutherford macro exclosure and partial macro treatment types saw an increase in biodiversity over the course of the field experiment, with the macro exclosure treatment ending up having the highest diversity for the site. The initial diversity in colonization for the treatments at HBOI was similar for all treatment types, except for the macro exclosure treatment (Figure 4A). At the end of the study, the HBOI site saw a further decrease in biodiversity for all treatment types, except for the micro exclosure treatment, which mimicked the species richness’ increasing slope, likely due to an increase in the number of species present. The combined pressure of macro- and micro-predation in the temperate zone had significantly decreased species richness (Figure 3).

At the conclusion of the experiment, the Shannon–Weiner values were largest for Arch Creek, the transition zone (Figure 4B). There was a decrease in biodiversity heading north to Rutherford, although there was a substantial overlap in biodiversity between the two sites. At Rutherford, the macro-predation exclosure treatment was the treatment with the highest biodiversity. The last site along the latitudinal gradient had the lowest biodiversity, with four of the treatments showing a drastic drop in diversity (Figure 4). At HBOI, the micro exclosure treatment was the only treatment that did not follow the trend, instead retaining a biodiversity that was extremely close to the Rutherford micro exclosure biodiversity and higher than the Rutherford open biodiversity. When HBOI was released from both micro- and macro-predation, the Shannon–Wiener value was within the values found for the next site to the south. Again, micro-predation appears to be a factor limiting biodiversity in the temperate zone of the experiment. The biodiversity at Curry Hammock is largely an artifact of the weak water flow at the site.

PERMANOVA was conducted using Bray–Curtis dissimilarities to examine the effects of site and predation exclosure treatment on the community composition, with 999 permutations (Table 3). Before interpreting the PERMANOVA results, a test for homogeneity of multivariate dispersions (PERMDISP) was performed using the same Bray–Curtis dissimilarity matrix to ensure that the differences in the location of the group centroids were not confounded by the differences in group variability. PERMDISP revealed a significant difference in the multivariate dispersion among the four sites (Pseudo-F = 19.499, df = 3, 78, p = 0.001). Pairwise PERMDISP tests indicated that all site comparisons were statistically significant (p ≤ 0.004), except for the comparison between Arch Creek and Rutherford (p = 0.251). These significant dispersion differences among sites suggest that the interpretation of the site effects from the PERMANOVA should be made with caution, as differences in the community composition between sites could be influenced by varying levels of within-group heterogeneity. PERMDISP was conducted for the treatment factor using the same Bray–Curtis dissimilarity matrix. The results of this test revealed no significant differences in the multivariate dispersion among the five treatment groups (Pseudo-F = 0.402, df = 4, 77, p = 0.8), suggesting that the variability within the treatment groups is similar. This finding increases confidence in the observed significant effect of treatment found in the PERMANOVA. The PERMDISP for the treatment within each site revealed a significant difference in the multivariate dispersion among the four sites (Pseudo-F = 3.825, df = 19, 62, p = 0.001).

PERMANOVA revealed that both site and predation exclosure treatment significantly influenced the community composition. The effect of the sites was stronger (Pseudo-F = 12.347, df = 3, 78, R² = 0.322, p = 0.001) than that of the predation exclosure treatment (Pseudo-F = 1.569, df = 4, R² = 0.0754, p = 0.001). Pairwise comparisons among the four sites were then performed to identify specific differences. All pairwise comparisons were statistically significant at p ≤ 0.002. Site explained approximately 32.2% of the total variation in community composition. However, even accounting for the site-specific differences, treatment still had a statistically significant influence on the community structure. Further pairwise tests for the nested effect of treatment within each site were completed. The comparison of open vs. macro (Pseudo-F = 4.734, df = 7, 78, R² = 0.570, p = 0.021) and open vs. micro (Pseudo-F = 3.361, df = 7, 78, R² = 0.475, p = 0.039) were significant across all sites. Macro was also nearly significantly different from partial macro (Pseudo-F = 3.644, df = 7, 78, R² = 0.495, p = 0.061) and was significantly different from partial micro (pseudo-F = 2.78, df = 7, 78, R² = 0.481, p = 0.016). Micro was not significantly different from either partial exclosure.

The variability in the F and R² values across sites for the same treatment comparison suggests that the impact of predation exclosure treatments is not uniform across all locations. Local environmental conditions and/or other site-specific factors likely modify how treatments affect the community. A larger F-value suggests a greater difference between groups. For example, the open vs. cage comparison is strongest in Arch Creek (F = 1.66), followed by Rutherford (F = 1.42) and HBOI (F = 1.26), and was weakest (and not significant) in Curry Hammock (F = 0.75) at the end of the study. The comparison between open, macro, and micro were frequently significant in Arch Creek, Rutherford, and HBOI. Arch Creek also had the most significant treatment differences, suggesting that its communities might be more sensitive to the predation exclosure treatments compared with the other sites (Table 3).

In general, there was a trend of increasing percent coverage toward the north. Coverage started low and remained lowest at Curry Hammock, the southernmost site. Arch Creek had the highest initial colonization. The two northernmost sites, Rutherford and HBOI, both had a relatively low initial colonization and a high final percent coverage (Figure 5). The percent cover was low where the species richness was low and was high where the species richness was high, except at HBOI, which had the highest percent coverage by predation exclosure treatment type. This deviation is likely due to shelled specimens leaving behind hard structures upon death and the percent coverage data unable to separate live versus dead specimens.

Hierarchical clustering was used to assess differences in the community assemblages between the groups, treatments nested within the sites. Overall, the sites clustered along the latitudinal gradient (Figure 6). Clusters were based on the relatedness of the species composition of the treatments within the sites. The open treatment type, which allowed for full predator access, followed the expected latitudinal gradient. Overall, the predation exclosure treatment types clustered within the specific site. Only one treatment did not cluster with the other treatments within its site. Removal of all predation at HBOI, the northernmost site, resulted in an epifaunal community that was more similar to that in the site to the south than to that in treatments within the site.

The PCoA confirmed the hierarchical clustering of the treatments within the sites, with site remaining as the primary clustering criterion. The first two dimensions of the PCoA accounted for 28.5 of the goodness of fit, the amount of variation captured. PCoA, which included secondary foundation species in the analysis with the epifaunal community, revealed that Curry Hammock did not have any overlap with the other sites (Figure 7). Points closer together on the plot are more similar in community composition according to the Bray–Curtis dissimilarity matrix used. Rutherford’s ellipse overlapped with both HBOI, the site to the north, and Arch Creek, the site to the south. The standard deviation ellipse for HBOI was entirely nested within the ellipse for Rutherford and did not overlap with the ellipse from Arch Creek. This indicates that while Rutherford exhibited a broader range of community compositions, the entire observed variation of the HBOI communities was contained within that larger range. This included the screen treatment type for HBOI that fell outside its standard deviation ellipse, but still within the Rutherford ellipse. For the site to the south, Arch Creek, Rutherford’s macro exclosure clustered closely due to more shared species. The standard deviation ellipse for Rutherford was positioned such that it largely encompassed or significantly overlapped with the ellipse for Arch Creek, suggesting a high degree of similarity in the community composition between these two groups, despite their different geographical locations.

Within the Arch Creek site, distinct patterns of community dispersion were observed across treatment types (Figure 7). Samples from the open treatment type clustered relatively closely together, indicating a higher degree of similarity in their community composition. In contrast, samples from the macro treatment type showed a greater spread, suggesting increased variability among these samples. This pattern of increasing dispersion continued with the micro treatment type, which exhibited the largest spread of points, reflecting the highest level of heterogeneity in community composition within that predation exclosure treatment at Arch Creek. These differences in the within-group variability were visually supported by the sizes of the corresponding standard deviation ellipses, which expanded progressively from open to macro to micro treatments at this site (Supplementary Figure S1).

A global test for the partial db-RDA was performed with predation exclosure type as the conditional variable. The model score was significant (F = 1.93, df = 18, 59, p < 0.001), indicating that the constrained model significantly explains the variation in species composition. This species composition did not include the secondary foundation species in the model. The analysis revealed a total inertia of 59.13, with 5.03% (2.97) explained by the conditional term (predation exclosure type) and 35.19% (20.81) explained by the environmental constraints. The remaining 59.78% (35.35) represented unconstrained inertia. The eigenvalues for the constrained axes ranged from 6.16 (CAP1) down to 0.35 (CAP18), while the unconstrained axes showed eigenvalues ranging from 4.41 (MDS1) to 0.66 (MDS8). The model’s explanatory power, as indicated by the adjusted R², was 0.178, meaning that approximately 17.8% of the variation in species composition can be explained by the constrained environmental variables after accounting for the number of predictors in the model. This value suggests that the included environmental factors explain a moderate portion of the variability observed in the species data.

The constrained global model (adjusted R² model), using the Bray–Curtis matrix, had the constrained inertia decrease to 19.34% (11.43), and the unconstrained inertia in turn increased to 75.63 (44.72), which included three environmental variables (Figure 8). When added to the model, salinity explained 8.65% of the adjusted variation (adjusted R² = 0.0865). This addition was statistically significant (F = 8.29934, p < 0.001). Latitude further improved the model fit, increasing the adjusted R² to 0.1654. Its contribution was also statistically significant (F = 8.18772, p < 0.001). The final variable added, Disturbance.Docks, increased the adjusted R² to 0.1712, with a statistically significant contribution (F = 1.52661, p = 0.044). Ultimately, the model including salinity, latitude, and Disturbance.Docks explained an adjusted 17.83% of the variation in species composition (adjusted R² = 0.1783). The eigenvalues for the constrained axes were 5.86 (CAP1), 4.93 (CAP2), and 0.65 (CAP3), representing the amount of variance explained by each axis. The eigenvalues for the unconstrained axes, which represent the remaining unexplained variation, ranged from 5.60 (MDS1) to 0.72 (MDS8). Further permutation tests were conducted to evaluate the significance of each canonical axis individually, holding all previous axes as conditions. The first constrained axis (salinity) was statistically significant (F = 9.69, df = 1, 74, p < 0.001), as was the second constrained axis (latitude; F = 8.15, df = 1, 74, p < 0.001). However, the third constrained axis (Disturbance.Docks) was not found to be statistically significant (F = 1.07, df = 1, 74, p = 0.299). These results suggest that the primary gradients of the environmental variation influencing species composition are captured by the first two CAP axes.

The inclusion of secondary foundation species in the environmental matrix for the partial db-RDA was revealed to be a significant model after a permutation test for the significance of the overall constrained model was performed using 999 permutations (F = 1.93, df = 18, 59, p < 0.001), with a total inertia of 59.13 and an adjusted R² of 0286. The conditional term, treatment (predation exclosure type), accounted for 5.03% (2.97) of the total variation. The constrained environmental variables explained a substantial portion of the variation, accounting for 45.65% (26.99) of the total inertia. The eigenvalues for the 21 constrained axes ranged from 8.05 (CAP1) down to 0.33 (CAP21). The unconstrained axes showed eigenvalues ranging from 4.03 (MDS1) to 0.54 (MDS8).

The constrained global model with secondary foundation species in the environmental variables, using the Bray–Curtis matrix, included four variables in the model, each of which improved the model’s explanatory power. Addition of the secondary foundation species barnacle resulted in an adjusted R² of 0.119 (F = 11.40, p < 0.001), salinity increased the adjusted R² to 0.213 (F = 10.15, p < 0.001), the secondary foundation species sponge resulted in an adjusted R² of 0.252 (F = 4.86, p < 0.001), and the secondary foundation species oyster brought the adjusted R² to 0.279 (F = 3.80, p < 0.001) (Figure 9). This model explained 30.09% (17.79) of the total inertia (59.13). The conditional term accounted for 5.03% (2.97) of the total inertia. The remaining 64.88% (38.37) represented unconstrained inertia. The four variables, along with the conditional treatment term (predation exclosure), collectively explained an adjusted R² of 0.286, suggesting that these variables are particularly important in structuring the community patterns. The eigenvalues for the four constrained axes were 7.97 (CAP1), 6.03 (CAP2), 2.67 (CAP3), and 1.12 (CAP4). The unconstrained axes showed eigenvalues ranging from 5.45 (MDS1) to 0.60 (MDS8). Further permutation tests were conducted to evaluate the significance of each canonical axis individually, holding all previous axes as conditions. All four constrained axes were found to be statistically significant. The first constrained axis (barnacle) was highly significant (F = 15.17, df = 1, 73, p < 0.001). The second (salinity; F = 11.47, df = 1, 73, p < 0.001), the third (sponge; F = 5.08, df = 1, 73, p < 0.001), and the fourth constrained axis (oyster) also showed a significant contribution (F = 2.14, df = 1, 73, p = 0.003). These results indicate that each of the four identified environmental gradients captured distinct and significant patterns in species composition. Moreover, the addition of secondary foundation species to the environmental variables increased the ability to model the system.

The association between secondary foundation species and predation changed along the latitudinal gradient. The relationship between secondary foundation species and associated species was not uniform across latitudes, suggesting context-dependent ecological dynamics. To investigate the relationship between species richness and the abundance of secondary foundation species, Spearman’s rank correlation was run for each treatment type. Overall, the strength and direction of the relationship between species richness and secondary foundation species varied a considerable amount based on the different predation exclosure treatment types (Figure 10). A few of the pairwise comparisons were significant, which is likely a result of the low power that occurred due to replicates, specifically at Rutherford and HBOI, being lost by the final visit of the experiment. Regardless, the site-specific relationship strength and patterns between the secondary foundation species and the associated species can be observed. At the southernmost site, Curry Hammock, the two secondary foundation species present—sponge and barnacles—although not significant, almost exclusively had positive relationships with species richness. At the next site to the north, Arch Creek, the patterns regarding the relationship between secondary foundation species and species richness were more varied. Sponges appeared to have a negative relationship with species richness in the fully enclosed predation exclosure treatment type, while the treatments with reduced predation had a positive relationship between species richness and sponge cover. This pattern, for the most part, was reversed when oysters were the secondary foundation species. There did not appear to be any discernible pattern between richness and barnacle cover.

In the northern half of the experiment, the association with secondary foundation species became more varied and species-dependent. At Rutherford, sponges had a very strong positive relationship with species richness for the micro exclosure treatment type. Oysters, on the other hand, appeared to have the opposite correlation, reducing the species richness when all predators were excluded and increasing the species richness in the open and macro exclosure treatment types. Here, barnacles either had no correlation or had a negative correlation with species richness. At the northernmost site, HBOI, there was no effect of sponges or oysters on the species richness when all predators were permitted (open treatment), but a positive relationship when small and large predators (micro exclosure) were excluded. Oysters also had a positive relationship when only larger predators (macro exclosure) were excluded. Barnacles at the HBOI site typically had a negative relationship across the majority of the treatment types. The barnacle percent coverage was very high at the HBOI site; thus, colonization by other benthic species would be limited. Moreover, overall, barnacles exhibited a trend where, in the south, their correlation was positive and generally became negative the further north, albeit the correlations were insignificant. It is important to note that these are not fully matured secondary foundation species, and as such, their effectiveness is not fully realized.

This study revealed significant latitudinal patterns in the predation intensity and the secondary foundation species interactions within mangrove ecosystems. Lower latitudes exhibited stronger predation pressures and a more pronounced role of secondary foundation species in mitigating the effects of predation. The experimental results support the complex nature of the ecological interactions along the latitudinal gradient, revealing diverse effects of predation and secondary foundation species. This pattern aligns with the latitudinal biodiversity gradient theory (Schemske et al., 2009) and the predation hypothesis in which predation in the tropics enhances species richness by limiting competitive dominants and opening up niches for colonization (Pianka, 1966; Shoemaker et al., 2020).

In addition to the abiotic factors that contribute to structuring ecosystems, biological interactions such as predation, facilitation, and competition all directly influence the biodiversity of an ecosystem. Analyses of the epifaunal community (excluding secondary foundation species) indicated that species similarities are significantly affected by latitude, supporting the species richness–energy hypothesis (Osman and Whitlatch, 1978; Platnick, 1991; Gaston, 2000). In harsh environments, such as mangrove ecosystems, many species live at abiotic extremes and gain resilience from their association with a secondary foundation species, which not only provides refuge from predation but can also give a species a competitive advantage.

The differential effects of foundation species could be linked to their adaptations to local environmental conditions and their ability to modify habitats. These findings support existing knowledge on latitudinal patterns in ecological interactions while challenging simplistic views of uniform ecological processes across large spatial scales. The study highlights the situational-dependent role of secondary foundation species in mediating predation and in shaping the community composition, potentially influencing ecosystem resilience and the responses to environmental changes. The findings indicate that community structure and function are dynamic across latitudes.

The northernmost site, HBOI, experienced high mortality after recruitment, coupled with a reduction in subsequent successful recruitment for all species, highlighting the importance of protection from micro-predation. Full predation exclosures provide refuge for epifauna in their early life stages, as micro-predation can have significant impacts on the recruitment, diversity, and percent coverage of marine invertebrate communities (Osman et al., 1992; Nydam and Stachowicz, 2007; Freestone et al., 2011). However, the predation dynamics change over time, with species that escape micro-predation then becoming vulnerable to larger predators. At the HBOI site, the inclusion of these soft-bodied associated species in the full predation exclosure was the main difference in the species composition between predation treatment types. Barnacles and other calcified species, which are abundant in the north, have shorter periods of time when they are vulnerable to being preyed on by small invertebrates compared with ascidians and other non-calcified species (Osman and Whitlatch, 2004).

As the effects of predation change over the latitudinal gradient, so too do the effects of each particular secondary foundation species. Secondary foundation species exhibited varying degrees of influence on the community structure and the predation dynamics across latitudes. In the more tropical regions, these species played a crucial role in mediating the predation effects, providing refuge for prey organisms. This protective function was less pronounced in temperate zones, suggesting a latitudinal gradient in the importance of facilitative interactions.

At Arch Creek, sponges significantly enhanced the species richness in open treatments, indicating their role in facilitating biodiversity. Barnacles exhibited a shift from positive to a negative relationship with species richness along the latitudinal gradient. Oysters had a positive relationship with species richness when predation limited the species richness. Sponges exhibited a difficult to discern pattern, which appeared to be dependent on the latitudinal gradient and predator access. At the northernmost site, HBOI, sponges and oysters showed the same general positive relationship with species richness when both macro- and micro-predation were removed. These observations suggest latitudinal shifts in the key processes between predation and foundation species that contribute to community patterns. In tropical mangrove ecosystems, species facilitated by a secondary foundation species may experience positive recruitment and post-settlement survival, potentially gaining a competitive advantage. The importance of biogenic refuge from predation would be clearly visible in an ecosystem where predation rather than competition limits the species richness and diversity.

Secondary foundation species generally outcompeted other species for the limited space on the panels in the first few months of the field study. Space is often a limiting resource in marine benthic communities (Dayton, 1971). As secondary foundation species, by definition, have a high biomass within their ecosystems, this competitive advantage is expected (Thomsen et al., 2010). The panels allowed for the communities to age at the same rate, but also left the secondary foundation species competing for the same space as the other species. This competition negates some of the hypothesized positive effects of having a secondary foundation species present in an ecosystem, such as the formation of habitats and the amelioration of harsh environmental factors. It is important to note that these are not fully matured secondary foundation species, and as such, their effectiveness is not fully realized.

Superimposed over the predation effects were the facilitative effects of the secondary foundation species. Ecosystems are organized by a wide number of co-occurring effects over a large range of temporal and spatial scales. Our results revealed a latitudinal gradient in both the top-down predation control and the bottom-up influence of secondary foundation species, with species interactions varying based on the biographic location and the predation mechanisms. In regions with high predation intensity, the presence of secondary foundation species became increasingly important for the maintenance of biodiversity and ecosystem stability. This finding underscores the context-dependent nature of species interactions and their collective influence on the community structure.

The findings of this research are consistent with other studies that found species interactions changing across the latitudinal gradient (Schemske et al., 2009; Freestone and Osman, 2011; Freestone et al., 2011). The community response across the latitudinal gradient can be taxon-dependent, both for the secondary foundation species and the epifauna. Shifts in the top-down control of predation over the latitudinal gradient coincided with the shifts in the bottom-up control of the secondary foundation species. The interactions of the secondary foundation species with the associated species were dependent not only on where the ecosystem lay on the latitudinal gradient but also on the mechanism effects of predation.

Future research could explore the long-term dynamics of these communities, particularly focusing on how the roles of secondary foundation species evolve as they mature and reach size thresholds. In addition, investigating the potential impacts of climate change on these latitudinal patterns and species interactions would provide valuable insights for the prediction and management of ecosystem responses to global environmental shifts.

Integration of secondary foundation species into the modeling of ecosystems will result in better predictions with regard to how shifts in taxa and the community structure will be altered with depending on how abiotic factors may change. This research supports the theory of latitudinal gradients in species diversity and interaction strength, highlighting the complex, latitude-dependent roles of predation and facilitation. These insights emphasize the importance of considering a biogeographic context in conservation strategies, particularly in mangrove ecosystems where predation pressure and secondary foundation species significantly influence the community dynamics.