# Conserving the beauty of the world's reef fish assemblages

Matthew Mclean, David Mouillot, Juliette Langlois, Suchinta Arif, Sonia Bejarano, Nicolas Casajus, Graham Edgar, Ulysse Flandrin, François Guilhaumon, Aaron Judah, et al.

► To cite this version:

Matthew Mclean, David Mouillot, Juliette Langlois, Suchinta Arif, Sonia Bejarano, et al. Conserving the beauty of the world's reef fish assemblages. Proceedings of the National Academy of Sciences of the United States of America, 2025, 122 (25), pp.e2415931122. 10.1073/pnas.2415931122. hal-05140021

HAL Id: hal-05140021

https://hal.science/hal-05140021v1

Submitted on 7 Nov 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
1 Conserving the beauty of the world's reef fish assemblages

2

3 Matthew McLean¹,*, David Mouillot², Juliette Langlois², Suchinta Arif³, Sonia

4 Bejarano⁴, Nicolas Casajus⁵, Graham J. Edgar⁶, Ulysse Flandrin², François

5 Guilhaumon⁷, Aaron B. Judah⁸, Nicolas Loiseau², M. Aaron MacNeil⁹, Eva Maire²,¹⁰,

6 Rick D. Stuart Smith⁶, and Nicolas Mouquet²,⁵,*

7

8 Affiliations

9

¹Department of Biology and Marine Biology, Center for Marine Science, University of North Carolina Wilmington, Wilmington, NC 28403. USA.

10

11 ²MARBEC, Univ Montpellier, CNRS, Ifremer, IRD, Montpellier, France.

12 ³Department of Biology, Dalhousie University, Halifax, NS B3H 4R2

13 ⁴Leibniz Centre for Tropical Marine Research (ZMT), Fahrenheitstraße 6, 28359 Bremen, Germany

14

15 ⁵FRB - CESAB, 34000 Montpellier, France.

16 ⁶Institute for Marine and Antarctic Studies, University of Tasmania, Hobart, Tasmania 7001, Australia.

17

18 ⁷IRD, la Réunion, France.

19 ⁸Department of Oceanography, University of Hawai'i at Mānoa, Honolulu, HI, United States

20 ⁹Ocean Frontier Institute, Department of Biology, Dalhousie University, Halifax, NS B3H 4R2

21 ¹⁰Lancaster Environment Centre, Lancaster University, Lancaster, LA1 4YQ, UK

22

23 * These authors contributed equally to this study

24

25 **Corresponding authors:** Matthew McLean: [mcleamj@gmail.com](mailto:mcleamj@gmail.com) Nicolas Mouquet:

26 [nicolas.mouquet@cnrs.fr](mailto:nicolas.mouquet@cnrs.fr)
## 27 Abstract

On shallow rocky and coral reefs, cultural and recreational values, like aesthetics, are critical aspects of Nature's Contributions to People (NCP) that support human well-being and provide billions of dollars in tourism revenue. Quantifying the aesthetic value of reef ecosystems and uncovering the conditions that enhance it could support NCP-based management. Here, we combine a global dataset of reef fish surveys, species-level aesthetic values, and causal modeling to assess the global status and drivers of reef fish assemblage aesthetic value. We find that aesthetic value is inherently linked to species richness, displaying a latitudinal gradient with peaks in the tropics, but varies strongly with the presence of exceptionally beautiful or less-beautiful species. Sea surface temperature, primary productivity, human gravity, and protection status are the strongest drivers of assemblage-level aesthetic value. Protection against human impacts consistently enhances aesthetic value by boosting taxonomic and phylogenetic diversity, and this effect is greatest in species-rich, tropical ecoregions. Economic development has little influence, indicating that low-income countries are not constrained from maintaining beautiful fish assemblages. Our results therefore suggest that marine protected areas (MPAs) can support multiple NCPs simultaneously, particularly in developing tropical countries. While we highlight the effectiveness of MPAs, given the low level of marine protection globally and the sensitivity of aesthetic value to environmental conditions, the beauty of the world's reefs appears severely threatened. Aesthetic value should be immediately integrated into reef conservation and management plans.

49

## 50 Significance Statement

Ecosystems provide critical services to human societies, the most well-known including food provisioning and coastal protection. However, ecosystems also provide cultural services that are essential to human well-being and economic livelihood. The aesthetic value of ecosystems creates an intimate connection between people and nature and supports tourism. We evaluated the distribution and drivers of the aesthetic value of fish assemblages on shallow reefs worldwide. We found that protected areas enhanced aesthetic value by harboring greater species richness and exceptionally beautiful species, which was strongest in tropical regions. We also found that aesthetic value was not related to countries' level of wealth. Protected areas can therefore
60 support fisheries, tourism, and human well-being simultaneously, and tropical,
61 developing countries can benefit most.

62

63 **Introduction**

64 Evaluating how and where ecosystems support human livelihood is critical for
65 achieving sustainable development goals. Nature's Contributions to People (NCP) is
66 a growing framework for assessing links between ecosystems and human society.
67 While NCP encompasses the whole spectrum of interactions between nature and
68 people (1, 2), previous studies have focused mainly on material goods, ecosystem
69 productivity, or invasive species resistance. Few studies have quantified non-material
70 dimensions like cultural importance, wellbeing, or the perceived aesthetics of
71 biodiversity (3). However, such non-material NCP can have a substantial role in
72 conservation biology (4, 5). Society's willingness to engage in conservation action is
73 not solely driven by what nature can offer in terms of material benefits, but also by how
74 people experience the natural environment (6). Conservation values often stem from
75 the emotional connection that humans have with nature (7), as there is a clear
76 association between affective factors (emotions that impact human decision-making)
77 and people's willingness to safeguard species and ecosystems (8, 9).

78 Non-material value refers to the cultural relationship people have with nature
79 and the impact nature has on people's mental state, creativity, and overall sense of
80 well-being (2). This underscores the importance of preserving and promoting natural
81 environments for the benefit of human health and happiness. Among non-material
82 NCP, aesthetic value is one of the most direct links people have with nature (5).
83 Quantitative measures of perceived species aesthetics (hereafter 'aesthetics') are now
84 becoming available, making comparisons possible with other NCP metrics (10, 11).
85 Yet estimates of aesthetics at the biological assemblage level (a relevant scale for
86 conservation) are still scarce (5) and do not account for the spatial distribution of
87 organisms or environmental characteristics (12, 13). To include biodiversity aesthetics
88 in the NCP agenda and to increase our capacity to understand its drivers and protect
89 it, we need broad-scale assessments across diverse social and environmental
90 conditions.

91 Shallow rocky and coral reef fish assemblages offer a unique opportunity to
92 quantify aesthetics. Among ecosystems emblematic for their aesthetic value, reefs
host exceptional biodiversity, including over 6,000 fish species (14), many that have
evolved vivid coloration and body shapes used in camouflage and communication (15,
16). Reef fishes provide at least two key NCPs: i) their aesthetic value is highly
important to human cultural heritage and is an iconic example of how biodiversity can
be experienced as beautiful, and ii) they have strong economic importance as a vital
food source (17) and by supporting aquarium and tourism industries (18). Shallow
reefs are indeed important destinations for scuba diving and snorkeling, attracting
millions of tourists worldwide and generating an estimated $36 billion per year in
revenue (18). Although quantitative aesthetic values are now available for more than
2,000 rocky and coral reef fish species, how these values scale up to fish assemblages
around the world has not yet been assessed. Moreover, the environmental and human
conditions that enhance (or degrade) fish assemblage aesthetic value have not been
identified.

Here, using global, standardized data from the Reef Life Survey (RLS) program (19), we estimated and mapped the aesthetic value of fish assemblages from 3,526 sites worldwide and identified aesthetic hotspots and coldspots. We then assessed the potential environmental and anthropogenic drivers of aesthetic value using causal inference and Bayesian regression models. Finally, we disentangled the pathways by which marine protected areas (MPAs) may affect aesthetic value and analyzed geographic variability in MPA effectiveness. Our study provides the first quantitative evaluation of the potential to conserve the aesthetics of these iconic and endangered ecosystems, providing direction for future management.

115

## 116 Results

117 *From fish species to assemblage aesthetic values*

