The purpose of the model is to predict coral population dynamics as a function of hydrodynamic (i.e. waves and cyclones) and thermal disturbances, grazing pressure, larval connectivity, sedimentation (i.e. sand import and export), interspecific competitive interactions, and benthic community diversity (species richness and functional diversity). Time series defining disturbance regimes, sand cover and the diversity and number of external coral larvae are imposed and need to be defined before launching simulations. The grazing regime is imposed but can also be determined by activating the feedback process linking reef rugosity (created by colonies) to herbivore fish density to grazing pressure. Patterns in species cover, colony size distributions, recruitment rates and rugosity are used to understand the model’s dynamics and its accuracy.

The model consists of grid-cell agents each representing 1 cm², so that the benthic community is represented at a scale of organization smaller than the colony (equivalent to that of a polyp, although polyp size varies among species by several orders of magnitude) (Figure 1). During a simulation, an agent can be temporally part of a coral colony (798 species), a patch of algae (i.e. macroalgae, allopathic macroalgae, Halimeda spp., turf, articulated coralline algae or crustose coralline algae), sand, or bare substratum. Each agent is characterized by 33 variables that describe where the agent is in space (its position is fixed), its species identity (i.e. one of the 798 coral species or six functional groups of algae) and related functional characteristics, its age, the colony’s planar area and identification number, if it is bleached or was grazed recently, et cetera (Appendix 2—table 1). Coral colonies and patches of algae are entities composed of multiple agents sharing the same variable values (except for their spatial coordinates) and changing their state simultaneously during certain processes. For instance, dislodgement is simulated by converting all the agents forming the dislodged colony into barren ground; a turf algae overgrowing a colony is simulated by converting the coral agents constituting the overgrown part into turf, but conserving the information about the colony (i.e. identification number, size, species, growth form).

The size of the reef and the length of a time step are changeable. We defined a 25 m² reef for our simulations (i.e. 250,000 agents), which is usually the scale at which benthic communities are assessed in detail (e.g. Holbrook et al., 2018; Torda et al., 2018). We defined a 6-month time step because the empirical data we used to calibrate the model were collected biannually. We acknowledge that six months is a coarse time step, potentially preventing the simulation of subtle dynamics, for instance changes triggered by mild disturbances. We opted for this time step considering that (i) corals grow slowly (<180 mm yr⁻¹), and (ii) their reproductive cycles, and (iii) thermal and hydrodynamic disturbance regimes are seasonal. It is, however, possible to define shorter periods (i.e., 3- and 4-month time steps) as time steps in the model.

The model estimates three-dimensional colony surface areas using geometric formulae (Appendix 2—table 5) to determine the number of larvae produced in each colony (Appendix 2: §7.2.1.1), and (optionally) the rugosity of the reef (Appendix 2: §7.1.2.2). The model also accounts for colony and algae heights in overtopping processes (Appendix 2: §7.5.2.2 and §7.5.3.2). Algae have constant heights and colony heights equal to the radius of the colony planar area, assuming the latter is circular (Appendix 2—table 18).

Each time step includes the following consecutive processes: (i) grazing—patches of agents are randomly selected and grazed until a certain proportion of the reef is reached, (ii) coral reproduction—locally and regionally produced larvae attempt to settle, (iii) thermal disturbance, which, if triggered, eventually causes colonies to bleach and/or die; (iv) dislodgement and fragmentation—the effect of waves and cyclones on certain colonies, (v) growth—each living agent, selected in a random order, attempts to convert its neighbouring agents within a certain radius to its own state, (vi) sedimentation—barren ground agents are converted to sand and vice versa until the desired sand cover is reached, (vii) algae invasion—the remaining ungrazed, barren-ground agents are converted into algae agents (Figure 2) (note that the process differs from species invasion). The order at which processes iii, iv, and v happen must be defined beforehand.

During each time step, the model exports response variables: the cover of each benthic group (coral, algae, and sand), the planar area of each colony present per species, the number of recruited coral larvae m⁻² species⁻¹, and the reef rugosity (in cases when rugosity-grazing feedback is activated). The first two variables are collected after processes iii, iv, and vii, and the third and fourth variable after process vii. There are six variables imported each time step; their values respectively determine the cover to be grazed, the intensity of waves or cyclone, and of thermal stress, the number of external larvae m⁻² entering the reef, the order that reproduction, bleaching, and wave or cyclone events happen, and the cover of sand to be achieved. We present a complete schedule that includes additional model-related processes in Appendix 2: §3.

Basic principles: the model combines agent-based, trait-based, and demographic approaches to simulate coral reef community dynamics in imposed environmental scenarios. The model captures fundamental principles in ecology: (i) biodiversity influences ecosystem resilience and (ii) functioning, (iii) disturbance regimes filter species and mediate interspecific competition, (iv) interspecific functional differences (or strategies) mediate competitive exclusion and coexistence, and (v) source-sink dynamics regulate species coexistence in metacommunities.

Emergence: The dynamics of the benthic community emerge from species traits and the imposed disturbance regime (waves, cyclone, thermal stress), larval connectivity, sedimentation, and grazing pressure intensity (which can also emerge from reef rugosity). Cell agents do not make decisions and their behavior results from imposed deterministic or probabilistic rules.

Interactions: Agents on the edge of coral colonies and patches of algae interact with one another when competing for space. The outcome of a coral-coral interaction is determined by its specific pairwise outcome probabilities—the probability of coral-algae interactions are the same for all coral species, and algae-algae interactions result in a stand-off except when competing against crustose coralline algae. Branching and plating species also have the capacity to overtop other colonies and algae depending on their size (Appendix 2: §7.5).

Stochasticity: The model draws success or failure outcomes each time a patch of algae is grazed, a larva attempts to settle and survive the first 6 months, an agent tries to convert another living agent (Appendix 2: §7.1, 7.2, 7.5), and a colony is thermally stressed (i.e. bleaching and bleaching-induced mortality; Appendix 2: §7.4). Each of these random events is based on a probability of success specific to the process and the species involved. During initialization, colonies are created and placed randomly in space; their sizes are drawn from right-skewed frequency distributions (Appendix 2: §5.1). Finally, grazing and larval settlement happen randomly in space.

Collectives: Collective behaviour of agents happens when a colony is (i) dislodged—agents sharing the same colony (i.e. coral and algae agents growing on a colony) are converted to barren ground (Appendix 2: §7.3), (ii) bleaches or dies—the coral agents of the colony are converted to a bleached or dead state, respectively (Appendix 4: §7.1, 7.2), or (iii) reproduces—the number of larvae or gametes produces by a colony depends on certain coral traits and the size and age of the colony (Appendix 2: §7.2).

Observation: Four types of data are collected during a simulation: (i) percentage cover of each taxon, (ii) number of recruits for each coral species m⁻², (iii) planar area of each colony species⁻¹, and (iv) optionally, the rugosity created by the coral colonies (Figure 2).

