The PNSAV is a protected natural area that follows the coastline, and it covers an extension of 65,516 hectares in the central region of the Veracruz continental shelf (Diario Oficial de la Federación [DOF], 2012). The SAV includes 50 reefs (Liaño-Carrera et al., 2019) which are naturally divided in two coral reefs subsystems: north and south ( Figure 1 ). The north subsystem comprises 31 reefs and is found close to the city of Veracruz up to the municipality of Boca del Río. The south subsystem has 19 reefs, located in front of the community of Antón Lizardo, Veracruz. Both subsystems are geographically divided by the Jamapa River mouth, there are two other rivers: La Antigua River on the north and the Papaloapan River on the south ( Figure 1 ) (Krutak, 1997). The reefs in the north are located above the 40 m isobath, whereas the reefs in the south are above the 50 m isobath (Pérez-España and Vargas-Hernández, 2008). In general, the PNSAV coral reefs have an approximate length between 0.3 and 3.2 km, an area coverage from 1 to 19 km², and a body height of between 20 and 45 cm (Allende-Arandía et al., 2016).

Following the hierarchical scheme of Ortiz-Lozano et al. (2009), on a seascape scale, the PNSAV was divided into two large regions, associated with the contours of depth and the water currents. The inner shelf (IS) region is located between the 10 m and the 25 m depth contour and is defined by the presence of marine currents that move parallel to the coast with latitudinal and seasonal variations (Salas-Pérez and Granados-Barba, 2008). The outer shelf (OS) region presents depths > 25 m up to the limit of the continental shelf and includes waters free from the shelf, where the water currents exhibit only local movements, not parallel to the coast (Ray and Hyden, 1992). Therefore, overall quality of each shelf, showed different potential effects caused by the outflow of Veracruz city sewage (Castañeda-Chávez and Lango-Reynoso, 2021).

As part of the park monitoring program (Caracterización ecológica y monitoreo del Parque Nacional Sistema Arrecifal Veracruzano), visual counts were carried out between 2006 and 2021. The monitoring covered four reefs in the north subsystem: Blanquilla, Verde and Sacrificios (inner shelf), and Anegada de Adentro (outer shelf). In the south subsystem it covered three reefs: Blanca and Enmedio (inner shelf) and Santiaguillo (outer shelf). Every reef in this study forms a platform, where the reefs Santiaguillo, Blanca, and Verde show healthier characteristics, such as high coral recruitment, cover and richness, high fish biomass and richness, a low prevalence of coral diseases, and a low turf growing rate. In contrast, reefs like Enmedio and Sacrificios exhibit opposite conditions (Pérez-España et al., 2012). The conditions of recently surveyed reefs are not yet described, but none are emergent reefs (Liaño-Carrera et al., 2019) ( Figure 1 ).

In each reef, the observations and counts were made at two depths: shallow (3-5 m) and deep (10-15 m), except in Sacrifices reef, the closest to mainland, where the corals did not appear below 7 m in our sampling site. For each reef and depth, five fixed transects of 10x4 m were made parallel to the reef slope and always maintaining the same depth. The visual censuses were carried out always by the same expert (HPE) to maintain the observer error constant. All fishes, from the water column to the bottom, were identified, counted, and the total length of every individual was estimated (see Supplementary Material Figure S1 ). The periodicity of the censuses varied throughout the monitoring regime, first it was bi-monthly, then it was reduced to three times a year, and lastly to an annual sampling. Finally, the number of individuals of each identified species, their size, depth, reef, month, and year were recorded in a database.

The modal length was determined for each species, and it was used as the species length. Length-weight relationship was used to calculate the biomass (Froese and Pauly, 2019):

where a and b are allometric constants taken for each species from FishBase (Froese and Pauly, 2019), L is modal length value observed and W estimated weight. The units of length and weight in FishBase are centimeters and grams, respectively. The biomass of each species per transect is, then, the same as the estimated weight of the total number of individuals (Nb):

Biomass was used instead of abundance as it better represents the metabolism, secondary production, organic matter, and energy flows (Villéger et al., 2008; Cucherousset and Villéger, 2015; Rigolet et al., 2015; Cresson et al., 2019).

A functional matrix was built for 149 species registered in the PNSAV between 2006 and 2021, which included information on their life history, trophic ecology, and habitat use as established in previous studies of reef fish assemblages (Stuart-Smith et al., 2013; D’agata et al., 2014; Mouillot et al., 2014; Villéger et al., 2017; McLean et al., 2021). Six traits were selected: size, position in the water column, gregariousness, activity schedule, diet breadth, and trophic level ( Supplementary Material Table S1 ). All these traits have been widely used to represent biological and environmental interactions, including the influence that the species have on the ecosystem processes, and their responses to environmental change (Mouillot et al., 2013a; Mouillot et al., 2013b; McLean et al., 2019).