Our assemblage-level estimates of aesthetic value are indirect and based on previous
work by Tribot et al. (13) and Langlois et al. (10) These authors directly quantified reef
fish assemblages' aesthetic values through online surveys with >2,000 respondents
who ranked fish assemblages by their perceived beauty. The authors similarly
quantified individual species' aesthetic values based on >10,000 survey respondents
(5, 13). These studies found no influence of socio-cultural backgrounds of survey
respondents, allowing a robust estimation of species and assemblage aesthetic
values. Using linear models, Tribot et al. (13) showed that assemblages' aesthetic
values are primarily determined by the combination of species richness and species
aesthetic values – assemblages with more species are considered more beautiful,
however, certain species can increase or decrease the overall value (13). Using model
coefficients from Tribot et al. (13), we calculated: 1) the aesthetic value of each fish
assemblage - the total aesthetic value based on both species richness and species
aesthetic values (weighted by species abundances), and 2) the deviation from the
value expected based purely on species richness - how much higher or lower an
assemblage's aesthetic value is than expected given how many species there are (see
Methods). We weighted species aesthetic values by abundances because highly
beautiful species will not contribute importantly to human perception of aesthetic value
if they are not commonly observed and vice versa. Worldwide, we found that
assemblage-level aesthetic values increase with species richness as expected (Fig.
1B) but that assemblages with the same species richness vary over 5-fold due to
compositional differences (Fig. 1C). Even though species-rich assemblages are more
likely to contain species with high aesthetic values, assemblage-level aesthetic values
are unevenly distributed along the species richness gradient with most values being
below predictions at low species richness (Fig. 1B).

143

144 **Global distribution of fish assemblage aesthetic values**

Assemblage aesthetic values closely mirrored global patterns in species richness,
increasing toward the equator and decreasing toward the poles (Fig. 2A). The
Maldives, Seychelles, Vanuatu, Ningaloo, and Western Sumatra ecoregions are
aesthetic hotspots with the highest average values (Fig. 2A). Conversely the
ecoregions with the lowest average values were the North and East Barents Sea, East
Greenland Shelf, Channels and Fjords of Southern Chile, Oyashio Current, and North
Patagonian Gulfs (Fig. 2A). The deviations from the values expected given species
richness (i.e. residual aesthetic values not explained by species richness), followed a
similar pattern (r = 0.76 between assemblage aesthetic values and residual values
across all ecoregions; Fig. 2C). However, several ecoregions deviated significantly in
both positive and negative directions. For instance, Vanuatu and the Banda Sea both
contained ~70 species per site on average, yet Vanuatu had a much higher average
aesthetic value than expected while the Banda Sea had a much lower value than
expected. Interestingly, there was also high variation within many tropical ecoregions,
meaning that neighboring sites could have contrasting aesthetic values, likely driven
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
by local variation in species composition. The Great Barrier Reef (GBR) had many sites with aesthetic values both substantially lower and higher than expected based on species richness (Fig. 2C). Looking closer revealed that GBR sites with high aesthetic values were more characterized by Acanthuridae, Cirrithidae, and Balistidae, whereas low-value sites where more characterized by Caesionidae, Carangidae, and Lutjanidae (SI Appendix, Fig. S1).

## Drivers of fish assemblage aesthetic values

Following the structural causal model (SCM (20)) framework, we generated a Directed Acyclic Graph (DAG) that was used to guide statistical model building (see Methods). We built individual Bayesian regression models for each anthropogenic and environmental predictor variable, controlling for necessary covariates through the SCM selection framework. From each model, we extracted the posterior distribution (effect size) for the predictor variable of interest and then aggregated all posterior distributions into a single visual plot. Several predictor variables in this study were assessed using ordination axes (see Methods). Benthic composition was assessed by applying a Principal Components Analysis to benthic substrate proportions recorded from the survey transects and using the first two PC axes as predictor variables (which primarily characterized gradients in live coral vs macroalgal cover; SI Appendix, Fig. S2). In total, we assessed 10 potential predictor variables (SI Appendix, Table S1), resulting in 10 regression models. All models had excellent performance, with $R^2$ values of ~0.88 (for all) when considering both random and fixed effects (of which fixed effects accounted for up to 0.58).

Assemblage aesthetic values were primarily driven by environmental gradients in sea surface temperature and net primary productivity (Fig. 2B), with temperature having a strong positive effect (effect size and 90% credible interval = 0.46 [0.42, 0.50]), and productivity having a strong negative effect (-0.15 [-0.17, -0.14]). Assemblage aesthetic values were also influenced by benthic composition (PC1 axis) and depth, with higher aesthetic values being associated with a higher proportion of live coral cover (PC1: 0.046 [0.033, 0.059]; SI Appendix, Fig. S2), and shallower depths (-0.030 [-0.038, -0.023]). Beyond environmental factors, aesthetic values were also shaped by human pressure and conservation status. Human gravity had a negative impact on aesthetic values (-0.047 [-0.061, -0.034]), while restricted-take (0.049 [0.038, 0.060]) and no-take (0.075 [0.063, 0.086]) MPAs had positive effects.
At the country level, Human Development Index (HDI) did not have an important influence, but countries with higher fisheries dependency had higher aesthetic values, although this effect was variable (Fig. 2B). It is important to note, however, that alternative model structures (that included a spatial Gaussian Process) provided slightly weaker effects for some variables (see SI Appendix, Supporting Methods and Fig. S3-4).

We also examined the drivers of the deviation from the values expected based on species richness. Across nearly all predictor variables, we found similar results, but with lower effect sizes (Fig. 2B,D). This means that environmental and anthropogenic factors largely influenced aesthetic value through species richness, but not exclusively. Sea surface temperature, net primary productivity, benthic composition, and no-take MPAs remained the most influential drivers (Fig. 2D).

206

## Effect of no-take MPAs on fish assemblage aesthetic value

After identifying the overall positive effect of no-take MPAs on assemblage aesthetic value, we examined whether this effect varied among ecoregions (see Methods). We found that the effect of no-take MPAs varied predictably across ecoregions, with a clear unimodal latitudinal pattern (Fig. 3). Tropical ecoregions, with much higher species richness, had substantially stronger MPA effect sizes than temperate regions. Thus, the capacity for MPAs to enhance fish assemblage aesthetic value increases toward the equator with increasing regional species richness.

To disentangle the effects of MPAs we also tested the strength of each potential causal pathway one by one by blocking each pathway and quantifying the subsequent decrease in the overall MPA effect (see Methods). We expected MPAs to affect aesthetic values through fish diversity (taxonomic, phylogenetic, or functional), taxonomic composition, trophic composition, and benthic composition, with the sum of these effects approximating the total effect of MPAs. We found that the positive effect of MPAs on aesthetic values was primarily explained by increasing taxonomic diversity, with secondary contributions from phylogenetic diversity and taxonomic composition. Functional diversity, trophic composition, and benthic composition had nearly no contribution to the positive effect of MPAs on aesthetic value (Fig. 4). However, the individual contributions did not account for the entire MPA effect, meaning MPAs also enhance aesthetic values through other mechanisms.

227
228 **Family contributions to MPA effect**

229 The previous analysis revealed that taxonomic diversity was the main pathway by
230 which MPAs enhanced assemblage aesthetic values, with less contribution by
231 taxonomic composition. Although aesthetic values depend highly on taxonomic
232 composition (10), this is not the primary mechanism by which MPAs enhance aesthetic
233 values. To fully understand this result, we examined how the abundance of each family
234 varied between fished sites and no-take MPAs alongside families' average aesthetic
235 values. We found that among families that dominated in MPAs, there was a mix of
236 both high and low aesthetic values (Fig. 5). For instance, families considered less
237 beautiful such as Latridae (Trumpeters), Lethrinidae (Emperors), and Carangidae
238 (Jacks and Pompanos) were more dominant in MPAs than fished sites, which would
239 reduce aesthetic values below expectations based on species richness (Fig. 5).
240 However, MPAs were also characterized by some highly beautiful families such as
241 Acanthuridae (Surgeonfishes, Unicornfishes, and Tangs), Zanclidae (Moorish Idols),
242 and Pomacanthidae (Angelfishes) (Fig. 5). Hence, MPAs protect families considered
243 both beautiful and less beautiful, so altering taxonomic composition is not the primary
244 pathway by which MPAs increase aesthetic values of fish assemblages.
245

245

246 **Discussion**

247 Combining a global dataset of reef fish surveys with species-level aesthetic values,
248 we assessed the global status and drivers of aesthetic value in reef fish assemblages.
249 We found that aesthetic value is intrinsically linked to species richness but varies
250 among reefs of similar richness due to differences in species composition, i.e., from
251 species individual aesthetic effects. Some species have exceptionally high aesthetic
252 values (e.g. Blue-striped Angelfish, Comet) and their presence increases assemblage
253 aesthetic value well above the value expected from species richness alone. As
254 previously highlighted (10), these species cluster in families and phylogenetic lineages
255 found primarily on species-rich tropical reefs. This is a result of the evolution of color
256 and signaling (i.e. through social behavior and camouflage) being an important driver
257 of the richness of highly colorful fish species (15, 16). Consequently, assemblages with
258 aesthetic values substantially higher than expected occur in species-rich tropical
259 regions hosting families such as Chaetodontidae (Butterflyfishes), Pomacanthidae
260 (Angelfishes), and Acanthuridae (Surgeonfishes) associated with the exploitation of
coral reef niches. Temperature and productivity had the strongest effects on aesthetic value, which is unsurprising given their role in shaping biogeographic patterns in species richness (21). However, we also reveal that human gravity and protection status had important effects on assemblage aesthetic values, meaning conservation actions can enhance and preserve the overall beauty of shallow reefs, which is key for Nature's Contribution to People (NCP).

