Functional diversity in reef fish assemblages in the Parque Nacional Sistema Arrecifal Veracruzano, Mexico: Temporal and spatial changes

The Parque Nacional Sistema Arrecifal Veracruzano (hereafter called PNSAV) is the largest coral reef extension in the central region of the Gulf of Mexico. These reefs are unique since they have developed near a coastal environment that is directly influenced by the discharges of Veracruz city, the rivers located on the continental shelf, and the port of Veracruz. This study evaluates the functional diversity, in terms of richness, evenness, and divergence, of the PNSAV fish community. We were interested in quantifying any similarities or differences in functional diversity metrics when one examines reef fish assemblages on a single reef or joint reef subsystems; thus, is there a difference based on scale? A total of 297 fish assemblages were observed in seven PNSAV reefs between May 2006 and April 2021. Significant differences were found in the Functional Richness of the assemblages between subsystems, years, and reef-depth interaction, but none were found among the reefs, or between seascapes. The Functional Diversity presented annual mean values between 0.83 (sd= 0.085) and 0.90 (sd= 0.068) and did not show statistical differences between years, seascape, or reefs. In contrast, statistical differences were found between subsystems, and depth level and the seascape-depth interaction. The annual mean Functional Evenness values ranged between 0.34 (sd= 0.128) and 0.44 (sd= 0.060), and significant differences were detected between years, reef, and reefs-depth level interaction, but no difference were found between subsystems. Reef-fish diversity was greater within the north coral reef subsystem than the southern of the PNSAV. There were no overall tendencies for increased functional diversity throughout the time during this study.


Introduction
Functional diversity describes the ecosystem functioning based on the species' performance in the environment, their evolutionary responses to environmental variability, and the interactions among coexisting species (Nagelkerken and Munday, 2016;Villeǵer et al., 2017;Goldenberg et al., 2018). The species' functional roles are measured by unique combinations of biological or ecological traits (i.e., morphological, physiological, and behavioral) that likely enhance their survival (Dıáz and Cabido, 2001;Petchey and Gaston, 2006;Violle et al., 2007;Cadotte et al., 2011); thus, by measuring them we gain insights into the complexity and the structuring processes of communities.
The range and types of functional traits used enable resource managers to address several ecological and management issues (Villeǵer et al., 2017), yielding successful results across various habitats and taxa (Brandl and Bellwood, 2014;Mouillot et al., 2014;Ramıŕez-Ortiz et al., 2017;Rincoń-Dıáz et al., 2018;Palacios-Salgado et al., 2019;Teittinen and Virta, 2021). Studies of the functional structure of fish communities can detect small changes in the community structure, which makes them a good tool for conservation (Freitas and Mantovani, 2018) and can complement traditional descriptors in biogeographical studies to assess ecosystem health, resilience, and responses to environmental disturbance more robustly (Tilman et al., 1997;Bellwood et al., 2004;Cadotte et al., 2011;Mouillot et al., 2011).
Reef fish functional diversity studies have highlighted the importance of fish in the resilience of ecosystems (Micheli and Halpern, 2005;Mouillot et al., 2014), the regulation of food webs and nutrient recycling (Holmlund and Hammer, 1999;Duffy et al., 2016) from tropical to temperate ecosystems (Stuart-Smith et al., 2013;Quimbayo et al., 2019), while also showing that significant changes can occur even within the same ecoregion (Duffy et al., 2016). Thus, the spatial scale one examines an ecosystem using this approach may influence or bias our understanding of how that ecosystem functions.
The Sistema Arrecifal Veracruzano (SAV) is the largest reef extension in the central region of the Gulf of Mexico (GoM) (Spalding et al., 2001). Although in 1992 it was declared a National Park (Parque Nacional Sistema Arrecifal Veracruzano, PNSAV) (Diario Oficial de la Federación [DOF], 1992), it is under substantial and continuous stress due to potential negative anthropogenic factors like harbor activity, hydrocarbon spills, artisanal fishing, tourist activity, and discharges of urban pollutants (Peŕez-España and Vargas-Hernańdez, 2008). Because of that, research information on its structural and biological communities' changes has been a priority.
Hence, this study evaluates the fish functional diversity of the PNSAV reef community from 2006 to 2021, considering its subsystems (southern and northern) and seascapes (inner shelf [IS], and outer shelf [OS]), each reef (seven different reefs), their depth level (shallow and deep), and the year in which they were sampled. We specifically hypothesized that the functional diversity of reef fishes would be greater within the southern coral reefs than the northern, due to their relative remoteness to the pollution of the city; hence making the subsystem the most consistent and important scale, with particular patterns occurring within each one. Consequently, we propose that a more accurate approach would be to define geographical subsystems at specific depths in which the traits are divergent from one another, and that protection and management policies should concentrate on areas with low trait diversity/redundancy.