The initial composition of the benthic community (i.e. the cover of coral species, algae, barren ground, and sand) is defined by the user and is imported from a comma-delimited text file. The space is filled first by creating circular coral colonies randomly in space. The colony diameters are drawn from skewed distributions that we defined using empirical data (E. H. Meesters and R. P. M. Bak, personal communication, May 2017) and as a function of the trait colony maximum diameter. Circular patches of algae (314 cm²) are then created and the remaining agents are converted into barren ground and sand (Appendix 2: §5).

Predefined time series (recorded in the text files) of input data are used to define the environmental context of the reef (Figure 2). At each time step, the model imports values for the corresponding period of the (i) surface grazed (%), (ii) number of external larvae settling, (iii) intensity of waves of cyclones (in dislodgement mechanical threshold, a dimensionless measure of the mechanical threshold imposed by waves and cyclones), (vi) thermal stress intensity (in degree-heating weeks), and (v) sand cover (%).

We used data collected between November 2001 and July 2011 in three sites located in Martinique in the Caribbean: Fond Boucher (14° 39′ 21.07″ N, 61° 09′ 38.98″ W), Pointe Borgnesse (14° 26′ 48.74″ N, 60° 54′ 12.72″ W), and Ilet à Rats (14° 40′ 58.04″ N, 60° 54′ 1.18″ W). The data were collected biannually (once per dry and wet seasons) by the Observatoire du Milieu Marin Martiniquais (OMMM) for the program Initiative Française pour les REcifs COralliens (IFRECOR). These data describe the benthic, macroinvertebrate, and fish communities at the species or genus levels, as well as sand cover for each site and at each sampling time (Appendix 3—figure 1, 2). We downloaded values of degree-heating weeks for the corresponding location from the US National Oceanic and Atmospheric Administration data server ERDDAP (Environmental Research Division's Data Access Program; coastwatch.pfeg.noaa.gov/erddap) (Appendix 3—figure 3). We identified cyclone tracks using the National Oceanic and Atmospheric Administration Historical Hurricane Tracks website (coast.noaa.gov/hurricanes).

We modelled thermal stress by inputting at each time step the maximum degree-heating week value found for the corresponding period. We represented the intensity of hydrodynamic regimes by inputting values of the dislodgement mechanical threshold (Madin and Connolly, 2006). We imposed a constant value in the absence of cyclone and a lower value when Hurricane Dean affected the reefs in August 2007 (its intensity changed from Category 1 to 2 while passing over Martinique). We chose threshold values arbitrarily considering wave exposure and cyclone intensity. Because of this uncertainty, we defined three different hydrodynamic regimes that we included in the calibration procedure for each site (Appendix 3—figure 5).

To estimate the percentage of the reef grazed at each time step, we first defined models predicting grazing intensity (i.e. percentage cover maintained in a cropped state) as a function of herbivorous fish and urchin density (we did not activate the rugosity-grazing feedback process). We defined these models using the empirical data from Williams and Polunin (2001) and Sammarco (1980) for fish and urchins, respectively (Appendix 3—figure 6). We then used these models and the population densities of Acanthuridae spp., Scaridae spp., and sea urchins measured in the three sites to predict their respective grazing regimes. Finally, we defined three additional similar regimes of different intensities, which we included in the calibration procedure (Appendix 3—figure 8).

The model adjusts the amount of sand cover (i.e. by removing or adding sand patches) at each time step according to the observed cover measured in each site (Appendix 3—figure 4). Having no information about larval connectivity at the three sites, we set the number of larvae m⁻²to 700 during each reproductive time period (i.e. once a year). This number corresponds to our estimate of competent larvae arriving on a hypothetical reef 20 km from an upstream reef having a 50% coral cover (Appendix 2: §7.2.1.2). The number is realistic considering that the distance separating the three sites from other coral communities is lower, but the average coral cover in the French West Indies is on average <40% (Wilkinson, 2008).

We calibrated the model for each site independently. We selected 12 parameters, for which we defined between two to five potential values (Appendix 3—table 1). We defined an algorithm to explore the parameter space optimally. The algorithm first selects the centroid, the most extreme values, and the values situated at mid-distance between the centroid and the extremes. A simulation with each parameter value is launched and replicated five times. We measured the fit between the empirical and simulated cover time series using an objective function. The objective function measures the performance of a given run by calculating the Euclidian distance between the empirical and simulated cover time series (averaged over five replicates), averaged over all the taxa (Appendix 3: §3.2). Performance is thus a positive value, with smaller values indicating higher performance (lower difference between simulated and empirical values). The algorithm then selects the 10 runs providing the best performance and generates for each of them the five closest (using Gower’s distance metric; Gower, 1971) and untested parameter combinations. The algorithm then launches these new simulations and repeats the procedure once more.

To compare the performance of model runs to a null expectation, we generated a null distribution of performance values for each empirical dataset by randomizing cover values within each row and calculating the distance from the original datasets.

Coral colonies grew ipso facto at their species-specific growth rate at low population density. However, as space filled up, colonies began constraining each other spatially and their growth rates decreased until eventual stasis (Appendix 5—figure 4).

For a single coral population, the different patterns of recruitment observed among three functionally distinct species (Appendix 5—figure 6) results from the interaction of several factors: (i) individual colony fecundity determined by its planar area, species-specific polyp fecundity, corallite area (polyp size), growth form, sexual system, and mode of larval development; (ii) the distribution of colony size in the population, which depends on maximum colony diameter; and (iii) the amount of surface available for larval settlement. Weedy (Agaricia tenuifolia) and stress-tolerant species (Echinophyllia orpheensis) produced bell-shape recruitment patterns (Appendix 5—figure 8). Recruitment rate was initially low because populations were composed of small, low-fecundity colonies, but the rate increased as colonies grew and became more fertile (Appendix 5—figure 8). Recruitment subsequently decreased as space became saturated. In contrast, recruitment rate for competitive species (Acropora gemmifera) was initially high and only decreased as cover occupancy increased. This pattern is essentially due to a higher vegetative growth rate associated with a population initially composed of fewer but larger, more fecund colonies (Appendix 5—figure 8).

In both the Western Atlantic and Eastern Pacific communities (we compared the two functionally distinct coral communities in the rest of the analysis), the competitive species dominated the coral community under low wave exposure (Appendix 5—figure 10, 11, 14). The success of the competitive species was due mainly to two interacting processes—with a higher vegetative growth rate, competitive species (i) overcame free space before other species, and (ii) enhanced recruitment by achieving large colony size rapidly. Higher wave exposure reduced the cover of competitive species because colonies were dislodged at a certain colony size, which reduced recruitment rate and provided other species with more available space to grow and recruit.

In the Western Atlantic community, increased availability of space favoured the weedy species (Madracis pharencis) over the stress-tolerant species (Orbicella annularis), principally because of the former’s brooding mode of larval development (twice a year), faster growth rate, and high wave-resistance of its growth form (digitate). In contrast, species coexisted in the Eastern Pacific community under the highest-intensity disturbances (Appendix 5—figure 14). Both the competitive (Pocillopora elegans) and stress-tolerant species (Porites lutea) recruited more than weedy species (P. damicornis) due to their spawning mode of reproduction; spawning species received three times more larvae from the regional pool (Appendix 2: §7.2.1.2). However, the weedy species recruited twice as frequently and is slightly more aggressive than the other two species. The stress-tolerant species has a massive growth form, which conferred higher resistance to waves compared to the other two branching species. Nonetheless, this advantage barely compensated for its slower growth rate and lower colony fecundity.

Population(s) recovered to pre-disturbance cover after only one year, regardless of the intensity of the event. This recovery is faster than most dynamics observed in real reef systems and arises because we imposed a constant and high number of larvae (7000 m⁻²) coming from the regional pool. In reality, larval supplies are reduced because a strong bleaching disturbance would also affect the surrounding reefs (Hughes et al., 2019). Another reason for this outcome was that recruitment preceded growth in the model (Figure 2), which inflated the former process because more space was available for settlement.

Low larval connectivity influenced the two coral communities differently (Appendix 5: §5). In the Western Atlantic, the weedy species thrived under zero to moderate larval input (0, 66, 700 larvae m⁻²) while the other two species went locally extinct (Appendix 5—figure 17). The weedy species produced ready-to-settle larvae twice a year, while the other two species reproduced annually and only a portion of their larvae were able to settle because of their time to motility (Appendix 2: §7.2.1.1.d). In contrast, the stress-tolerant species dominated in the Eastern Pacific community (Appendix 5—figure 18) due to its higher wave-resistance compared to the other two branching species.

Under the highest larval connectivity (7000 and 35,000 larvae m⁻²), the competitive species dominated in both communities, principally because of their higher growth rates, spawning mode of reproduction, and their capacity to overtop smaller colonies.

Population dynamics were similar between the two communities (Appendix 5—figure 20, 25). The total coral cover corresponded approximately to the imposed percentage of reef grazed, and the remaining ungrazed part of the reef was occupied by algae (also, ungrazed coral agents potentially exist). We observed no hysteresis because we did not implement feedback processes. Turf dominated the algae community in all simulated grazing regimes, despite having the highest palatability among algae (Appendix 2—table 3). The success of turf was due to its much higher growth rate compared to other algae (Appendix 2—table 19).

Coral recruitment rates at the steady state were the highest under medium grazing pressure (50%) because under lower and higher grazing intensities, space was saturated by turf and coral colonies, respectively. Under the lowest grazing pressure, most of the colonies were ≤100 cm² in surface area (Appendix 5—figure 22) and coral populations were rescued by external larval input. At intermediate pressures (30% and 50%), the competitive species dominated the coral community mainly because of their higher external larval input compared to the brooding (weedy) species, and their higher growth rate than the stress-tolerant species. Having a high growth rate was particularly important under low grazing pressure because this trait compensated better for the cover lost in competition with turf, which wins all its interactions with corals in the model (Appendix 2—table 19).

Under higher grazing pressures (70% and 90%), there were more coral-coral and fewer coral-algae interactions, which changed coral species dominance. The Western Atlantic community was dominated by the weedy species, followed by the stress-tolerant species and the competitive species was competitively excluded, mainly because of its lower aggressiveness and highest vulnerability to waves (Appendix 5—figure 10). In the Eastern Pacific community, the stress-tolerant species dominated the coral community and slowly outcompeted the other two species mainly because of its much higher wave resistance (Appendix 5—figure 27).

We collected trait values from coraltraits.org (Madin et al., 2016a) and other sources from the primary literature (Appendix 1—table 1). We limited our analysis to zooxanthellate scleractinian coral species, because our main focus is on species forming typical tropical-reef habitats. We first assembled a total of 828 coral species after correcting for nomenclature using the World Register of Marine Species (marinespecies.org) as a reference. We then removed 30 species that had only categorical trait values available (e.g. growth form, model of larval development, sexual system) and for which we did not have phylogenetic information. We included the categorical traits ‘coloniality’ and ‘sexual system’ to increase the number of predictors in the random-forest imputation procedure, because these traits have been defined for many species and have different degrees of phylogenetic conservation (i.e. sexual system is conserved, whereas coloniality has evolved several times within different lineages; Baird et al., 2009b; Kitahara et al., 2010).

For most species, the outcome of direct interactions between coral colonies is determined using the probability of species-pair interactions from Precoda et al. (2017) (Appendix 2: §7.5.2). We also defined for each species an aggressiveness ranking value, which is used to determine if a colony can overgrow another species’ colony. Aggressiveness values are only used for the species not considered by Precoda et al. (2017).

We defined aggressiveness ranking values combining six lists of coral species whose interspecific dominance relationships were assessed (Appendix 1—table 2). We first ranked the species in each study based on the metrics the authors used. Because each list had at least one species in common with another list, we could combine the six lists and rank all species considered from the least to the most competitive. We used the iterative partial rank-aggregation pivoting algorithm (IPRAPA) developed and explained in detail by Swain et al. (2017) (the authors used the method to rank 110 symbiodinium phylotypes aggregated from 35 reports based on their thermotolerance). The IPRAPA implements the ranking method based on consensus-based, Borda-rank aggregation (developed in voting systems) and updates ranks by a pivot element (i.e. a species shared between input lists). When more than one species is present in both lists, the species with the least uncertainty is designated as pivot element. In case of equal uncertainty, we selected the first species of the first list. We repeated the process 10 times. Each iteration provides an updated ranking score from the previous iteration. Ranking scores are decimal values comprised between 0 and 100, which represent the lowest and highest aggressiveness, respectively.

We selected the final ranking scores of a particular iteration based on two metrics: (1) the percentage of species-pair associations that were reversed after the aggregation compared to the original ranking in each list (i.e. Ppa, a change in dominance to equality or vice versa was not counted as reversed), and (2) the Kendall’s τ, which is a measure of ordinal association between two quantities (τ = 0 if all the species-pair associations are discordant, and τ = 1 if they are preserved). Finally, we measured the match between the final global ranking and the individual ranking of each study with the Spearman rank ρ. By ranking species in each list separately using the ranking score, we can see that the procedure conserves most of the initial ranking in each list (Appendix 1—figure 1).

Related code

Traits_and_imputation/Rscripts/ aggressiveness.R (Carturan et al., 2020).

Appendix 1—table 2

References	Metrics	Original no. of taxa	Final no. of species	Remarks
Lang, 1973	Ranking number based on number of subordinates	27	22
Sheppard, 1979	CI	26	22
Logan, 1984	Number of subordinates	17	15	Paired-interactions showing opposite result between field and lab experiment were not considered
Dai, 1990	Classification in five categories based on CI	76	76
Abelson and Loya, 1999	Ranking number based on frequency of wins and losses and composition of losing and winning species	33	30
Connell et al., 2004	% of winning interactions	13	11
Total number of species:		147	116

Appendix 1—figure 1

Download asset Open asset

Traits_and_imputation/Rscripts/phylogeny_coral.R and phylogeny_coral_and_functional_traits.R (Carturan et al., 2020).

We combined phylogenetic information with species functional traits to improve the prediction of missing trait data. The reasoning is based on the following assumptions: (1) closely related species are generally more functionally similar because of phylogenetic conservation of traits and (2) species having similar values for certain traits might have similar values for other traits because of functional trade-offs (Darling et al., 2012).

The phylogenetic information is included in the form of eigenvectors obtained from doing a principal component analysis on the distance matrix representing the phylogenetic branch length separating each species in the phylogeny (Swenson, 2014a). We first calculated the distance matrix of each of the 1000 trees using the function cophenetic from the R package stats (R Development Core Team, 2017). We averaged the 1000 matrices to obtain a final distance matrix. We then did a principal components analysis using the function prcomp from the R package stats. We selected the first nine eigenvectors using the broken stick method and used them as predictors (Appendix 1—figure 2). We retained the nine principal components because each describe different levels of the phylogeny (Swenson, 2014b).

Appendix 1—figure 2

Download asset Open asset

Among the different statistical methods available to perform imputation of missing data in functional-trait datasets (Penone et al., 2014; Schrodt et al., 2015), we chose the random-forest algorithm (Breiman, 2001) because it can handle highly dimensional datasets, does not rely on distributional assumptions, and is particularly appropriate for modeling complex interactions and non-linear relationships among variables. It is one of the best-performing techniques and allows the inclusion of phylogenetic information (Penone et al., 2014). We used the R package missForest 1.4 (Stekhoven and Bühlmann, 2012) to do the imputation. We set the maximum number of iterations to 10 and the number of trees generated by iteration to 100 (default values). Performance is evaluated with the normalized root mean squared error for numerical traits and proportion of falsely classified entries for categorical traits (good performance leads to a value close to 0, and bad performance to a value close to one). The method stopped after five iterations. It produced little error (normalized root mean squared error = 0.0985; proportion of falsely classified = 0.0735), and conserved the shape of the trait distributions (Appendix 1—figure 3, Appendix 1—figure 4, Appendix 1—figure 5).

Related code

Traits_and_imputation/Rscripts/imputation_traits_missForest.R (Carturan et al., 2020).

Appendix 1—figure 3

Download asset Open asset

Appendix 1—figure 4

Download asset Open asset

Appendix 1—figure 5

Download asset Open asset

coralreef2/src/coralreef2/Agent.java (Carturan et al., 2020).

coralreef2/src/coralreef2/InputData/ FunctionalTraitData.java (Carturan et al., 2020).

Appendix 2—figure 1

Download asset Open asset

We decided to implement coral growth, spatial competition, recruitment, their response to thermal and hydrodynamic disturbances and grazing pressure because they are fundamental ecological processes that shape coral communities and their dynamics. Additionally, these processes have been well studied, and associated functional traits have been collected on numerous species.

Executing these processes simultaneously, like in reality, is not possible because grid-cell agents can only be activated one after another. We consequently had to decide upon the order at which these processes were executed. We placed grazing, thermal and hydrodynamic disturbances before growth so they would affect the growth and interactions between coral colonies and algae. We implemented the possibility to decide on the order of occurrence between disturbances and coral reproduction because the timing of incidence of these events can be critical for coral population dynamics. Once all the organisms grow and compete, sand is imported or exported to a desired % cover. We did not implement a process of sedimentation (associated to hydrodynamic disturbances, for instance) due to insufficient empirical data. We thought it necessary to potentially represent sand cover because sand can occupy a large proportion of the reef and is an unstable substrate that prevents organisms from growing and settling. We placed this sedimentation process between growth and algae invasion to secure a sand cover close to the desired value and to constrain the remaining available space for the last process. Finally, we implemented algae invasion because we assumed that available non-grazed surface would be occupied by an algae within a several-month period in a real reef.

Related code

coralreef2/src/coralreef2/CoralReef2 Builder.java and ContextCoralReef2.java (Carturan et al., 2020).

Appendix 2—figure 2

Download asset Open asset

The model combines three fundamental approaches: (i) a complex agent-based approach (Breckling et al., 2006; Grimm, 2019; Grimm and Railsback, 2005)—the dynamics observed at higher levels of organization (e.g. population percentage cover, recruitment rates, colony size distributions, rugosity) emerge from agent-related processes happening at the lower scales (e.g. conversion of a barren ground agent to a new coral recruit agent, dislodgement of a single colony); (ii) a functional trait-based approach (Madin et al., 2016b; McGill et al., 2006)—species diversity and dynamics are determined by mechanistically linked trait-process associations; (iii) a demographic approach (Edmunds et al., 2014; Tuljapurkar and Caswell, 1997)—the dynamics of a population depends on its demographic structure (e.g. colony size distribution) because the size of a colony influences its capacity to reproduce, compete and resist disturbances. We combined these three approaches by implementing evidence-based trait-process mechanistic associations at an appropriate spatial scale (e.g. single agent for a larvae recruiting, all the agents forming a colony during a bleaching event) and accounting for size-process interactions (e.g. proportion of fecund polyps in a colony, colony shape factor for dislodgement).

The model captures several fundamental principles: (i) species diversity influences ecosystem resilience (i.e. the diversity-stability relationship; Ives and Carpenter, 2007; McCann, 2000; Nyström, 2006) and (ii) functioning (i.e. the diversity-ecosystem functioning relationship; Brandl et al., 2019; Loreau, 2000), (iii) disturbance regimes shape communities by filtering species and mediating interspecific competition (Kraft et al., 2015a; Sommer et al., 2014), (iv) interspecific functional differences (or strategies) mediate competitive exclusion and coexistence (Kraft et al., 2015b; Vellend, 2010), (v) source-sink dynamics regulate species coexistence in metacommunities (Amarasekare and Nisbet, 2001; Loreau and Mouquet, 1999). Note that there are many mechanisms leading to species coexistence (e.g. niche partitioning, spatial heterogeneity, facilitation; Adler et al., 2013; Chesson, 2000) that we did not implement in the model.

All the model outputs (i.e. population % cover, number coral recruits m⁻², colony size distributions and eventually rugosity) emerge from the processes implemented and depend on the imposed environmental conditions (i.e. larvae connectivity, grazing pressure, hydrodynamics and thermal disturbance regime, sand input). The outputted percentage of sand cover (if included) results mostly from the imposed sand cover (but is it possible that not enough space is available to reach the desired cover). If the rugosity-grazing feedback process is activated (Appendix 2—figure 2), the percentage of reef grazed results from the sum of the % imposed (if the user decides to maintain a certain imposed grazing regime) and the % emerging from the complexity of the habitat created by coral colonies (see §7.1.2.2). Cell agents do not make decisions and their behavior results from imposed deterministic or probabilistic rules.

Not implemented.

Not implemented.

Not implemented.

Not implemented.

Not implemented.

Agents on the edge of coral colonies and patches of algae interact directly with one another when competing for space. The outcome of a coral-coral interaction is determined by its specific pairwise outcome probabilities (or aggressiveness if the probability is not available for the species pair), the probability of coral-algae interactions are the same for all coral species, and algae-algae interactions result in a stand-off except when competing against crustose coralline algae. Branching and plating species also have the capacity to overtop other colonies and algae depending on their size (see §7.5). Agent also directly interact by occupied space, which prevents potential new recruits to settle.

The model draws success or failure outcomes in several processes, based on probabilities that represent specific ecological phonemes or details. We used (i) grazing probabilities to represent the difference of palatability among algae functional groups (see §7.1.2); (ii) settlement probability to represent the differences of settlement suitability among substrate types (see §7.2.1.3); (iii) a surviving probability for the first few months of the life of new coral recruits in order to control recruitment rates without implementing the causes of new recruit mortality (see §7.2.1.3); (iv) coral-algae, algae-algae and published species-pair probabilities of interaction outcomes to represent the non-transitive and inconsistent nature of direct interactions between these organisms (see §7.5); (v) species-specific bleaching probabilities and a bleaching-induced probability of mortality to account for the complexity and variability of bleaching responses observed within populations (see §7.4); (vi) species-specific positively skewed density distributions of colony diameters to create coral colonies during the initialization of the reef (see §5.2.2). Finally, grazing and larval settlement and the placement of colonies and algae during the initialization of the reef happen randomly in space.

Coral colonies are collectives of coral agents emerging from the growth of individual coral agents (new recruits). Agents being part of a same colony (i.e. alive, bleached and dead coral and algae on a colony) behave collectively when dislodged (see §7.3.1.2.a); living coral agents behave collectively when bleaching, dying from bleaching (see §7.4) and reproducing (see §7.2.1.1). Plating and branching corals can overtop other colonies and algae if their colony is large enough (see §7.5.2.2.b).

The model exports data that ecologists typically collect on the field to describe and understand patterns of community dynamics: (i) the percentage cover of each taxon (collected initially, after hydrodynamic and thermal disturbances and at the end of each time step and recorded in output/PercentageCover/PercentageCover_...csv), (ii) the mean number of recruits for each coral species m⁻² (collected at the end of each time step and recorded in output/NumberRecruits/NumberRecruits_...csv), (iii) the planar area of each colony (collected initially, after hydrodynamic and thermal disturbances and at the end of each time step and recorded in output/ColonyPlanarArea/ColonyPlanarArea_... csv), and optionally (iv) the rugosity of the reef created by coral colonies and the associated % cover grazed due to the fish community (collected during the grazing process and recorded in output/RugosityCoverGrazed/RugosityCoverGrazed_...csv).

The principal objective of the model is to simulate community dynamics as a function of species composition in different environmental scenarios. The composition and initial % cover for each coral species and algae functional group is imported from data/Initial_benthic_composition.csv and communities can be selected using their unique communityNumber (or ‘Community Number’ in the GUI). It is possible to define a % cover of bleached and dead coral species. Once the communityNumber is selected and the corresponding % cover imported, the grid is filled by creating agents. First, coral colonies are created at random locations, ordering the species from largest to smallest maximum colony diameter (this avoids having small colonies being enclosed by larger ones). Colonies are created by selecting an empty cell from the list of all available cells. A colony radius is then drawn from a species-specific size distribution (see §5.2) and used to select all the neighboring cells present in the list. These cells are removed from the list, and a cell-agent is created in each one of them. Upon creation, each agent is associated to the context and grid, and is given its x and y coordinates; all agent forming a same colony share the same other variables (Appendix 2—table 1). The procedure continues until the cover occupied reaches the defined value and is repeated for each coral species. Similarly, patches of algae agents (of a 10 cm radius) are then created in the remaining available cells, one functional group at a time, until the inputted cover is reached. The remaining empty grid-cells are filled by creating barren ground agents and finally patches of sand are created by converting patches of barren grounds agents. The initial % cover of sand is imported from data/Disturbance_sand.csv and it is associated to a certain environmental scenario, which can be selected using its unique disturbanceScenarioNumber (or ‘Disturbance Scenario’ in the GUI; see §6).

Placing coral colonies randomly in space is probably a limitation as aggregation of colonies influence coral reef dynamics (Brito-Millán et al., 2019; Eynaud et al., 2016). Providing the option to define the initial spatial arrangement of the colonies could be part of future development of the model. So far, only Wakeford et al. (2008) implemented coral colonies spatial arrangements according to empirical data.

The number of agents created (i.e. the size of the grid) is defined by the parameters reef_height and reef_width (in cm).

coralreef2/src/coralreef2/CoralReef2 Builder.java

coralreef2/src/coralreef2/InputData/ BiodiversityData.java,

Disturbance_bleaching.java, Disturbance_cyclone.java, Disturbance_grazing.java,

Disturbance_larvalConnectivity.java, Disturbance_priority.java, Sand_cover.java (Carturan et al., 2020).

Empirical_datasets_Martinique_calibration/Rscripts/

Benthos_data_2001_2010_backup.R, cyclones_bleaching.R, Martinique_algae.R,

Martinique_corals.R, Martinique_dead_substratum.R, Martinique_herbivores.R,

Martinique_Sites.R, species_selection_model.R (Carturan et al., 2020).

Appendix 3—figure 1

Download asset Open asset

Appendix 3—figure 2

Download asset Open asset

Appendix 3—figure 3

Download asset Open asset

We manipulated the empirical data so it could be used in our 6-month time step model. Biological data were collected at the end of the dry and wet seasons, which respectively span from December to May and from June to November. Measures of benthic cover were used in the calibration by being compared to the simulated cover values. We adjusted temporally the empirical cover data when necessary, so each measure would fall exactly 6 months after the previous one.

Sand cover was used to determine the target percentage at each time step (Appendix 3—figure 4). Small patches of barren ground agents were converted to sand or vice versa in case sand needed to be added or removed.

Appendix 3—figure 4

Download asset Open asset

The value of degree-heating weeks to impose during a time step is the maximum DHW value recorded during the corresponding period (Appendix 3—figure 5).

The intensity of hydrodynamic regimes is defined by the dislodgement mechanical threshold (DMT), a dimensionless measure of the mechanical threshold imposed by waves and cyclones (Appendix 2: §7.3.1; Madin and Connolly, 2006). A lower value represents a more intense disturbance. For instance, Madin and Connolly (2006) estimated that cyclone Rona (category three) imposed a DMT approximately equal to 18 at the crest and 88 at the back of Lizard Island in February 1999. In comparison, they also estimated the DMT imposed by a ‘moderate cyclone’ at the same locations: around 38 and 162 at the crest and back of the reef, respectively. We had no information about the intensity of the hydrodynamic regimes other than the category of the hurricane Dean. Considering that the three sites are not directly in contact with the Atlantic Oceanic currents (Fond Boucher and Pointe Borgnesse are in the Caribbean side and Ilet à Rats is protected in a bight), we first arbitrarily attributed DMT values to the different disturbance regimes: 140 for waves, 120 for tropical storms (TS), 100 for hurricanes of category 1 (H1), 80 for H2, 60 for H3, 40 for H4 and 10 for H5. The sites having potentially different exposure to both wave and cyclone, we defined two additional hydrodynamic disturbance regimes by (1) subtracting 10 and (2) adding 10 to the DMT values at each time step (respectively cyclone_model1 and cyclone_model3 in Appendix 3—figure 5).

Appendix 3—figure 5

Download asset Open asset

We defined four models predicting the proportion of the reef maintained in a grazed state as a function of herbivore densities (Appendix 3—figure 6). Among the different species present in the datasets, we only considered Acanthuridae spp. and Scaridae spp. and sea urchins, as they are considered the most important herbivores in reefs (Steneck, 1988). We determined the fish grazing pressure using Williams and Polunin (2001)’s empirical data. They conducted field surveys in 19 Caribbean reefs and analysed the relationship between the percentage cover cropped (i.e. covered by either turf, crustose coralline algae or bare substratum but not by macroalgae) and the density of Acanthuridae spp. and Scaridae spp. present. (Other grazers such as the sea urchin Diadema spp. were not present in high enough abundance to influence their results.) We combined abundances of the two genera and established a least-square linear regression between the total fish biomass and the percentage cover grazed (we considered ‘cropped’ as equivalent to ‘grazed’ in our model). The fish densities measured in Williams and Polunin (2001)'s study do not reach the maximum values observed in the three Caribbean sites (Appendix 3—figure 2). In order to avoid predicting percentage values above 100, we defined two asymptotic models approximating the linear regression but plateauing respectively at 90% and 70% for higher fish densities (Appendix 3—figure 6A):

where S_G = the surface of the reef grazed and D_F = the density of herbivore fish (g.m⁻²). We excluded the possibility for the reef to be grazed at 100% because such a scenario is not realistic even at high fish densities (Paddack et al., 2006; Williams et al., 2001).

Similarly, we modeled urchin grazing using data from Sammarco (1980), who experimentally manipulated the density of the sea urchin Diadema antillarum (excluding other types of grazers) and analysed the relationship with algal cover (the author did not specify which functional groups the term includes) in a reef at Discovery Bay (Jamaica). We defined a least-squares linear regression model to predict the ‘percentage of reef grazed’ (obtained by subtracting the percentage cover of algae from 100) with the density of D. antillarum (Appendix 3—figure 6). We log-transformed urchin density to meet the assumptions of normality and equal variance of the linear model (Appendix 3—figure 7). Both urchin species D. antillarum and Echinometra viridis were present in the three Caribbean reefs. We determined the urchin grazing pressure by pooling their respective abundances, assuming a similar function of the two species. This could potentially be a limitation as the two species have little niche overlap (McClanahan, 1999).

The total surface grazed is obtained by adding up the surface respectively grazed by fishes and urchins, which provides the grazing regimes grazing_model1 and grazing_model2 in Appendix 3—figure 8. We defined regimes grazing_model3 and grazing_model4 by subtracting 10% and 20% to grazing_model2 respectively.

Related code

Empirical_datasets_Martinique_calibration/Rscripts/Prep_for_model.R, Prep_for_model_Output_comparison.R (Carturan et al., 2020).

Appendix 3—figure 6

Download asset Open asset

Appendix 3—figure 7

Download asset Open asset

Appendix 3—figure 8

Download asset Open asset

We calibrated the model for each site independently and then conducted a between sites comparison of the parameter values that maximized model fit. A common approach for model calibration and validation is to divide the empirical data into a training and a test set. The model parameters are calibrated with the training set and the model generalizability is evaluated with the test set. Alternatively, resampling procedures can be used to train and test the model on the whole dataset iteratively (e.g. k-fold cross-validation). The time series we used were too short to be split and still allow for proper parameter calibration and model testing. Alternatively, we could have calibrated the model with one site and tested the model on the other sites. But we feared that model testing would not be satisfactory because of the paucity of the data describing the environmental context of the sites. By calibrating the model with the three datasets independently and then comparing the selected parameter values, we can dissociate between the parameters that can be generalized from those that need to be adjusted to a specific site. Alternatively, between-site differences in the parameter values offer the opportunity to question how relevant is the implementation of the related processes, if and how the implementation could be improved and if additional processes should be considered.

All the possible combinations of parameter values are generated. Each parameter combination is interpretable as a point, whose parameter values are its coordinates in the parameter space. The algorithm we defined first selects the centroid, the most extreme parameter points, and the points situated at mid-distance between the centroid and the extreme points (1). A simulation with each parameter point is launched and replicated five times for each of the three Caribbean sites (2). The fit between the empirical and the simulated cover time series is measured using an objective function. The algorithm then selects the ten points providing the best performance and samples around each of them to select the five closest (and not yet tested) points in the parameter space (3). Steps (2) and (3) are repeated one more time. Distance between points in the parameter space is measured using the Gower’s distance metric (Gower, 1971).

The objective function measures the performance of a given run by calculating the Euclidian distance between the empirical (coverEmp) and simulated cover time series (coverSim, averaged over five replicates), averaged over all the taxa:

where j = a specific taxon, Ntaxa = the last taxon, t = a specific time step, t_end = the last time step. Smaller values indicate better performance, zero being the minimum value and indicating perfect match.

In order to compare the performance of the runs to a second reference value, we generated for each empirical dataset a null distribution of performance values. These values are measures of performance between the empirical datasets and a randomized version of themselves. We randomized by rows (i.e. time periods) to keep the total taxa cover identical to the real data.

Being limited by the number of simulations we could perform in a reasonable time, we decided to calibrate the model in several rounds of simulations, each round being composed of one or several of the following steps: (1) select a limited number of parameters, define several possible values (instead of ranges for continuous parameters) and generate all the possible combination of parameter values; (2) launch the previously described sampling algorithm; (3) analyze the results of the calibration to identify the parameter values providing the worst performance; (4) analyze the simulated time series of the runs having the best performance to identify the processes and associated parameters that could improve the performance of the model and decide on the new values to add; (5) generate a new set of all the possible combinations of parameter values with the new parameter values and without those being associated with the worst performance; finally repeat (2 , 3) and (4) and eventually (5) until results are satisfactory. Eventually (7), launch additional simulations with specific parameter values selected by the modeller. This ‘hands-on’ calibration approach allows exploring specific parts of the parameter space potentially omitted by the algorithm or influencing specific processes that the modeller estimates to be relevant to act upon after having analysed previous results of the calibration. This human-directed search served to accelerate and guide the process, providing a faster convergence to the best parameter values.

We conducted in total four rounds of simulations. In round one, we selected only the parameters related to the different disturbance regimes (i.e. bleaching, cyclone and grazing models) and those whose value could not be found or estimated from the literature (i.e. growth rate reduction interaction, otherProportions, prob cover crustose coralline algae, ratio overtop colony, Appendix 3—table 1). We limited the number of possible parameter values to a maximum of five and generated the 4860 possible parameter points, from which our algorithm sampled 357 points for each site. This first round of simulations showed that none of the fittest parameter points at one site were sampled in another site. This could have been due to either the sites having different intrinsic characteristics leading to different local optima, or our algorithm was omitting certain parts of the parameter space. We consequently selected 40 of the parameter points providing best performance in each site and launched simulations with the subset of these points that were not run in a given site. This resulted in implementing round two with 41, 49 and 38 additional parameter points for Fond Boucher, Pointe Borgnesse and Ilet à Rats, respectively. None of these additional runs showed a better performance (Appendix 3—figure 9), confirming that the divergence in certain parameter values we observed between sites are due to intrinsic environmental differences. Analysis results at this stage revealed a consistently poorer performance of the simulations implementing values 0.01 and 0.001 of the parameter otherProportions (Appendix 3—figure 10). We consequently excluded these values in the next simulations. In addition, inspecting surface covers of the best fitted runs in each site revealed a poor match of the populations of algae. Grazing being one of the most important processes controlling these populations, we added more values to algae grazing probabilities. We then generated 2592 additional combinations of parameter values from which the algorithm sampled 202 of them (round three). The analysis of the simulated cover at this stage revealed a consistent incapacity of crustose coralline algae to maintain its population in all three sites. We decided to improve its persistence by enhancing its competitiveness against other algae (i.e. by adding the value 0.1 to prob cover crustose coralline algae) and reducing its palatability (i.e. by adding the value 0.05 to prob grazing crustose coralline algae). We launched the fourth round of simulations by generating 72 new parameter points with these new values and a subset of the other parameter values that consistently provided a better performance (i.e. bleaching model: 3; cyclone model: 1, 2, 3; grazing model: 3, 4; growth rate reduction interaction: 2, 4, 8; ratio overtop colony: 2; prob grazing macroalgae: 0.7; prob grazing allopathic macroalgae: 0.3; prob grazing Halimeda: 0.3, 0.5; prob grazing articulated coralline algae: 0.7). We generated a total of 7716 different parameter combinations from which we sampled 672, 680 and 669 parameter points for Fond Boucher, Pointe Borgnesse and Ilet à Rats respectively.

Performance values are bounded between 28 and 10 (where lower values are better) and are all below the lower 95% confidence limit of the null distribution of performance values (Appendix 3—figure 9). This shows that despite the model complexity, the high number of parameters, and the uncertainty around them, the model outputs population dynamics significantly closer to the real data as compared to randomly generated ones. The best performance values converged toward 10 in the three Caribbean sites (minimum values with standard error are 10.93 ± 3.677, 10.89 ± 2.872, 10.39 ± 3.119 for Fond Boucher, Pointe Borgnesse and Ilet à Rats, respectively) and were obtained during round three and four of the calibration (Appendix 3—figure 9). The first 257 points in round one were sampled uniformly in the parameter space, and expectedly show the highest ranges of performance values. The ranges from runs 258 to 357 narrow down drastically as sampling happened around the points providing best performance but the lower limits stayed constant. In round two, running in one site the combinations of parameter values that provided best performance in other sites increased the lower limit slightly for Fond Boucher, importantly for Pointe Borgnesse, and increased the higher limit for Ilet à Rats. The addition of alternative values for the probabilities of grazing the different groups of algae in round three lowered the lower range limits in all three sites. Finally, our attempt to support crustose coralline algae populations by increasing its competitiveness and reducing its palatability in round four, did not change further the range of performance values (Appendix 3—figure 9), as crustose coralline algae cover remained quasi absent even in the runs having the best performance (Appendix 3—figure 14, Appendix 3—figure 15, Figure 3).

For six parameters, the same value was the most commonly observed among the best fitted runs in the three sites: bleaching model (2), grazing model (4), otherProportions (0.0001), ratio overtop colony (2, Appendix 3—figure 10), and to a lesser extent prob grazing macroalgae (0.7) and prob grazing articulated coralline algae (0.7, Appendix 3—figure 11).

Three parameters have contrasting values between sites: growth rate reduction interaction (eight for Fond Boucher and Ilet à Rats and two for Pointe Borgnesse), prob grazing allopathic macroalgae (0.3 for Pointe Borgnesse and either 0.3 or 0.5 for the other two sites), prob grazing Halimeda (0.5 for Pointe Borgnesse and 0.3 for Ilet à Rats). The differences observed for prob grazing allopathic macroalgae and prob grazing Halimeda (Appendix 3—figure 11) could suggest that the accurate value for both parameters lies between 0.3 and 0.5. The small growth rate reduction value calibrated at Pointe Borgnesse compared to the two other sites could be due to the difference of turf abundance. Probably less turf-coral direct interactions occur at Pointe Borgnesse because turf is less abundant and shares its dominance with allopathic macroalgae (Figure 3; Appendix 3—figures 14 and 15). Turf has the fastest growth rate among algae and coral species and is the most competitive algae when competing with corals. As a result, a higher growth rate reduction value might be necessary for coral populations to persist in Fond Boucher and Ilet à Rats compared to Pointe Borgnesse. These differences in the values providing best performance for these three parameters could be caused by contrasting algae growth rates between sites due to either different species composing the functional groups or differences in grazing regimes that we did not capture or other processes not implemented (e.g. effects of nutrient concentrations on growth rate).

Lastly, values considered for three parameters did not influence performance: cyclone model, prob cover crustose coralline algae and prob grazing crustose coralline algae. Two reasons might explain why none of the three hydrodynamic disturbance regimes simulated influence performance: (i) the coral species present were either too resistant or (ii) they did not reach a colony size large enough to be dislodged. Colony dislodgement due to cyclones and waves is determined from colony planar area and growth form (Appendix 2: §7.3.1). Several of the dominant species in the three sites are massive (i.e. P. astreoides, O. faveolata, O. annularis, O. franksi) and could withstand even the hardest regime (i.e. model 3 in Appendix 3—figure 5). The other relatively abundant species have more vulnerable growth forms and their fragility increases with colony size (i.e. branching: M. mirabilis, P. furgata; digitate: M. decactis; laminar: A. agaricites). Possibly, spatial competition in the model prevented these species from reaching colony sizes large enough to be dislodged even under the most intense regime.

Analyzing the distribution of the combinations of parameter values with their associated performance in the parameter space reveals one unique global optimum in each site as the parameter combinations having the highest performance are clustered together (Appendix 3—figure 12). The parameters that influence performance the most in all sites are grazing model, and otherProportions, growth rate reduction interaction (respectively arrow number 3, 5 and 4 in Appendix 3—figure 12). The probabilities of grazing the different algae clearly dissociate runs generated during round one and two from the ones generated during rounds three and four: in the latter, higher values of prob grazing macroalgae and lower values of prob grazing allopathic macroalgae, Halimeda, articulated coralline algae and crustose coralline algae were implemented (respectively arrow number 8, 9, 10, 11 and 12 in Appendix 3—figure 12). When compared all together, the combination of parameter values providing best performance in sites Fond Boucher and Ilet à Rats overlap (Appendix 3—figure 13). Pointe Borgnesse performance optimum is distant from the other sites due to opposite calibrated values for growth rate reduction interaction (Appendix 3—figure 10, arrow number four in Appendix 3—figure 13).

Appendix 3—figure 9

Download asset Open asset

Appendix 3—figure 10

Download asset Open asset

Appendix 3—figure 11

Download asset Open asset

Appendix 3—figure 12

Download asset Open asset

Appendix 3—figure 13

Download asset Open asset

Appendix 3—figure 14

Download asset Open asset

Appendix 3—figure 15

Download asset Open asset

Our goal was to define an intrinsic (i.e. independent of environmental conditions) index of bleaching susceptibility expressed as a function of bleaching resistance traits. We established a statistical model between an empirical measure of bleaching susceptibility (dependent variables) and functional-traits (independent variables) on a subset of the 798 coral species for which the dependent variables had been measured. We then used the model to predict the bleaching susceptibility index value for the 798 coral species.

The dependent variable we selected is the bleaching response index (taxon-BRI), which represents the species-specific average percentage cover that bleached or died during an event. It was obtained for 374 taxa (304 if only considering species level values and after correcting for species names) from 2036 records concerning 316 sites, between 1982 and 2006 by Swain et al. (2016b). 65% of the Taxon-BRI value is determined by intrinsic factors (i.e. coral biology), 6% by extrinsic factors (i.e., environmental conditions) and 29% by measurement uncertainty (Swain et al., 2016b).

Coral bleaching is a complex process involving numerous resistant traits. We initially selected five traits based on the number of species for which the traits were measured and on the degree to which they are implicated in the bleaching process: (i) colony maximum diameter, (ii) growth rate, (iii) microscopic reduced scattering coefficient, (iv) growth forms and (v) corallite area (Carturan et al., 2018). We used the imputed traits dataset in order to avoid having missing predictor values.

We generated numerous beta regression models, defined the confidence set (i.e. the subset of models being the most supported by the data), averaged the latter and predicted the index values for the 798 species.

We first converted the taxon-BRI’s interval from [0100] to]0,1[. We then logit-transformed the taxon-BRI_]0,1[, as appropriate when modeling a proportional dependent variable (Warton and Hui, 2011). We log-transformed the numerical traits in order to reduce the skewness of their distributions and improve linearity of their association with taxon-BRI_]0,1[ (Appendix 4—figure 1).

The independent variables do show significant covariation (Appendix 4—figure 2): colony maximum diameter with corallite area (Spearman r_s = −0.21, p<0.001), growth form (Spearman r_s = −0.18, p<0.01), growth rate (Spearman r_s = 0.25, p<0.001) and microscopic reduced scattering coefficient (Spearman r_s = 0.33, p<0.001); corallite area with growth form (Spearman r_s = 0.55, p<0.001), growth rate (Spearman r_s = −0.73, p<0.001); and growth form with growth rate (Spearman r_s = −0.75, p<0.001). These covariations need to be considered when defining the full model (see §1.3). (We produced Appendix 4—figure 2 using the pairs.panels function from the R package psych 1.8.3.3; Revelle, 2017).

Appendix 4—figure 1

Download asset Open asset

Appendix 4—figure 2

Download asset Open asset

Proportional data often display asymmetric distributions, which contradict the normality assumption required in classic regression models. Contrastingly, beta regression models allow for the prediction of response variables having any unimodal distribution in the interval]0,1[ (i.e. ‘beta’ distributions; Ferrari and Cribari-Neto, 2004). Beta density distributions are defined by two parameters: the mean μ and the precision parameter ϕ and are mathematically expressed as Β(μ,ϕ). Obtaining the final beta regression model requires following several steps: (i) define the ‘global’ model, (ii) remove the observations (species) that are too influential on the model parameters (based on the Cook’s distance), (iii) chose the appropriate link function for μ, (iv) chose the appropriate link function for ϕ, (v) define the confidence set of models (i.e. data dredging), and eventually (vi) average these models. These steps are described in detail below:

Removal of influential observations

We used the Cook’s distance (Cook, 1977) to identify observations (i.e. species) being the most influential on the model parameters. Observations with a Cook’s distance superior to one should be removed (Johnson and Omland, 2004 and reference therein). We conserved all the observations because none of them had a large Cook’s distance (Appendix 4—figure 3). We calculated the Cook’s distance using the cooks.distance function from the package stats (R Development Core Team, 2017).

Appendix 4—figure 3

Download asset Open asset

Selection of the link function for μ

Different link functions can be used to model the distribution of the response variable (e.g. logit, loglog, etc.). We created a global beta regression for each of the link functions available in the function betareg (Appendix 4—figure 4). We then selected the link function based on (1) fit maximization—assess with the adjusted R² and the Akaike information criterion (AIC)—and (2) the distribution of the residuals (Ferrari and Cribari-Neto, 2004; Zuur et al., 2009). Johnson and Omland (2004) argued that AIC should be used over the adjusted R² in model selection because the latter does not account for the complexity of the model (i.e. the parsimony principal). However, in our situation, all the model candidates have the same complexity (i.e. number of predictors). It is consequently relevant to also consider this measure of fit. We assessed the homoscedasticity of the residuals by comparing Pearson residuals against predicted values (Appendix 4—figure 5). Finally, we assessed the distribution of the residuals by plotting the Pearson residuals against normal quantiles (Appendix 4—figure 6).

Each version of the global model satisfies the assumptions of homogeneity (Appendix 4—figure 5) and normality (Appendix 4—figure 6) of the residual distribution and no model seems to perform better in these aspects. We chose the model implementing the cloglog link function as it provides the highest fit (R² = 0.23; Appendix 4—figure 4).

Appendix 4—figure 4

Download asset Open asset

Appendix 4—figure 5

Download asset Open asset

Appendix 4—figure 6

Download asset Open asset