Although aesthetic values depend highly on species composition, we highlight that MPAs primarily enhance fish assemblage aesthetic values by increasing taxonomic diversity. People's perceptions of aesthetic value are known to depend strongly on species richness, with greater preference for more diverse assemblages (13). Assemblages with higher taxonomic and phylogenetic diversity are also more likely to host species considered exceptionally beautiful (Fig. 1). Yet, differences in taxonomic composition did not contribute strongly to the positive effect of MPAs on aesthetic value. Many heavily fished families such as Carangidae (Jacks and Pompanos), Sparidae (Sea Breams and Porgies), and Lutjanidae (Snappers) have low average aesthetic values, and these groups benefit from protection against fishing (22, 23). Although MPAs also host many beautiful families, the fact that protection benefits species considered less beautiful implies that differences in taxonomic composition are not the primary means by which MPAs enhance fish assemblage aesthetic values. This is an intuitive result, but it has important implications for conservation, because it demonstrates that MPAs can potentially serve multiple NCPs simultaneously. MPAs have the potential to boost populations of fished species, while simultaneously creating hotspots of beauty by supporting greater species diversity. We thus find important support for MPAs as broad conservation tools that can potentially benefit both fisheries and tourism simultaneously.

We also found that fish assemblage aesthetic value is higher on reefs with greater coralline algae and live coral cover. Coral cover creates complex three-dimensional habitat that can enhance fish species diversity through greater niche availability, as well as by increasing refuge for juveniles and prey species (24, 25). Many exceptionally beautiful species are also facultative or obligate coral users/feeders, notably butterflyfishes and wrasses (26). The positive relationship between live benthic habitat and aesthetic value is particularly important because it
means that managers can potentially enhance the beauty of reef fish assemblages by protecting living coral or through coral restoration, albeit at very localized scales. Unfortunately, the influence of benthic habitat also means that future coral mortality will strongly degrade fish assemblage aesthetic value, and major declines in overall reef beauty are to be expected in coming years. Surprisingly, benthic composition was not identified as a major pathway by which MPAs increase aesthetic value. However, it is important to recognize that this is the average result across all reefs surveyed worldwide and results are context-dependent at the local scale. In many cases, MPAs and openly fished sites have similar benthic composition, and MPAs may be deliberately placed in degraded areas to promote recovery (27). However, given that live coral cover boosts the aesthetic value of reef fish assemblages, managing for coral-dominated reefs is a key way managers can indirectly conserve aesthetic value, and placing MPAs around thriving coral is likely to produce the greatest overall aesthetic benefits (28).

Interestingly, countries' human development index had no influence on fish assemblage aesthetic value, while those wither greater fishing dependency hosted more beautiful assemblages. This is another important result, because it indicates that the ability to protect and enhance aesthetic value is not constrained by wealth or socioeconomic development. Wealthier countries could have greater capacity and infrastructure to create and enforce MPAs, yet both high and low-income countries can host exceptionally beautiful fish assemblages. Moreover, given the greater potential for MPA benefits in species-rich, tropical regions and the fact that human development and GDP decrease toward the equator (29), MPAs are promising tools to safeguard multiple NCPs in developing tropical countries where socio-economically viable or acceptable.

Although our results highlight the effectiveness of MPAs in conserving the beauty of the world's reef fish assemblages, only 7% of global coral reefs are partially or fully sheltered by MPAs (30). Moreover, previous studies have shown that most MPAs are ineffective due to illegal fishing, poor management, improper placement, or other issues (31). Given the low level of marine protection globally and the sensitivity of aesthetic value to environmental conditions and coral cover, the beauty of the world's reefs appears severely threatened, and immediate action is needed to
integrate aesthetic value into reef conservation and management plans. Biodiversity
worldwide is deteriorating at rates unprecedented in human history, yet humans
continue to depend on biodiversity for recreation, cultural inspiration, and material
goods. In response to this challenge, the Kunming-Montreal Global Biodiversity
Framework set out to ensure that by 2030 at least 30% of marine and coastal areas
of particular importance for biodiversity and ecosystem functions are conserved and
managed through protected areas and that 30% of degraded areas are under effective
restoration. Our findings suggest that investments in expanding and strengthening
MPAs or restoration will pay off particularly when effectively protecting reef systems
where both material (e.g., fish biomass for consumption and nutrition) and non-
material (aesthetic value) contributions can be secured.

Despite the importance of our findings, our study had several important
limitations. First, we only considered fish assemblages, which are one component of
reef ecosystems. Other organisms from hard corals and seaweeds to crabs and
nudibranchs all contribute directly to reef aesthetic values. Given the aesthetic appeal
of hard and soft corals themselves and the invertebrates that live within them, many
reefs with depauperate fish assemblages are likely to still be considered beautiful by
many people. Comparing patterns in aesthetic values across taxa should reveal
important insights into the covariation of beauty and whether beautiful fish species are
associated with beautiful species of coral or other organisms. We also removed
Pleuronectiformes (Flatfishes) and Syngnathiformes (Seahorses and Pipefishes) that
could contribute importantly to aesthetic value, particularly species like Pygmy
Seahorses, which are highly sought after by recreational divers (32). Additionally, by
only considering ray finned fishes (Actinopterygii) we miss the potential contributions
of sharks and rays, which also attract major tourism enthusiasm (11, 32). The
disproportionate presence of very large fishes within MPAs should also add to the
aesthetic values quantified here and should be incorporated into more comprehensive
estimation of people's perception of reef fish assemblages. People may also be drawn
to large schools of fishes, and although we weighted aesthetic value by abundance,
we cannot fully capture this effect, but this is an interesting avenue for future research.
It has also been demonstrated that some taxonomic groups are more likely to avoid
divers (33), and diver surveys do not account for uneven species detectability.
However, identical protocols are used in all RLS surveys and because human interest
in reef fish aesthetics is mainly through diving and snorkeling, our data truly reflect the potential human perception of beauty. We also acknowledge that appreciation of beauty can be subjective and vary from person to person depending on socio-cultural background (e.g. (34)). However, despite these potential differences, the studies on which we based our aesthetic evaluation (Tribot et al. (13), Langlois et al. (10)) found no effect of the public's socio-cultural background on their aesthetic preferences (despite using very large pools of more than 2000 and 13000 participants, respectively, although it should be noted that these surveys were only available in two languages). While we did not assess water clarity in this study, realized aesthetic appeal will also depend highly on visibility, and assemblages in clear, oligotrophic waters likely have greater potential for aesthetic value. Lastly, while we find that taxonomic diversity is the main pathway by which MPAs boost aesthetic value, the absolute difference in richness between fished sites and MPAs is low (~5 species on average), and evidence for MPAs enhancing richness at local scales is mixed (35). The causal pathways we tested did not fully account for the total effect of MPAs, and the importance of taxonomic composition may vary regionally. For instance, higher aesthetic value in MPAs could be driven by greater abundances of Acanthurids in one region but by Pomacanthids in another region, and a broad assessment of taxonomic composition may not fully account for this. The finding that benthic composition was not a primary pathway by which MPAs enhance aesthetic value is also likely influenced by the limited and simplified habitat date available for this study. Further study, particularly at local scales, is needed to verify whether MPAs have the potential to boost aesthetic value by harboring greater coral habitat.

By shifting towards a Nature's Contributions to People framework, we can recognize the critical role that aesthetic value plays in shaping our relationship with nature. This recognition can help us develop more comprehensive conservation and management strategies that recognize and enhance the full range of benefits that shallow reefs provide to people. By conserving and restoring the aesthetic value of shallow reefs, we can not only protect these ecosystems but also contribute to the well-being of the people that rely on them.
389 **Material and Methods**

390 **Reef Life Survey data**

391 We used fish assemblage data from the Reef Life Survey (RLS) – a monitoring program that uses standardized protocols to assess the composition and diversity of fishes, invertebrates, and benthos on shallow rocky and coral reefs worldwide (19).
392
393 During RLS surveys, trained divers identify, count, and estimate the sizes of all fishes encountered along two 50-m x 5-m blocks per transect. The fundamental survey unit in RLS is therefore the individual transect. Usually, two transects are surveyed within the same patch of reef (within 200 meters of each other), an area referred to as a site.
394
395
396
397
398 To maximize spatial coverage and avoid repeat surveys, we selected the most recent surveys conducted across RLS sites (2006-2019). This resulted in 7,013 surveys from 3,526 sites, encompassing 2,655 fish taxa. Fish taxa included unidentified families and genera, for which precise aesthetic values are not available. We then filtered out such taxa for which we did not have aesthetic information. The final list of species used in this study is based on Langlois et al. (10) and included the ray-finned fishes (Actinopterygii) from the RLS database, minus the orders Pleuronectiformes (14 species) and Syngnathiformes (31 species) that were removed because of their unusual morphologies that did not fit with the protocol used by Tribot et al. (5). After filtering the RLS fish assemblage data to species for which we had aesthetic information, we ended up with a total of 2,270 fish species.
399
400
401
402
403
404
405
406
407
408
409