Information for the full species list was extracted primarily from FishBase (Froese and Pauly, 2019) using the ‘rfishbase’ package in R (Boettiger et al., 2012; R CoreTeam, 2022) and from Rincón-Díaz et al. (2018). Alternative sources were used whenever there was a gap in information, combining published research and reports, or information of the closest species (i.e., within the same genus) from a similar climatic region. These sources are referenced in Supplementary Material Table S1 . The continuous traits (size, diet breadth and trophic level) were weighted by their maximum value, thus obtaining scaled values between [0, 1] for each species. This process avoids the bias introduced by the different scales and magnitudes of the traits.

We followed the framework proposed by Maire et al. (2015) to build and select a representative functional space. The ‘qual_funct_space’ function was applied in R to evaluate the quality of the functional space, first producing a functional dendrogram using the UPGMA algorithm and multidimensional functional spaces (MDS) for two to 10 dimensions using a principal coordinates analysis (PCoA), with Cailliez’s correction for negative eigenvalues (Cailliez, 1983). The dendrogram and PCoA were computed from a Gower’s distance matrix because it allows the use of mixed data (Gower, 1971).

We considered a high-quality functional space when the distance between each pair of species was congruent with the initial functional distance. In this sense, the mean square deviation (mSD) was used to measure the deviation between the initial functional distance (Gower´s distance) and the scaled distance in the functional space (Cophenetic distance for dendrogram and Euclidian distance for MDS). The mSD is 0 when the functional space perfectly represents the initial distance and increases as pairs of species become less represented in the functional space (Maire et al., 2015).

Indices for functional richness (FRic), divergence (FDiv) and evenness (FEve) (Villéger et al., 2008) were calculated from the highest quality functional space obtained, using the {ape} and {geometry} libraries (Villéger et al., 2008; Mouillot et al., 2013a; Laliberté et al., 2015). FRic is defined as the amount of the functional space that is occupied by species, regardless of their biomass (Villéger et al., 2008). Low FRic values indicate that some of the potentially available resources are not being exploited and could result in an ecosystem productivity reduction (Villéger et al., 2010; Córdova-Tapia and Zambrano, 2015). FEve is the homogeneity in the distribution of the biomass of the species in the functional space (Mason et al., 2005; Villéger et al., 2008). Low FEve values implies that the species and their abundances are not homogeneously distributed in the functional space, which could reduce the productivity of the community as functional niches are underutilized and others overexploited. Low values may suggest an increased opportunity for potential invaders to establish themselves in the area (Mattingly et al., 2007; Córdova-Tapia and Zambrano, 2015). Finally, FDiv is a measure of functional similarity between the dominant species of a community (Mason et al., 2005: Villéger et al., 2008). The greater the FDiv in the community, the greater the functional difference between the dominant species of the community, which in theory favors the efficiency of the ecosystem by limiting competition and making better use of resources (Mouillot et al., 2007). The three indices were scaled by the maximum possible value, given all the species present in the set of assemblages so that they range from 0 to 1; hence, they are unitless and easily interpretable (Mouillot et al., 2013a).

The methodology to evaluate the percent contribution of each functional trait in the calculation of functional diversity indices proposed by Rincón-Díaz et al. (2018) was followed. First, we used linear regressions to identify functional traits that explained variation in Functional Diversity (FD). We then compared the coefficients of determination (R²) using the full pool of traits with the R² obtained when dropping each trait from the calculations, allowing us to identify the contribution of each trait to the FD, and to explain variations in FD metrics (Stuart-Smith et al., 2013). Low R² represented a significant gap in information and indicated traits that contributed more to variation in FD (Stuart-Smith et al., 2013).

Temporal hotspots of functional redundancy were identified following Rincón-Díaz et al. (2018) by plotting kernel density maps calculated from species locations in the functional trait space by years. Density maps were calculated in the PAST Program version 3.08 (Villéger et al., 2010) using a Gaussian kernel and locating nearby species within a radius of 0.02 within the trait functional space. Hotspots are formed by species with similar, but not equal functional traits.

We explored the main and joint effects of scale on the multivariate fish assemblage structure (based on abundances of individual species) by performing an unconstrained permutational multivariate analysis of variance (PERMANOVA, Anderson et al., 2008) using the ‘adonis’ function from the {vegan} R package. The factors we examined were subsystem (southern and northern), seascapes (inner shelf [IS], and outer shelf [OS]), reefs (seven in total), depth levels (shallow and deep), and years (including all interaction terms). Differences within the factors were considered significant at p ≤ 0.05.