Study area
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 Federacioń [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ıó. The south subsystem has 19 reefs, located in front of the community of Antoń 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 (Peŕez-España and Vargas-Hernańdez, 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 2 , and a body height of between 20 and 45 cm (Allende-Arandıá 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-Peŕez 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-Chavez and Lango-Reynoso, 2021).

Sampling
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 (Peŕez-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 (Villeǵer et al., 2008;Cucherousset and Villeǵer, 2015;Rigolet et al., 2015;Cresson et al., 2019).

Data information
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;Villeǵer 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 Rincoń- Dıáz 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.

Functional space
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ś 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).

Functional diversity indices
Indices for functional richness (FRic), divergence (FDiv) and evenness (FEve) (Villeǵer et al., 2008) were calculated from the highest quality functional space obtained, using the {ape} and {geometry} libraries (Villeǵer 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 (Villeǵer 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 (Villeǵer et al., 2010;Coŕdova-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;Villeǵer 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: Villeǵer 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).

Sensitivity of functional diversity metrics to chosen traits
The methodology to evaluate the percent contribution of each functional trait in the calculation of functional diversity indices proposed by Rincoń-Dıáz 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).

Kernel density maps
Temporal hotspots of functional redundancy were identified following Rincoń- Dıáz 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 (Villeǵer 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.

Main and interaction effects
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.

Functional space and sensitivity of functional diversity metrics to chosen traits
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%).

Functional diversity indices
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 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).

Kernel density maps
Fish assemblages of the PNSAV exhibited variation of FD through time ( Figure 5). Hotspots in the functional trait space were reduced  The quality of functional spaces used to compute functional diversity metrics for the PNSAV fish assemblages sampled between 2006 and 2021. The top-left panel shows the mean squared deviation for two to ten functional spaces and dendrogram. The remaining ten panels show the correlation between the initial Gower's distances and the standardized distances in each of the ten functional spaces and dendrogram. The mean Squared-Deviation (mSD) equals the average deviation between Gower's distance and Euclidean distance for multidimensional spaces and Cophenetic distance for dendrogram, mSD can be viewed like an average error or quality space metric. Each point represents a species.
after 2014 and never returned to their initial states ( Figure 5).

Discussion
Despite the PNSAV being a protected area for thirty years, and the relevance of functional diversity metrics to describe ecosystem functioning, responses to environmental variability and the interactions among coexisting species, this is the first study that evaluates the functional diversity of the PNSAV's reef fish community. We found that reef fish´s functional structure varied with the spatial scale of analysis, based on the effect of the geographical position of the subsystem (north or south) in the FRic and FDiv. Thus, our results indicate that, rather than a site priority approach to protection and management, it would be better to define geographical subsystems at specific depths in which the traits are divergent from one another. Protection and management should concentrate on areas where: a) low trait diversity or redundancy is detected (Trindade-Santos et al., 2022), and b) populations are highly aggregated in time or space (Edgar et al., 2008), while considering human activities like tourism and sport fishing (Comisioń Nacional de Áreas Naturales Protegidas [CONANP], 2017;Tuohy et al., 2022).
While it would be expected that the less diverse sites were located closer to human activities, we found the opposite: the south subsystem (further away from the activities) exhibited lower values of FRic than the north subsystem (more urbanized area). Usually, niches tend to be exploited more efficiently on sites with suboptimal conditions, leading to a lower redundancy due to the usually reduced competition among fewer species. Thus, species are more likely to occupy more functions at less altered sites, leading to higher FRic values (Stuart-Smith et al., 2013;McLean et al., 2019). The opposite pattern that we found could, then, be explained by factors other than human activities alone.
One of such factors could be the coral coverage (structural complexity; Richardson et al., 2017). Higher FRic values were found in sites where Peŕez-España et al. (2015) reported a high coral coverage found at 7 m or deeper. The fish assemblages in these  Boxplot with spatial variation in functional diversity (FD) metrics by subsystem, depth level, and seascape region; values of FD metrics are express as a mean rather than median. Functional diversity metrics were calculated with four coordinate axes that explained 73.87% of the cumulative total variation of fish species in the functional trait space. Functional richness (FRic), divergence (FDiv) and evenness (FEve), Inner Shelf (IS), and Outer Shelf (OS).

FIGURE 4
Boxplot with spatio-temporal variation in functional diversity (FD) metrics by reef and years. Values of FD metrics are express as mean rather than median. Functional diversity metrics were calculated with four coordinate axes that explained 73.87% of the cumulative total variation of fish species in the functional trait space. Functional richness (FRic), divergence (FDiv), and evenness (FEve).
habitats appear to have a broader range of traits and greater abundances, which would likely create a higher complementarity of niches and a more efficient resource exploitation (Richardson et al., 2017;Yeager et al., 2017). The depth level is another variable that seems to have an important effect, either direct and/or as a causation of other conditions that were identified to be important in the community functional structure (coral coverage, less impacted sites, reduce discharge of continental waters). Other studies focusing on deep-sea functional communities have reached a similar conclusion, both directly (Mindel et al., 2016) and indirectly (Carrington et al., 2021;Myers et al., 2021). Bringing these ideas together, habitat filtering tends to be the guiding force in fish assemblages, having stronger effects based on the characteristics of local sites (Richardson et al., 2017;Valleé et al., 2019). Deep and exposed reefs tend to be dominated by solitary largebodied species (Quimbayo et al., 2017), as in our study, while being more functionally heterogeneous than shallow reefs (Medeiros et al., 2021). In contrast, low structural complexity reefs with low coral coverage tend to favor the presence of medium and large-bodied species, usually in high densities (Quimbayo et al., 2017). This fish sizecoral coverage relationship could be supported by the positioning and discharge of the Jamapa, Papalopan and La Antigua rivers, where the difference of salinity and sedimentation between the shallow and deeper areas enhance the food availability for demersal species with smaller sizes and higher population densities. These discharge patterns may not influence deeper reefs (Peŕez-España et al., 2015;Castañeda-Chavez and Lango-Reynoso, 2021) as the coral coverage would be less affected (Fabricius, 2005).
Although significant changes were detected for the FRic and FEve throughout the 14-year time course of this study, there was not a clear increase or decrease in neither indices, but rather two cycles of increasing and decreasing density spanning a four-year window. Similar patterns were found by Peŕez-España et al. (2015) for species richness and species density. These similarities are likely due to the strong correlation between functional richness and number of species in our data (rho = 0.855, p< 0.001; Supplementary Material Figure S2) (Villeǵer et al., 2010). The changes in FRic and FEve, but not in FDiv across the 14 year time frame implies that dominant species have extreme functional trait values (remained functionally different). In a multivariate context, functional divergence relate to abundance is distributed within the volume of functional trait space (Villeǵer et al., 2008). These species are defined as well as specie vertices (Supplementary Materials Figure  S3 and Table S2). This could be related to how local conditions favor certain traits (Jeppesen et al., 2010), creating a structured reef assemblage where only one dominant species group thrives and which 'guides' the reef food web (McLean et al., 2019;Robinson et al., 2019).
However, a change in functional space did occur, as most traits shifted from a superficial redundant community (density of superficial species) to a benthic one and later becoming more evenly balanced in the range of traits. The tipping point coincided with the development of the marine port of Veracruz in 2014 (Jimeńez-Badillo et al., 2014), where the Sacrificios reefs were nearby and the most affected by the development. Some studies had detected alterations in community structure at sites close to construction areas (Gehrke et al., 2002;Hang Limbu et al., 2021), with the impacts being especially problematic on the more mobile fishes. In conjunction with the development itself, the increase in average rainfall from 2013-2016 (Comisión Nacional del Agua Temporal variation of functional redundancy in fish assemblages of the Parque Nacional Sistema Arrecifal Veracruzano (PNSAV). Maps show boundaries of the functional trait space thru time with species location in the trait space as dots, and hotspots of species with similar trophic function as reddish regions. Functional trait spaces were plotted with the first two axes of the principal coordinate analysis for graphical representation where coordinate 1 and coordinate 2 explained 23.98% and 22.57% of the cumulative total variation respectively. The scale bar gives an estimate of the number of species per area based on species distances within a radius of 0.097 within the trait functional space. Gradients of functional traits are shown for the first and second coordinate axes only for this graphical representation.
[CONAGUA], 2019) likely increased the runoff, thus increasing sedimentation and altered salinity conditions which may be additional confounding issues (Edinger et al., 2000;Risk, 2014). These alterations may enhance the dominant functional species traits per site as documented in other species and habitats were changes in salinity shape the functional community structure (Costa e Silva et al., 2019;Payan-Alcacio et al., 2020;Borba et al., 2021). However, there is not enough in situ information available to back up this claim. Another caveat is that the effect of the environmental variability on the functional diversity of coral reef fishes is scaledependent and may affect each metric differently (Yeager et al., 2017).
The high values of FDiv documented in the PNSAV indicate specialized trait communities (Payań-Alcacio et al., 2021). Without a focus on which functional process dominates, it is difficult to delineate which trait or set of traits (e.g., locomotion, feeding, habitat use) would most likely be related to depth and food availability. Depth was shown to be important in our results, and feeding habits are a significant variable for community structure differentiation (Ferreira et al., 2001;Deus and Petrere-Junior, 2003;Carrington et al., 2021). It has been widely suggested that the primary mechanisms regulating reef fish community structure may vary as a function of reef type, reef context, and geographic region (e.g., Beukers and Jones, 1998;Friedlander et al., 2003;Darling et al., 2017;Wetmore et al., 2020). Thus, the high FDiv we documented could be an accumulation of other causes like anthropogenic activities, disruption of fauna due to vessel transportation (over 1700 vessels are registered in the PNSAV), different coral species habitats, and coral habitat fragmentation and bleaching (Richardson et al., 2017;Ortiz-Lozano et al., 2018).
From a management and environmental protection perspective, it is necessary to consider and understand key ecological indicators and processes like a functional ecology approach, as well as biological composition (Vora, 1997) and socio-economic components (Tuohy et al., 2022). In general, the current protected fishing areas of PNSAV (Reyna-Gonzaĺez et al., 2014) do not coincide with those with greater functional diversity. The great mobility of fishes implies that any disturbance can be dispersed between reefs that are inside and outside the shelf; differences between depths and FRic paired with low redundancy point towards a fragile ecosystem with high vulnerability. This study highlights the importance of incorporating functional diversity studies into the monitoring programs of natural protected areas.

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.

Ethics statement
Ethical review and approval was not required for the animal study because research was only observational.

Author contributions
AR-R, idea developer, processing, and data analysis. HP-E, idea developer, project responsible, funding administrator, data collection, processing, and data analysis. JP-A, GD-A, and MP, data analysis. AE-G and AO, processing and data analysis. VC-E, idea developer and data analysis. All authors contributed to the article and approved the submitted version.

SUPPLEMENTARY TABLE 2
List of species vertices shaping the convex hull of 149 reef fish species present in Parque Nacional Sistema Arrecifal Veracruzano between 2006 and 2021.

SUPPLEMENTARY FIGURE 1
Cartoon scheme shows hierarchical sampling design of Parque Nacional Sistema Arrecifal Veracruzano between 2006 and 2021. Both subsystem (north and south) was sampled at two different seascapes (inner shelf and outer shelf). The observations in reefs at different seascapes were made at two depth levels (shallows and deep).

SUPPLEMENTARY FIGURE 2
Spearman rank correlation analysis (rho) between metrics of fish functional diversity, and the number of species (Nb_sp). Position of 149 reef fish species of Parque Nacional Sistema Arrecifal Veracruzano community are represented along four functional axes. Vertices of the convex hull are in purple whereas non-vertices are in green circles.