410 **Aesthetic value of fish assemblages**

411 The next step was to estimate the aesthetic value of the 7,013 surveys in our dataset.
412
413
414 To do so, we relied on three previous studies. Tribot et al. (13) Conducted an online survey to directly quantify the aesthetic values of coral reef fish assemblages along a gradient of species richness. They then determined the influences of species richness and species composition on assemblage-level aesthetic value using a multiple regression model. In this model, species composition was included as the presence/absence of each species (0 or 1). This model showed that richness and species composition together explained nearly all variation (R²=0.94) in assemblage-level aesthetic value. The authors were then able to assess the contribution of each species to assemblage-level aesthetic value based on species individual regression coefficients, which they referred to as species “aesthetic effects” (i.e., how much
species presence or absence increased or decreased overall aesthetic value). In a separate study, Tribot et al. (5) also directly quantified individual species aesthetic values based on an online survey where respondents ranked individual species by their perceived beauty. Finally, Tribot et al. (13) showed that there was a strong correlation between species aesthetic effects at the assemblage level (coefficients from the regression model) and species individual aesthetic values from the survey rankings.

To quantify aesthetic values for our 7,013 RLS surveys, we therefore needed i) species richness and ii) species aesthetic effects. Species richness was simply calculated as the number of species per survey. However, because we were working backward (trying to predict assemblage-level aesthetic value), we did not have aesthetic effects for each species. Fortunately, because Tribot et al. (13) showed that species aesthetic effects were strongly correlated to species individual aesthetic values, we were able to estimate species aesthetic effects using species aesthetic values from Langlois et al. (10). Our procedure for estimating aesthetic value for each RLS survey was therefore i) convert species aesthetic values from Langlois et al. (10) to species aesthetic effects using the fitted model from Tribot et al. (13), ii) input species richness and species aesthetic effects into the regression model from Tribot et al. (13), i.e., Aesthetic Value ~ Sp. Richness + Sp. Aesthetic Effects. To account for uncertainty in the models fitted by Tribot et al. (13), we regenerated their original models as Bayesian regressions and used the resulting posterior distributions to generate predictions. For each parameter (i.e. regression coefficient), we extracted 4000 posterior draws, ran predictions, and took the median value of the resulting prediction distribution. However, while Tribot et al. (13) focused on species' occurrences only, we integrated species abundances, which can strongly influence assemblage-level patterns. To do so, we weighted species' aesthetic effects in the regression model by their log-transformed abundances. Thus, highly beautiful species would contribute less to assemblage aesthetic value if they are rare, etc. This procedure provided estimates of the total aesthetic value of all 7,013 surveys in the dataset.

We also calculated the expected aesthetic value based purely on species richness, which was simply calculated by inputting species richness only (without species aesthetic effects or abundances) into the model from Tribot et al. (13). Lastly, we calculated the deviation in total aesthetic value from the value expected based
456 purely on species richness, which was a measure of how much higher or lower each
457 assemblages' total aesthetic value was due to compositional effects (e.g., hosting an
458 abundance of exceptionally beautiful or non-beAUTIFUL species). To make this clear for
459 the reader, the total aesthetic value of each survey can be visualized as the value of
460 each dot along the Y-axis in Fig. 1b. Again, this is the value predicted by the
461 combination of species richness and species abundance-weighted aesthetic effects.
462 The aesthetic value expected based purely on species richness can be visualized as
463 the fitted white curve in Fig. 1b (the value predicted by species richness alone). The
464 deviation from the value expected based on species richness can be visualized as the
465 residuals between the fitted white curve and the actual blue dots in Fig. 1b (i.e., how
466 much higher or lower each dot is along the Y-axis in comparison to the fitted white
467 curve).

468

469 **Global maps**

470 RLS contains over 3000 sites, many of which are grouped within a few kilometers.
471 Visualizing these data on maps is difficult because many points are superimposed,
472 which can mask important values and trends. To generate a global map that accurately
473 depicted aesthetic value across sites without superimposition, we used Inverse
474 Distance Weighted spatial interpolation (IDW), following (36). IDW functions by
475 assuming that similarity between locations diminishes with increasing distance. We
476 first created a 1000x1000 cell raster grid and assigned existing sites to cells. We then
477 created a 150-km buffer around all empty cells and calculated the distance to existing
478 sites within the buffer. For each empty cell, the value was interpolated as the average
479 of all corresponding sites weighted by the inverse distance raised to the power of 1.6
480 (power values between 1 and 2 are recommended, with higher values decreasing the
481 influence of faraway points). This produced a 1000x1000 cell grid of aesthetic value
482 that was then mapped in base R with the center point of each grid cell plotted as a
483 point.

484

485 **Causal inference theory and Directed Acyclic Graphs (DAGs)**

486 It is often accepted that cause-and-effect relationships can only be derived from
487 randomized experiments where confounding variables are carefully controlled and
488 only variables of interest are manipulated (20, 37). Consequently, studies based on
489 observational data usually conclude that they can identify correlations, but that
490

correlations do not necessarily represent cause-and-effect relationships. However,
recent integration of causal inference theory within ecology and biodiversity science
has improved our ability to estimate cause-and-effect relationships from observational
data (20, 38).

491

492

493

494

Here we apply a recently emerging and widely applicable causal inference
framework, the Structural Causal Model (SCM) (20, 37), to determine the global
drivers of aesthetic value in reef fish assemblages. SCM uses directed acyclic graphs
(DAGs) to visualize the causal structure of a system under study, and subsequently
guide covariate selection required to answer causal queries. DAGs are used to
represent a researchers' assumptions about the causal structure of a system or
process under study, and must be justified, for example through published literature,
expert opinion, and/or scientific consensus. A DAG consists of a set of nodes
(variables) that are connected to each other by arrows. These arrows represent
hypothesized causal relationships between variables, pointing from cause to effect,
with causes preceding their effects. Importantly, DAGs must include all variables (both
measured and unobserved) required to depict a system or process under study,
including common cause variables, defined as variables that affect two or more
variables already included in a DAG (37)

495

496

497

498

499

500

501

502

503

504

505

506

507

508

During a series of meetings with reef ecology experts from around the world (BiodivERsA REEF-FUTURES consortium), we constructed a DAG depicting hypothesized cause-and-effect relationships between several anthropogenic and environmental covariates (SI Appendix, Fig. S5, Table S1-S2). First, we included all variables of interest (i.e., predictors), that is those for which we wished to assess the potential causal effect on fish assemblage aesthetic value. Next, we added common cause variables, which are variables that affect two or more of the variables already in the DAG. For example, although we were not interested in the causal effect of latitude on aesthetic value, we included it in our DAG as we hypothesized that latitude influences several variables of interest such as sea surface temperature, net primary productivity, and the human development index. We further scrutinized our DAG for potential missing links and added these links where necessary to depict the overall causal structure of our study system. We continued to discuss and modify the DAG until the majority of participants reached a consensus and were satisfied with the links and structure.

509

510

511

512

513

514

515

516

517

518

519

520

521

522
Once the DAG was finalized, we tested DAG-data consistency. Simply put, a DAG usually implies many independencies (e.g., X is independent of Y) and conditional independencies (e.g., X is independent of Y, given Z) that should be consistent with the observational data, given that both the observational data and DAG are representative of the data-generating process. If DAG-data consistency is ensured, this provides support for the asserted structure of the DAG itself. On the other hand, failed independencies indicate potential incorrect or missing links. The Dagitty framework, and the corresponding R package *dagitty* (39), conveniently tests DAG-data consistency via a series of partial correlations and allows users to specify the threshold for determining whether an independence has passed or failed. This threshold can be specified according to either p-values or Pearson correlation coefficients. We chose to use a Pearson correlation coefficient of 0.30, which most literature defines as a “weak” relationship (40). We opted against p-values because our database has over 7000 observations, and a Pearson correlation as low as 0.03 will yield a “significant” p-value with a sample size this large. Our DAG resulted in 170 conditional independencies. Initial DAG-data consistency tests revealed some failed independencies that we scrutinized and used to update and improve the DAG, i.e., the failures indicated we had missed important links that we then added. For example, we initially included a link from latitude to human development index (on the basis of (29, 41)), but no link from sea surface temperature to human development index. DAG-data checks indicated that a correlation existed between sea surface temperature and human development index even after accounting for their potential co-correlation with latitude. Thus, we added a direct link from sea surface temperature to the human development index, which is supported by (41–43). We continued to scrutinize and refine our DAG until all conditional interdependencies were validated.