For the 149 reef fish species described, with a mix of six categorical, ordinal, and continuous traits, the best functional space (FS) had four dimensions (mSD= 0.003) ( Figure 2 ). The quality of this functional space was slightly higher than the quality of spaces with five, six and seven dimensions, but was almost three times higher than the quality of the two-dimensional space (mSD= 0.0088) and 19.3 times higher than the quality of the best functional dendrogram (mSD= 0.058, built with the UPGMA clustering algorithm). We did not choose five, six or seven axes because: a) four axes have been shown to effectively capture the FD variation in assemblages, while requiring less computation time (Maire et al., 2015), and b) the resulting FS had the lowest root mean square deviations.

Species’ traits contributed differently to explain the variability in fish FD indices ( Table 1 ). Activity schedule, gregariousness and water column position were the most important traits to predict FD metrics (lowest R² values; Table 1 ). The activity schedule explained the greatest variation in FDiv (88%) and FEve (65%).

A total of 297 fish assemblages were observed from seven PNSAV reefs between May 2006 and April 2021. Significant differences were found in the FRic of the assemblages between subsystems (p= 0.001), years (p= 0.01), depth level (p= 0.02) and reef-depth interaction (p= 0.01), but not among the reefs (p= 0.23), or between seascapes (p= 0.11) ( Table 2 ). We found that the mean FRic was higher in the north (0.26, sd = 0.084) than in the south subsystem (0.22, sd = 0.089) with a p< 0.001 ( Figure 3 ). Although there were no significant differences between reefs alone (p = 0.23), the reefs that had the highest mean values of FRic were Verde, Anegada Adentro and Blanquilla with 0.28 (sd= 0.077), 0.27 (sd= 0.087), and 0.25 (sd= 0.087), respectively, while the lowest mean value was found in Blanca reef with 0.20 (sd= 0.087; Figure 4 ). In general, the mean values of FRic did not exceed 0.28, which implies that, on average, the observed fish did not occupy more than 28% of all the theoretical FS for this system. Temporally, the FRic of the community varied significantly (p = 0.01). The highest annual FRic mean values were found in 2007 (0.31, sd= 0.062), followed by 2013 (0.27, sd= 0.089), after of which a rapid drop occurs reaching the minimum value observed in 2014 (0.18, sd= 0.068). In more recent years, values of FRic recovered gradually until reaching (0.25, sd= 0.057) in 2017 and remain without major changes until 2021.

The FDiv presented annual mean values between 0.83 (sd= 0.085) and 0.90 (sd= 0.044; Figure 4 ), without statistical differences between years (p = 0.44; Table 2 ). In contrast, statistical differences were found between subsystems (p= 0.03), depth level (p= 0.01) and the seascape-depth interaction (p= 0.01), opposite to the reef and seascape p= 0.048 and p= 0.67, respectively; Table 2 ). We found a higher mean FDiv in the north (0.90, sd= 0.052) than in south subsystem (0.88, sd= 0.072; Figure 3 ).

The annual mean FEve values varied between 0.34 (sd= 0.128) in 2013 and 0.44 (sd= 0.06) in 2017, with significant differences between years (p= 0.02), reefs (p= 0.02) and reef-depth level interaction (p= 0.02), but no differences were found between subsystems (p= 0.47), seascapes (p= 0.09) and depth-level alone (p= 0.09; Table 2 ). The reefs with the highest mean FEve values were Sacrificios (0.42, sd= 0.11), Verde (0.41, sd= 0.10) and Enmedio (0.40, sd= 0.10), while the lowest was from Blanquilla reef (0.35, sd= 0.11; Figure 4 ). There was no consistent trend throughout the period in FEve ( Figure 4 ); rather, it had an erratic behavior with significant increases and decreases.

Fish assemblages of the PNSAV exhibited variation of FD through time ( Figure 5 ). Hotspots in the functional trait space were reduced after 2014 and never returned to their initial states ( Figure 5 ). The configuration of functional redundancy shows two conspicuous hotspots of species’ density with similar traits combinations. The hotspot located in upper left quadrant (ULQ1) include species like Synodus intermedius, Bodianus rufus, Malacoctenus triangulatus, Microspathodon chrysurus and Elacatinus jarocho, which are characterized by being solitary or in small groups of small to medium sized individuals, intermediate trophic levels, and a strong association with the bottom, many inhabit small cavities, corals or live buried in the sand. The second hotspot is denser and is situated in the lower left quadrant (LLQ2) including species like Mycteroperca bonai, M. phenax, M. venenose, Epinephelus adscensionis, Calamus calamus and Serranus subligarius. These species are solitary with high trophic levels, are demersal and predator species of medium to large sizes with crepuscular and nocturnal behaviors, and a broader diet than ULQ1. In 2014, the hotspot ULQ1 density increased while LLQ2 density decreased. After 2014, the density of species with similar trait combination decreases drastically in all trait space suggesting a reduction in fish functional redundancy.