Marine reserves lag behind wilderness in the conservation of key functional roles

Although marine reserves represent one of the most effective management responses to human impacts, their capacity to sustain the same diversity of species, functional roles and biomass of reef fishes as wilderness areas remains questionable, in particular in regions with deep and long-lasting human footprints. Here we show that fish functional diversity and biomass of top predators are significantly higher on coral reefs located at more than 20 h travel time from the main market compared with even the oldest (38 years old), largest (17,500 ha) and most restrictive (no entry) marine reserve in New Caledonia (South-Western Pacific). We further demonstrate that wilderness areas support unique ecological values with no equivalency as one gets closer to humans, even in large and well-managed marine reserves. Wilderness areas may therefore serve as benchmarks for management effectiveness and act as the last refuges for the most vulnerable functional roles.


Supplementary Figure 4. Partial dependence plots of biomass and biodiversity indices for different categories of areas.
Fitted variations were predicted using 8 environmental explanatory variables and the "management" as the human variables as predictors in the BRT models. The left y-axis is the percentage of variation from the maximum value for each community aspect. The percentage of the maximum value is independent of the range and the unit of each index and thus comparable between indices. The fitted levels of fish community aspects for different managements: i) exploited areas (< 3h) (red), ii) the no-take small MPAs (black), iii) the no-entry large MPA (black), iv) the traditionally managed areas < 5 hours travel time (grey), v) the traditionally managed areas at 19 hours travel time (grey) and vi) the wilderness areas (> 20h) (blue). Figure 5. Contributions of explanatory variables for different fish community aspects using the simplified "Management" model. The contributions of each explanatory variable (%) from "Management" simplified BRT models are given for a) the Total Biomass (g.m -²), b) Species Density (Number of species per transect), c) Herbivores biomass (g.m -²), d) Functional Richness (FRic), e) Apex biomass (g.m -²) and f) Biomass-weighted functional diversity (Rao entropy) converted to equivalent number of species. The variable "Management" is highlighted in light grey.

Supplementary Figure 6. Functional space with PCoAs axes characterization.
Ecological characterization of the functional space with five most important functional traits related to PCoA axes 1 and 2: "Size" of fish, "Home-range" or the mobility of fish, 3) "Schooling" or the gregariousness of fish species, 4) "Level" or the vertical position in the water column and 5) the "Diet" with 7 categories (HD: herbivorousdetritivorous, HM:macroalgal herbivorous, IS: invertivorous targeting sessile invertebrates, IM: invertivorous targeting mobile invertebrate, PK: planktivorous, FC: piscivorous, and OM: omnivorous). The range of coordinates along PCoA axes for the different modalities of ordered qualitative functional traits were used to define the trait modalities characterizing PCoA axes. For the "Size", "Home range", "Schooling" and "Level" traits, the larger part of arrows indicate the highest values of the modalities. For the "Diet" trait, the range of values is indicated for each modality In black are raw data for "no-managed areas", in dark gray for "traditional areas" and in light grey for marine protected areas.

Supplementary Table 1. Parameters and performance of simplified BRT models.
Parameters and performance of simplified BRT models for Biomass of commercial fish (B), Biomass of apex predators (B APEX) and biomass of herbivores (B HERB.), Species density (S), Rao Equivalent of number of species (Rao Eq.) and Functional Richness (FRic). lr is the learning rate, tc is the tree complexity, N.trees is the number of optimal trees, R²Tr is the correlation coefficient based on training dataset, R²CV is the correlation coefficient from the k-fold cross-validation procedure with n folds (nf) equal to 20 and bag fraction (bf), SE R² CV is the standard error. D 2 is the cross-validated proportion of the total deviance explained.

Supplementary Table 4. Wilcoxon and non-parametric effect size pairwise comparisons between management levels.
Wilcoxon (W test and p-value) and non-parametric effect size (r) pairwise comparisons between the 6 levels of management for: i) exploited areas ("< 3h"), ii) the no-take small MPAs, iii) the no-entry large MPA, iv) the small traditionally managed areas under 5 hours of travel time, v) the traditionally managed areas at 19 hours of travel time and vi) the wilderness areas ("> 20h)".The results of the post-hoc Kruskal-Wallis test are indicated for each community aspect with letters. The results of the pairwise Wilcoxon tests (W value and its p-value) are indicated. The non-parametric effect size is the ratio of the z score from the Wilcoxon test divided by the square-root of the number n of transect (r = z/√ ). (ns), p>0.05; * p ≤ 0.05; ** p ≤ 0.01; *** p ≤ 0.001; *** p ≤ 0.0001. At a very small scale we may have some autocorrelation and thus pseudoreplication which inflates p-values.  Supplementary Table 5. Percentage of regional functional volume filled by the six different categories of areas.
The percentage of the regional functional volume for 352 commercial reef fish species have been computed as the ratio between the functional volume shared by 50% of the transects divided by the regional functional volume. n is the number of transect per location.