It is important to note that results from our DAG-based approach to causal inference are dependent on the structure of the DAG and the underlying causal assumptions. Other DAGs are possible that could lead to slightly different results. However, a key strength of this approach is that our DAG is based on a synthesis of scientific literature and discussions among reef ecology experts. Additionally, our causal assumptions are transparently published in our DAG, and readers can directly see how our assumptions translate to our models and results. Thus, readers can appreciate our results in the context of these assumptions, whereas most analyses
556
557
558
559

in ecology do not follow a formal procedure for selecting control variables that remove
statistical biases, including confounding and overcontrol bias.

Anthropogenic and environmental factors

560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588

We extracted data for all available anthropogenic and environmental factors depicted
in our DAG. Full details on all variables, such as their scale, their source, and whether
they were included as a predictor or control is provided in SI Appendix, Table S1.
Briefly, survey depth was recorded in-situ during RLS data collection using a diver
depth gauge. Sea surface temperature (SST) data came from the NOAA Coral Reef
Watch Daily 5km Satellite Heat Stress Monitoring database (44), which provides
mean, minimum, and maximum daily SST values from 1985 to the present. We
selected mean annual SST for each site over the last five years. Degree Heating
Weeks (DHW) is an indicator of acute heat stress that combines duration (number of
weeks) and intensity (degrees C above local expected maximum monthly
temperature) and is expressed as number of heating weeks (1 week at 1C over
expected = 1 DHW). We selected mean annual DHW over the last five years. DHW
data also came from the NOAA Coral Reef Watch Database. Net primary productivity
(NPP) came from the Vertically Generalized Production Model that estimates NPP
using satellite-derived surface chlorophyll concentrations, daytime SST, and
photosynthetically active radiation and is available from 2002-2019 (45). We selected
the mean annual NPP for each site over the last five years. Human gravity is a proxy
for anthropogenic pressure expressed as a function of human population size and
travel time to a reef and examines the amount of human pressure within the
surrounding 500km of a reef (46). Marine Protected Area (MPA) status was classified
as either no-take, where all fishing and extraction activities are prohibited, restricted-
take, where certain fishing gears or species removals are permitted, or openly fished,
where no regulations are in place. MPA status in the RLS database is usually recorded
at the time of surveys, although in some cases has been added using boundaries from
local MPA maps or based on local advice. Human Development Index (HDI) is an
integrative metric that reflects the level of human health, education, and wealth within
a country, and ranges from 0 to 1 (low to high). HDI data (from 2017) came from the
United Nations Human Development Indicators database (https://hdr.undp.org/data-
center) Fisheries dependency data came from (47), which integrates quantitative
589

indicators of the contribution of coastal fisheries to a countries' 1) economy, 2) employment, or 3) food security.

590

591

Benthic composition is assessed during RLS surveys using photo-quadrats of the substrate. Twenty evenly spaced photos are recorded along the survey transect and later processed using SQUIDLE+ software (https://squidle.org/). Points are randomly overlaid on photographs and users identify the organisms the points land on into broad functional groups or to the lowest taxonomic resolution possible. This provides percent cover data for benthic organisms such as corals, macroalgae, sponges, etc. Here, we binned the groups into nine broader benthic groups that can be consistently applied globally: coral, algae, coralline algae, seagrass, other sessile invertebrates (e.g., ascidians, sponges), coral rubble, rock, sand, and microalgal mats (e.g., cyanobacteria, encrusting algae films). Although we worked at the global scale, benthic data were primarily available for tropical sites (benthic data had not been processed and analyzed for many RLS locations at the time of analyses). We therefore only considered benthic composition from tropical sites and statistical models were adjusted to account for this (see Bayesian Regression Models below). We defined tropical sites as those that experience minimum monthly sea surface temperatures of 17 °C (48). Although the tropical sea surface temperature limit is usually defined at 20 °C, we used this broader definition of tropical oceans, thereby including locations where species of tropical affinity are present (48). Because we were interested in overall benthic composition across all categories, we used an ordination analysis. We first applied an arcsine transformation (inverse Sine of square-root transformed values) to the proportional benthic cover data, which were highly right skewed. This is similar to a log transformation but more appropriate for proportional data, which are bounded between 0 and 1 (49). We then applied a Principal Component Analysis (PCA) to the resulting transformed data. We then used the first two PCA axes as indicators of overall benthic composition and included them as predictor variables in our drivers models. Unfortunately, benthic data were unavailable for 36% of tropical surveys. A typical approach would be to delete all surveys where benthic data are unavailable. However, this would throw away a substantial amount of valuable surveys where data for all other covariates are available. Secondly, removing missing data can induce bias, for instance if missing data share common characteristics, e.g., missing data come from the most remote locations. In this case, removing surveys with missing data would systematically bias subsequent analyses. Therefore, we chose to impute

592

593

594

595

596

597

598

599

600

601

602

603

604

605

606

607

608

609

610

611

612

613

614

615

616

617

618

619

620

621

622
623

missing values using random forest, through the R package missForest (50). Full
details for this procedure are detailed in the Supporting Information.

624

625

626

## Species traits

627

To calculate the functional diversity of reef fish assemblages, we compiled a trait
database of six categorical traits that describe the ecological lifestyle of the 2,270
species included in this analysis. The traits selected were body size, diet, trophic level,
gregariousness, substrate preference, and water column position. Body size was
categorized into six ordered size bins: 0-7 cm, 7.1-15 cm, 15.1-30 cm, 30.1-50 cm,
50.1-80 cm and >80 cm. Diet was classified into seven categories: piscivores,
invertivores, corallivores, omnivores, microphages and grazers. Trophic level was
binned into five ordered categories of >2.5, 2.51-3.0, 3.01-3.5, 3.51-4.0, > 4.00.

628

629

630

631

632

633

634

635

Gregariousness, which describes schooling behavior, was classified into three ordered
categories: solitary, pairing or forming small schools, and schooling. Substrate
preference was classified into four categories: sand, coral, rock, or water column.

636

637

638

Finally, water column position was categorized into four ordered categories: benthic,
benthopelagic, pelagic site-attached and pelagic-mobile (moving between reefs). Trait
data was sourced from Stuart-Smith et al. (51), FishBase (52) and expert knowledge.

639

640

641

These methods resulted in a trait database of 2,460 species, where 90.16% of the
data was complete (2,218 species). For the 9.84% of missing trait data (242 spp.), we
imputed the missing values using missForest. To increase imputation accuracy, we
included phylogenetic position captured by the first four axes of a Principal
Coordinates Analysis (PCoA) applied to the ray-finned fish phylogeny, fishtree (53, 54)
using pairwise cophenetic distances (55). For all imputations, body size and trophic
level were converted to their continuous forms. The trait database was then filtered to
2,270 species based on the availability of aesthetic value information for the analysis.

642

643

644

645

646

647

648

649

650

## Assemblage metrics

651

We assessed taxonomic composition at the family level through an ordination analysis.

652

We first calculated the average abundance of each family in each site and then log-
transformed and Hellinger-transformed the abundances. The log transformation
reduces skew and dampens the influence of extremely abundant taxa, while the
Hellinger transformation reduces the influence of double zero similarity in ordination

653

654

655
656

analysis. We then applied a PCA to the transformed abundance data and used the first four axes as indices of taxonomic composition among sites (SI Appendix, Fig. S6).

657

658

We also assessed trophic composition through an ordination approach. We first calculated community-weighted mean values of trophic guild categories for each survey using log-transformed abundances. This results in a matrix of proportions of abundance comprised of each trophic guild for each site, e.g., 40% invertivores, 30% herbivores, 20% planktivores, 10% piscivores, etc. We again applied an arcsine transformation to the proportions because they were highly skewed and then applied a PCA (SI Appendix, Fig. S7). We used the first four PCA axes as indicators of overall trophic composition.

659

660

661

662

663

664

665

666

Because assemblage-level aesthetic values were calculated using both species richness and abundance, we used taxonomic, phylogenetic, and functional entropy to assess assemblage diversity. We applied the common mathematical framework based on Hill numbers, which generalizes diversity measures to allow for sensitivity to species abundances based on the parameter $q$. For all entropy analyses, we set $q = 1$, which weights species, functional groups, or phylogenetic lineages by their relative abundances. This corresponds to the exponential of Shannon entropy for taxonomic diversity, Rao's entropy for functional diversity, and Faith's PD for phylogenetic diversity (56, 57). For functional diversity, the parameter $\tau$ was set to the mean functional distance between all pairs of species (56). Phylogenetic diversity was calculated following (57) using the phylogenetic fishtree from (54) with updates by (53).

667

668

669

670

671

672

673

674

675

676

677

678

## Bayesian regression models

679

680

681

682

683

684

685

686

687

688

689

To estimate the potential causal influences of anthropogenic and environmental factors on fish assemblage aesthetic value, we used a series of hierarchical Bayesian regression models. In following with causal inference theory, we did not build one single statistical model to assess the potential influences of all predictor variables simultaneously, but rather built one individual model for each predictor variable. For each predictor variable, we used the DAG to determine the structure of the regression model. To do this, we used the online DAGitty interface (http://www.dagitty.net/). We individually set each predictor variable to the “exposure” in the DAG with aesthetic value set as the “outcome”. We identified the minimum adjustment set necessary to obtain the total causal effect of each predictor variable (SI Appendix, Fig. S8). The minimum adjustment set provides the list of additional covariates that must be included
in a regression model to satisfy the backdoor criterion and account for potential statistical biases. For each individual regression model, although multiple covariates were often included, we extracted only the effect size for the variable of interest. For instance, when estimating the total causal effect of sea surface temperature (SST), the DAG indicated the minimum sufficient adjustment set must include latitude as a covariate. Hence our regression model for sea surface temperature took the basic form: aesthetic value ~ SST + Latitude, yet only the effect of SST was derived from this model, with latitude included purely as a statistical control. For some variables, the DAG provided multiple options for minimum adjustment sets. Among the candidate sets, we chose the set that i) included the fewest covariates, and/or ii) included covariates for which we believed that data quality was the highest (e.g. measured in-situ or derived from robust remote sensing products at high spatial resolution). This approach was chosen to minimize potential error associated with high collinearity among covariates or variable measurement error.

Because each regression model had the same response variable (aesthetic value), the same model structure was always used. Aesthetic value is bounded between 0 and infinity and is highly log-skewed (many small and few large values). The maximum entropy distribution for modeling such data is either log-normal (equivalent to Gaussian applied on log-transformed values) or Gamma (58). All models therefore used the log-transformed values of aesthetic value as the response with a Gaussian likelihood distribution and an identity link function. All models included transect nested in site and site nested in country as random effects to account for the hierarchical structure of the RLS data (transect surveys are nested within sites and those sites are nested within countries), to account for the non-independence and spatial autocorrelation of observations belonging to the same site or country, and to allow for predictor variables at different spatial levels. For instance, depth is measured at the transect level, temperature is measured at the site level, and HDI is measured at the country level. Random effects ensure these predictor variables are assessed at the relevant spatial scale within the model. Additionally, each model included climatic zone - tropical or temperate - as a fixed effect. This fixed effect was included for two reasons, i) to account for the large intrinsic differences that exist in aesthetic value and anthropogenic/environmental factors between tropical and temperate locations, and ii) to account for the fact that benthic composition data were not available in temperate locations. Setting all benthic composition observations to 0 in temperate locations and
including climatic zone as a fixed effect rather than a random effect prevented temperate observations from contributing to the overall effect size of benthic composition. The full model structure and Bayesian notation are shown in the Supporting Information.

We additionally ran a second set of models in which the effect of species richness was “removed”, meaning the contribution of species richness to the total causal effect of any given predictor variable was deleted. This was done simply by including species richness as a control variable within each regression model. This assessed how anthropogenic and environmental conditions contribute to the residual value of aesthetic value not accounted for by species richness, i.e., the degree to which aesthetic value is higher or lower than expected for a given level of species richness. This approach is synonymous to first running a model of aesthetic value ~ species richness and extracting the residuals of that model to be used in a second series of models, but is the proper method for achieving this goal (see (59) for an explanation on the misuse of residuals as data).

All models were built using the R package *brms* (60) and were run with four chains of 4000 iterations. As indicated in the equations above, the prior distribution for all slope coefficients and intercepts was a normal distribution with mean 0 and standard deviation 3. This distribution was chosen to be non-informative yet regularizing - constraining the potential slope and intercept values to a realistic range to improve model running time (58). However, priors have little influence with very large datasets and sensitivity tests confirmed that our results were robust to both weakly informative and flat priors (SI Appendix, Fig. S9).

After each model was run, the posterior distribution for the variable of interest was extracted and finally all posterior distributions were regrouped to generate forest plots. All models were run until sufficient effective sampling size was achieved and Gelman-Rubin diagnostic statistics converged to 1.0. Model performance was checked for each model using the function *check_model* in the R package *performance* (61), which provided diagnostic plots of i) posterior predictive check, ii) linearity, iii) homogeneity of variance, iv) influential observations, v) collinearity, vi) normality of residuals, and vii) normality of random effects. All model diagnostic plots are provided in SI Appendix, Fig. S10-19. Model checking revealed high Variance
756 Inflation Factors (VIF) for the SST, NPP, DHW, and HDI models due to multicollinearity
757 with latitude. Unfortunately, given our DAG, no alternative minimum adjustment sets
758 were available that did not include latitude as a control variable. Therefore, to test the
759 degree to which multicollinearity influenced model results, we ran sensitivity tests
760 where we reran the models without latitude as a control variable. In all cases, effect
761 sizes showed zero or very little change whether latitude was included, and the test
762 models had low VIF, indicating that multicollinearity had little influence on model
763 results (SI Appendix, Fig. S20-23).

764 Lastly, we tested alternative model approaches. While we used random
765 intercepts for sites nested in countries, an alternative approach would be to include
766 random slopes. We tested models with random intercepts for country and site and
767 random slopes for country (however sites could not be nested in countries in this
768 configuration within the brms modeling package). However, we found that some
769 models produced singular fits, meaning that the data didn't support meaningful
770 variability for the random effects and that these random slopes models were potentially
771 too complex for our data, which are sparse in some locations. Thus, random intercepts
772 models were more parsimonious in that they accounted for spatial patterns and
773 heterogeneity in the data without being overly complex and causing convergence
774 issues. We also tested models that included a spatial Gaussian Process to explicitly
775 model spatial dependence in the data. These models provided very similar results to
776 our original models, although with some variables having weaker or more uncertain
777 effect sizes. These models and their results are fully detailed in the Supporting
778 Information.

779

780 **Estimating geographical variation in the effect of no-take MPAs**

781 To examine regional variation in the effect of MPAs, we re-ran the regression model
782 for the total causal effect of MPAs with the same minimal adjustment set as previous,
783 but we allowed the slope (i.e., effect size of MPA) to vary by ecoregion. This way, we
784 obtained an estimated effect of no-take MPAs for each ecoregion. The full model
785 structure and Bayesian notation are shown in the Supporting Information. After running
786 the model, we plotted the relationship between latitude (i.e., average latitude of sites)
787 and the estimated effect size of no-take MPAs across ecoregions (Fig. 3).

788
789 Disentangling the influence of no-take MPAs

790 We hypothesized that MPAs could enhance fish assemblage aesthetic value by
791 boosting diversity (taxonomic, phylogenetic, or functional), containing benthic
792 substrates more likely to host beautiful species (benthic composition), containing
793 trophic groups characterized by beautiful species (trophic composition), or containing
794 beautiful species that are sensitive to human impacts (taxonomic composition). To
795 disentangle these potential pathways, we again used our DAG to run a mediation
796 analysis. We ran a series of regression models where we blocked each corresponding
797 pathway individually and calculated by how much the total causal effect of MPAs
798 declined. For instance, to calculate the contribution of taxonomic diversity to the MPA
799 effect, we re-ran the same model for the total causal effect of MPAs, but this time with
800 taxonomic diversity included as a control. We then extracted the new effect size of
801 MPAs and subtracted this from the original (by subtracting the second posterior
802 distribution from the first). We did this for each pathway and recorded the resulting
803 value, which was considered the pathway's contribution to the overall MPA effect.
804

804

805 Examining how families differ between fished sites and no-take MPAs

806 We examined how family' abundances varied between fished sites and no-take MPAs
807 (i.e. MPA effect for each family) alongside the average aesthetic value of each family.
808 We first calculated the average abundance of each family at each site. To calculate
809 the MPA effect per family, we then ran an individual regression model for each family.
810 These models all had the same structure and used a Gaussian likelihood distribution
811 applied to log-transformed abundance values. Because abundances were site-level
812 averages, values were continuous and not integers, therefore permitting a Gaussian
813 distribution rather than a Poisson or Negative Binomial. Each model also included
814 Ecoregion as a random effect, under the assumption that sites within Ecoregions were
815 non-independent, particularly fished sites and MPAs. (see Spalding et al. (62) for
816 details on the Marine Ecoregions of the World). Ecoregion was also used in lieu of
817 country in these models, as no country-level predictor variables were included, and
818 ecoregion provided a finer-scale and more geographically balanced resolution. The
819 full model structure and Bayesian notation are shown in the Supporting Information.
820 We then calculated the average aesthetic value of all species in each family. To
821 understand how different families contribute to overall aesthetic value in MPAs, we
822 then identified the families with the greatest positive or negative MPA effects (top 25%
823

highest absolute effect sizes) and plotted them in descending order and colored them by their average aesthetic values.

824

825

Finally, while our results focus heavily on the potential influence of MPAs on aesthetic value, we would like to clarify that our global scale analysis does not compare spatially paired (i.e. adjacent) MPAs and fished sites. While we do consider random effects for sites, countries, and ecoregions (which account for the spatial distribution and non-independence of sites), our analyses assess average patterns across all no-take MPAs in the RLS database. A truly paired analysis would require a high number of directly adjacent MPA and fished sites, which is unfortunately not available at the global scale.

826

827

828

829

830

831

832

833

## Acknowledgments

834

835

836

837

838

839

840

841

We thank members of BiodivERsA REEF-FUTURES consortium for insightful discussions during the elaboration of our DAG. We thank the many RLS divers, scientific and government collaborators who assisted with field surveys, and the many funders of RLS field activities. Data from RLS used in the analyses are managed through, and were sourced from, Australia's Integrated Marine Observing System (IMOS)—IMOS is enabled by the National Collaborative Research Infrastructure Strategy (NCRIS). Lastly, MM would like to thank his son Owen for taking reasonable naps during paternity leave, otherwise this paper would have never been finished. Sadly, Owen is now two years old...

842

843

844

## Funding

845

846

847

848

849

850

This work was supported by 2017–2018 Belmont Forum and BiodivERsA REEF-FUTURES project under the BiodivScen ERA-Net COFUND program with the French National Research Agency. RDS-S was supported by the Pew Fellows Program in Marine Conservation at The Pew Charitable Trusts and an ARC Future Fellowship (FT190100599). Reef Life Survey data management is supported by Australia's Integrated Marine Observing System (IMOS)—IMOS is enabled by the National Collaborative Research Infrastructure Strategy (NCRIS).

852

853

## Data availability
854 All data required to recreate the analyses and figures used in this study are available
855 on the GitHub repository https://github.com/mcleamj/aesthetic_value, which has been
856 archived as a publication release at https://zenodo.org/doi/10.5281/zenodo.11551780

857

858 Code availability

859 All analyses were performed in R, and all codes are available on the GitHub repository
860 https://github.com/mcleamj/aesthetic_value, which has been archived as a publication
861 release at https://zenodo.org/doi/10.5281/zenodo.11551780

862

863

864

865

866

867

868

869

870

871

872

873

874

875

876

877

878

879

880

881

882

883

884

885
886 **References**

887 1. R. Hill, et al., Nature's contributions to people: Weaving plural perspectives. *One Earth* **4**, 910–915 (2021).

889 2. S. Díaz, et al., Assessing nature's contributions to people. *Science* **359**, 270–272 (2018).

891 3. J. Garcia Rodrigues, S. Villasante, I. Sousa Pinto, Non-material nature's
892 contributions to people from a marine protected area support multiple dimensions
893 of human well-being. *Sustain. Sci.* **17**, 793–808 (2022).

894 4. T. C. Daniel, et al., Contributions of cultural services to the ecosystem services
895 agenda. *Proc. Natl. Acad. Sci.* **109**, 8812–8819 (2012).

896 5. A.-S. Tribot, J. Deter, N. Mouquet, Integrating the aesthetic value of landscapes
897 and biological diversity. *Proc. R. Soc. B Biol. Sci.* **285**, 20180971 (2018).

898 6. M. Soga, K. J. Gaston, The ecology of human-nature interactions. *Proc. R. Soc. B* **287**, 20191882 (2020).

900 7. N. Cooper, E. Brady, H. Steen, R. Bryce, Aesthetic and spiritual values of
901 ecosystems: Recognising the ontological and axiological plurality of cultural
902 ecosystem 'services.' *Ecosyst. Serv.* **21**, 218–229 (2016).

903 8. D. L. Stokes, Things we like: human preferences among similar organisms and
904 implications for conservation. *Hum. Ecol.* **35**, 361–369 (2007).

905 9. B. Martín-López, C. Montes, J. Benayas, The non-economic motives behind the
906 willingness to pay for biodiversity conservation. *Biol. Conserv.* **139**, 67–82 (2007).

907 10. J. Langlois, et al., The aesthetic value of reef fishes is globally mismatched to
908 their conservation priorities. *Plos Biol.* **20**, e3001640 (2022).

909 11. L. S. Waechter, et al., The aesthetic value of Brazilian reefs: from species to
910 seascape. *Ocean Coast. Manag.* **247**, 106882 (2024).

911 12. G. E. Southon, A. Jorgensen, N. Dunnett, H. Hoyle, K. L. Evans, Biodiverse
912 perennial meadows have aesthetic value and increase residents' perceptions of
913 site quality in urban green-space. *Landsc. Urban Plan.* **158**, 105–118 (2017).

914 13. A.-S. Tribot, et al., Species diversity and composition drive the aesthetic value of
915 coral reef fish assemblages. *Biol. Lett.* **15**, 20190703 (2019).

916 14. G. Strona, et al., Global tropical reef fish richness could decline by around half if
917 corals are lost. *Proc. R. Soc. B* **288**, 20210274 (2021).

918 15. N. J. Marshall, F. Cortesi, F. de Busserolles, U. E. Siebeck, K. L. Cheney, Colours
919 and colour vision in reef fishes: Past, present and future research directions. *J.
920 Fish Biol.* **95**, 5–38 (2019).
921
922
923

16. P. Salis, T. Lorin, V. Laudet, B. Frédérich, Magic traits in magic fish: understanding color pattern evolution using reef fish. *Trends Genet.* **35**, 265–278 (2019).

924
925
926

17. J. N. Kittinger, et al., From reef to table: social and ecological factors affecting coral reef fisheries, artisanal seafood supply chains, and seafood security. *PloS One* **10**, e0123856 (2015).

927
928

18. M. Spalding, et al., Mapping the global value and distribution of coral reef tourism. *Mar. Policy* **82**, 104–113 (2017).

929
930

19. G. J. Edgar, R. D. Stuart-Smith, Systematic global assessment of reef fish communities by the Reef Life Survey program. *Sci. Data* **1**, 1–8 (2014).

931
932

20. S. Arif, M. A. MacNeil, Applying the structural causal model framework for observational causal inference in ecology. *Ecol. Monogr.* **93**, e1554 (2023).

933
934
935

21. J. E. Duffy, J. S. Lefcheck, R. D. Stuart-Smith, S. A. Navarrete, G. J. Edgar, Biodiversity enhances reef fish biomass and resistance to climate change. *Proc. Natl. Acad. Sci.* **113**, 6230–6235 (2016).

936
937
938

22. S. J. Campbell, G. J. Edgar, R. D. Stuart-Smith, G. Soler, A. E. Bates, Fishing-gear restrictions and biomass gains for coral reef fishes in marine protected areas. *Conserv. Biol.* **32**, 401–410 (2018).

939
940
941
942
943

23. J. Jouvenel, D. Pollard, Some effects of marine reserve protection on the population structure of two spearfishing target-fish species, Dicentrarchus labrax (Moronidae) and Sparus aurata (Sparidae), in shallow inshore waters, along a rocky coast in the northwestern Mediterranean Sea. *Aquat. Conserv. Mar. Freshw. Ecosyst.* **11**, 1–9 (2001).

944
945
946

24. V. Komyakova, P. L. Munday, G. P. Jones, Relative importance of coral cover, habitat complexity and diversity in determining the structure of reef fish communities. *PloS One* **8**, e83178 (2013).

947
948

25. D. Coker, N. Graham, M. Pratchett, Interactive effects of live coral and structural complexity on the recruitment of reef fishes. *Coral Reefs* **31**, 919–927 (2012).

949
950

26. A. J. Cole, M. S. Pratchett, G. P. Jones, Diversity and functional importance of coral-feeding fishes on tropical coral reefs. *Fish Fish.* **9**, 286–307 (2008).

951
952

27. A. Abelson, et al., Expanding marine protected areas to include degraded coral reefs. *Conserv. Biol.* **30**, 1182–1191 (2016).

953
954
955

28. A. Friedlander, E. Brown, P. Jokiel, W. Smith, K. Rodgers, Effects of habitat, wave exposure, and marine protected area status on coral reef fish assemblages in the Hawaiian archipelago. *Coral Reefs* **22**, 291–305 (2003).

956
957
958

29. M. Kummu, O. Varis, The world by latitudes: A global analysis of human population, development level and environment across the north–south axis over the past half century. *Appl. Geogr.* **31**, 495–507 (2011).
959
960

30. Q. Zhao, et al., Where marine protected areas would best represent 30% of ocean biodiversity. *Biol. Conserv.* **244**, 108536 (2020).

961
962

31. G. J. Edgar, et al., Global conservation outcomes depend on marine protected areas with five key features. *Nature* **506**, 216–220 (2014).

963
964
965

32. C. Cater, "Perceptions of and interactions with marine environments: Diving attractions from great whites to pygmy seahorses" in *New Frontiers in Marine Tourism*, (Routledge, 2007), pp. 49–64.

966
967
968

33. Y.-M. Bozec, M. Kulbicki, F. Laloë, G. Mou-Tham, D. Gascuel, Factors affecting the detection distances of reef fish: implications for visual counts. *Mar. Biol.* **158**, 969–981 (2011).

969
970

34. C. Pugach, H. Leder, D. J. Graham, How stable are human aesthetic preferences across the lifespan? *Front. Hum. Neurosci.* **11**, 289 (2017).

971
972

35. G. Ramírez-Ortiz, et al., Reduced fish diversity despite increased fish biomass in a Gulf of California Marine Protected Area. *PeerJ* **8**, e8885 (2020).

973
974

36. R. D. Stuart-Smith, et al., Tracking widespread climate-driven change on temperate and tropical reefs. *Curr. Biol.* **32**, 4128–4138 (2022).

975

37. J. Pearl, Causal inference in statistics: An overview. (2009).

976
977

38. J. Runge, Modern causal inference approaches to investigate biodiversity-ecosystem functioning relationships. *Nat. Commun.* **14**, 1917 (2023).

978
979
980

39. J. Textor, B. van der Zander, M. S. Gilthorpe, M. Liśkiewicz, G. T. Ellison, Robust causal inference using directed acyclic graphs: the R package ‘dagitty.’ *Int. J. Epidemiol.* **45**, 1887–1894 (2016).

981
982

40. P. Schober, C. Boer, L. A. Schwarte, Correlation Coefficients: Appropriate Use and Interpretation. *Anesth. Analg.* **126** (2018).

983
984

41. J. D. Sachs, A. D. Mellinger, J. L. Gallup, The geography of poverty and wealth. *Sci. Am.* **284**, 70–75 (2001).

985
986
987

42. R. G. Newell, B. C. Prest, S. E. Sexton, The GDP-temperature relationship: implications for climate change damages. *J. Environ. Econ. Manag.* **108**, 102445 (2021).

988
989

43. G. H. Park, D.-Y. Kim, J.-S. Kug, The Role of SST in Explaining the Economic Growth of Indian Ocean Countries. (2021).

990
991
992

44. G. Liu, et al., Reef-scale thermal stress monitoring of coral ecosystems: new 5-km global products from NOAA Coral Reef Watch. *Remote Sens.* **6**, 11579–11606 (2014).

993
994

45. M. J. Behrenfeld, P. G. Falkowski, Photosynthetic rates derived from satellite-based chlorophyll concentration. *Limnol. Oceanogr.* **42**, 1–20 (1997).
995
996

46. J. E. Cinner, et al., Gravity of human impacts mediates coral reef conservation gains. *Proc. Natl. Acad. Sci.* **115**, E6116 (2018).

997
998

47. M. Andrello, et al., Global mismatch between fishing dependency and larval supply from marine reserves. *Nat. Commun.* **8**, 16039 (2017).

999
1000

48. V. Parravicini, et al., Global patterns and predictors of tropical reef fish species richness. *Ecography* **36**, 1254–1262 (2013).

1001

49. R. R. Sokal, F. J. Rohlf, *biometry* (Macmillan, 1995).

1002
1003

50. D. J. Stekhoven, P. Bühlmann, MissForest—non-parametric missing value imputation for mixed-type data. *Bioinformatics* **28**, 112–118 (2012).

1004
1005

51. R. D. Stuart-Smith, et al., Integrating abundance and functional traits reveals new global hotspots of fish diversity. *Nature* **501**, 539–542 (2013).

1006

52. R. Froese, D. Pauly, FishBase. (2010).

1007
1008
1009

53. J. Chang, D. L. Rabosky, S. A. Smith, M. E. Alfaro, An R package and online resource for macroevolutionary studies using the ray-finned fish tree of life. *Methods Ecol. Evol.* **10**, 1118–1124 (2019).

1010
1011

54. D. L. Rabosky, et al., An inverse latitudinal gradient in speciation rate for marine fishes. *Nature* **559**, 392–395 (2018).

1012
1013
1014

55. V. J. Deastiani, V. A. Bastazini, V. D. Pillar, Using phylogenetic information to impute missing functional trait values in ecological databases. *Ecol. Inform.* **63**, 101315 (2021).

1015
1016
1017

56. A. Chao, et al., An attribute-diversity approach to functional diversity, functional beta diversity, and related (dis) similarity measures. *Ecol. Monogr.* **89**, e01343 (2019).

1018
1019

57. A. Chao, C.-H. Chiu, L. Jost, Phylogenetic diversity measures based on Hill numbers. *Philos. Trans. R. Soc. B Biol. Sci.* **365**, 3599–3609 (2010).

1020
1021

58. R. McElreath, *Statistical rethinking: A Bayesian course with examples in R and Stan* (Chapman and Hall/CRC, 2018).

1022
1023

59. R. P. Freckleton, On the Misuse of Residuals in Ecology: Regression of Residuals vs. Multiple Regression. *J. Anim. Ecol.* **71**, 542–545 (2002).

1024
1025

60. P.-C. Bürkner, brms: An R Package for Bayesian Multilevel Models Using Stan. *J. Stat. Softw.* **80**, 1–28 (2017).

1026
1027
1028

61. D. Lüdecke, M. S. Ben-Shachar, I. Patil, P. Waggoner, D. Makowski, performance: An R package for assessment, comparison and testing of statistical models. *J. Open Source Softw.* **6** (2021).

1029
1030

62. M. D. Spalding, et al., Marine Ecoregions of the World: A Bioregionalization of Coastal and Shelf Areas. *BioScience* **57**, 573–583 (2007).
1031

1032 **Figure legends**

1033

1034 **Fig. 1. Evaluating the aesthetic values of Reef Life Surveys. (A) Aesthetic effects**
1035 were computed for the 2,270 fish species in the RLS dataset. Digitally created
1036 drawings of three sets of species with high, intermediate, and low aesthetic
1037 values are shown as examples (CC0 1.0). (B) Relationship between the number
1038 of species in each survey (n=7,013 surveys across 3,526 sites) and the estimated
1039 assemblage aesthetic values. The thick line represents the expected aesthetic
1040 values computed with species richness only (see Methods). Two surveys with
1041 the same species richness but high or low assemblage aesthetic values are
1042 highlighted in red. (C) Distribution of the aesthetic effects among the species
1043 present in the two highlighted assemblages in panel (B). Boxplots show the
1044 median, 50%, and 75% quantiles.
1045

1045

1046 **Fig. 2. Global distribution and drivers of fish assemblage aesthetic value. (A)**
1047 Global map of assemblage aesthetic value. (B) Forest plot showing the effect
1048 sizes, 50% credible intervals (thick lines), and 90% credible intervals (thin lines)
1049 for all predictor variables of assemblage aesthetic value. (C) Global map of the
1050 deviation from the expected aesthetic value given species richness. (D) Forest
1051 plot showing the effect sizes, 50% credible intervals (thick lines), and 90%
1052 credible intervals (thin lines) for all predictor variables of the deviation value. In
1053 both panel (B) and (D), environmental variables are plotted in light gray and
1054 anthropogenic predictors are plotted in dark gray.
1055

1055

1056 **Fig. 3. Geographic variation in the effect of no-take MPAs on fish assemblage**
1057 aesthetic value. Ecoregions are plotted on the global map and colored by their
1058 no-take MPA effect size as predicted from the varying slope model (see Methods
1059 for details). Effect sizes for each ecoregion are also plotted against latitude
1060 (average latitude among sites within each ecoregion).

1061

1062 **Fig. 4. Contributions of the potential pathways by which no-take MPAs enhance**
1063 fish assemblage aesthetic value. Forest plot showing the effect sizes, 50%
1064 credible intervals (thick lines), and 90% credible intervals (thin lines) for the total
1065 causal effect of no-take MPAs and the contributions of each pathway to the total
1066 causal effect of no-take MPAs. Red points indicate negative effects, blue points
1067 indicate positive effects, and color intensity is scaled to effect size.
1068

1068

1069 **Fig. 5. Family contributions to the effect of no-take MPAs on fish assemblage**
1070 aesthetic value. Forest plot showing the effect sizes, 50% credible intervals
1071 (thick lines), and 90% credible intervals (thin lines) for the influence of MPA
1072 status on family abundances. Only the families with the highest absolute MPA
1073 effects (top 25%) are shown. Families are colored by the average aesthetic value
1074 of their species.

1075