Fast-growing species shape the evolution of reef corals
Fossil dataWe downloaded all fossil occurrences recorded for the order Scleractinia at the species level from the Paleobiology Database (PBDB – paleobiodb.org; accessed on 3 August 2021). This is the most comprehensive repository for palaeontological data in reef corals to date. Due to the nature of the data, no ethics approval was required. To minimize identification issues, we excluded taxa with uncertain generic and species assignments (i.e., classified as aff. and cf.) and only selected species that had accepted names. We also selected the variables classification and palaeoenvironment from the output options to facilitate taxonomic and environmental filters applied in downstream analyses. The full dataset consisted of 24,011 occurrences across 4235 species, spanning over 250 Myr of coral evolution from the Triassic to the present. Although our focus here lies on the Cenozoic, we used the complete fossil dataset (i.e., including all of the occurrences) to have estimates of the diversification dynamics in scleractinian corals throughout the whole timespan of their evolution.Evolutionary ratesWith the full palaeontological dataset, we estimated evolutionary rates through time in scleractinian corals using the Bayesian framework of the program PyRate (v3.0)12,36,37. This program uses fossil occurrence data to calculate the temporal variation in rates of preservation, speciation and extinction, while incorporating multiple sources of uncertainty12. At its core implementation, PyRate jointly estimates the times of origination (Ts) and extinction (Te) for each fossil lineage; the fossilization and sampling parameters that determine preservation rates (q); and the overall rates of speciation (λ) and extinction (μ) through time36. Recently, the program has been upgraded to include a reversible jump Markov Chain Monte Carlo (rjMCMC) algorithm to estimate diversification rate heterogeneity, which provides more accurate and precise estimates than other commonly used methods12. Therefore, despite the inherent bias of the fossil record (i.e., estimates are conditioned on sampled lineages), PyRate is a robust method to quantify speciation and extinction rates, and their respective temporal shifts, from fossil occurrence data.Extant taxa can also be included in the PyRate framework as long as they are also represented in the fossil record. This is done to extend the fossil geologic ranges to the recent times. Hence, the first step in our analysis was to identify which species in our dataset is still alive at the present. To do this, we matched the accepted species names in the PBDB dataset with those from the extant species dataset of Huang et al.38. Subsequently, we split our dataset into eleven independent subsets, with the goal of keeping each subset with an equal number of species. Each data subset included a random selection of species with their respective occurrences, which was enough to calculate Ts and Te (see below). This was done to avoid convergence issues, given the large size of our dataset and the consequent complexity of the model37. For each of our subsets, we generated fifty replicates by resampling the fossil ages from their temporal ranges to account for the uncertainty associated with the age of occurrences. We then used the maximum-likelihood test in PyRate to compare between three models of fossil preservation12: the homogeneous Poisson process (HPP; q is constant through time); the nonhomogeneous Poisson process (NHPP; q varies throughout the lifespan of a species); and the time-variable Poisson process (TPP; q varies across geological epochs). The latter model (TPP) was selected across all of our data subsets (Supplementary Table 2).After selecting the preservation model, we first focused on assessing the estimates of times of origination and extinction in each data subset, rather than using the full dataset to jointly estimate all parameters at once as in the original implementation of PyRate. This further reduced the complexity of the model and allowed for more precise parameter estimates. For each replicate in all of our data subsets, we approximated the posterior distribution of Ts and Te through a 50 million generation run of the rjMCMC algorithm under the TPP, sampling parameters every 40 thousand iterations. At the end of each run, we discarded 20% of the samples as burn-in and assessed chain convergence through the effective sample sizes of posterior parameter estimates, using the software Tracer39 (v1.7.1).From the results of this first set of models, we extracted the median estimates of Ts and Te across replicates, and we merged the estimates from the eleven independent data subsets. This merged data frame contained estimated times of origination and extinction for all coral lineages within our fossil dataset. We then used this merged Ts and Te data frame as input for another rjMCMC chain to finally estimate overall λ and μ through time, by applying the option -d in PyRate. In this option, Ts and Te for all fossil lineages are given as fixed values and, therefore, are not estimated by the model. The chain for this model was run for 100 million generations, sampling parameters at every 40 thousand iterations. Once again, we excluded 20% of the initial samples as burn-in and checked model convergence using Tracer. Finally, we calculated net diversification rates through time by subtracting the post burn-in samples of μ from λ.To explore the taxonomic idiosyncrasies in the evolutionary rates of reef corals, we selected the most abundant families on present-day coral reefs in terms of the number of colonies per area18 (Acroporidae, Agariciidae, Merulinidae, Mussidae, Pocilloporidae, and Poritidae). Altogether, species within these families account for ~40% of the total extant diversity in Scleractinia. These families also account for most of the occurrences in the PBDB fossil dataset (excluding extinct families, which are generally older and had little temporal overlap with extant ones): Acroporidae (1457 occ. in 165 spp.); Agariciidae (722 occ. in 89 spp.); Merulinidae (2464 occ. in 229 spp.); Mussidae (1146 occ. in 100 spp.); Pocilloporidae (615 occ. in 64 spp.); and Poritidae (1149 occ. in 91 spp.). Therefore, from our full dataset, we selected six independent ones encompassing all species in each of the selected families. We also selected only species that are classified as reef-associated within these families, since we were specifically interested in these environments. This selection had a negligible effect on the size of the individual datasets, given that the vast majority of fossil species within these families are reef-associated. In each family, we followed the same modelling steps described above to estimate μ and λ, and diversity trajectories. However, this time it was not necessary to split the datasets into subsets, given that each family has far less occurrences than the full dataset. We started by comparing models of preservation, which showed the TPP as the best supported for all families (Supplementary Table 3). Then we created fifty replicates by resampling fossil ages to accommodate the uncertainty associated with the time of occurrences. For each replicate, we ran the rjMCMC algorithm for 50 million generations under the TPP model, with a sampling frequency of 40 thousand iterations. We discarded initial 20% of the samples as burn-in, and assessed convergence through Tracer. We then combined all replicates, resampling 100 random samples from each replicate to assess the estimates of μ and λ through time for each family. Finally, we extracted diversity trajectories in each family for all of the replicates by applying the -ltt option in PyRate, which generates a table with estimated range-through diversity at every 0.1 Myr. From these trajectories, we calculated the mean difference in diversity (slope in species per 0.1 Myr) between subsequent time samples backwards from the present (i.e., diversity in time t was subtracted from diversity in time t-1) using the diff function in R (v4.0.3).As an alternative to PyRate, we also calculated the diversity dynamics of reef coral fossils using the R package divDyn40, which combines a range of published methods for quantifying fossil diversification rates. Differently from PyRate, the metrics applied in divDyn require that the fossil occurrences are split into discrete time bins. Therefore, these metrics treat the origination and extinction rates as independent parameters in each bin, while PyRate is designed to detect rate heterogeneity through a continuous time setting12. Our goal here, however, was not to compare models but to assess the robustness of our rate patterns and diversity trajectories using alternative methods. We divided our dataset into one-million-year time bins to have enough temporal resolution for rate calculations. To account for the uncertainty in the assignment of fossil ages, we created 50 binned replicates by sampling the age of each occurrence from a random uniform distribution, with bounds defined by the age ranges provided in the PBDB dataset. We then used the divDyn function to calculate the per capita rates of origination and extinction through time (based on the rate equations by Foote41) for all scleractinians (Supplementary Fig. 5a) and for reef-associated acroporids alone (Supplementary Fig. 5b). We also used the same procedure to generate range-through diversity curves for each of the six families selected previously, to compare with the curves generated by PyRate (Fig. 2a). Although the rate results differed between the PyRate (Fig. 1) and the divDyn (Supplementary Fig. 5) approaches, the general patterns remained unchanged. Rates are more volatile through time in divDyn estimates, with larger confidence intervals, which is expected from the metrics applied in the package12,42. Yet, we found the same peaks in extinction for Scleractinia: at the Cretaceous-Paleogene and Eocene-Oligocene boundaries, and at the Pliocene-Pleistocene (Supplementary Fig. 5a). The recent peak in speciation in Acroporidae was also detected, although less strong (Supplementary Fig. 5b). Despite these slight differences in rate estimates, the diversity curves reconstructed through divDyn (Supplementary Fig. 6) mirrored almost exactly the ones found with PyRate (Fig. 2), demonstrating that the overall macroevolutionary trends described herein (Figs. 1 and 2) are robust to methodological choices.Diversity-dependent modelsTo assess the effects of diversity dependency on the evolution of reef coral lineages, we implemented the Multivariate Birth-Death model (MBD)11 within the PyRate framework. This method was first described as the Multiple Clade Diversity Dependence model (MCDD)19, in which rates of speciation and extinction are modelled as having linear correlations with the diversity trajectories of other clades. At its original implementation, the MCDD was developed to assess the effects of negative interactions, where increasing species diversity in one group can suppress speciation rates and/or promote extinction in itself or in other ecologically similar clades19. However, the model also incorporates the possibility of positive interactions, where increasing diversity in one clade can correlate with enhanced rates of speciation or buffered extinction. Through further model developments43, the MCDD was updated to also include a horseshoe prior44 on the diversity-dependence parameters, which helped controlling for overparameterization and enhanced the power of the model to recover true effects43. More recently, this model took its current form as the MBD11, with the additional possibilities of including environmental correlates and setting exponential, rather than just linear, correlations.We first applied the MBD to estimate the diversity-dependent effects of individual extant coral families (i.e., the ones selected in the previous analysis; see Evolutionary rates) in their combined diversity trajectories. From the rjMCMC model results for individual families, we extracted estimates of Ts and Te in each of the fifty replicates and merged them across families. This merged dataset with fifty replicates of Ts and Te was then used as input for the MBD model, where we set the relative diversity trajectories of each individual family as predictors. We also included three key environmental predictors—paleotemperature, sea level and rate of sea-level change—to assess their influence in overall evolutionary rates. The paleotemperature data was obtained from Westerhold et al.45, and consists of global mean temperature estimates for the last 66 million years, averaged across 0.1 Myr time bins. Eustatic sea-level data was downloaded from Miller et al.46, and contains estimates of sea level for the last 100 million years in comparison to present-day levels, also split in ~0.1 Myr time bins. With this dataset, we calculated the average rate of sea-level change per million years, as measured from the absolute difference between subsequent sea-level values backwards in time (i.e., sea level in time t was subtracted from sea level in time t-1). These environmental factors were rescaled between 0 and 1 to maintain all predictors on the same relative scale.Under our MBD model, the speciation and extinction rates of all families combined could change through time and through correlations with the relative diversity of individual families or environmental factors. The strength and directionality (positive or negative) of the correlations are also jointly estimated for each predictor within the model11. We ran both linear and exponential correlation models (see formulas in Lehtonen et al.11) in each of our fifty replicates for 25 million generations, sampling parameters at every 25 thousand iterations. We then compared the linear and exponential models through the posterior harmonic means of their log likelihoods, which supported the exponential one as having a better fit. From the posterior estimates, we summarized the speciation (Fig. 3c) and extinction (Fig. 3d) correlation parameters (i.e., the strength of the effect) by calculating their median and 95% Highest Posterior Density (HPD) interval across replicates. Finally, we also summarized the effect of families on lineage turnover (Fig. 3e), which we conceptualize as the sum of the effects on speciation and extinction.The MBD model also provides posterior samples of the weight of the correlation parameters, which is estimated through the horseshoe prior11. In essence, this prior is able to reliably distinguish correlation parameters that should be considered noise from those that represent a true signal in the data11. The parameterization of the horseshoe prior contains local and global Bayesian shrinkage parameters44 from which shrinkage weights (w) can be calculated (see formulas in Lehtonen et al.11). These shrinkage weights associated with each correlation parameter in the MBD model vary between 0 and 1, with values closer to 0 representing noise and values closer to 1 representing a true signal. Through simulations, it has been shown that values of w > 0.5 indicate that the correlation parameter in question significantly differs from the background noise, being the correlation positive or negative43. However, as a conservative way to infer the weight of correlation parameters, here we use a value of w > 0.7 to detect significance. This value was calculated for each diversity-dependence parameter (speciation, extinction and turnover) from the median values drawn from the model posteriors.The spatial distribution of reef-associated taxa varied considerably throughout the Cenozoic, with biodiversity hotspots moving halfway across the globe47. Therefore, the best way to capture this dynamic biogeographic history in reef corals is by analysing global diversity patterns like we did in our main MBD model. However, to assess the robustness of our diversity-dependent results against the influence of geographic scale and site co-occurrences, we repeated all the modelling steps described above with two data subsets. First, we selected only fossil species that have occurrences in the Indo-Pacific Ocean (i.e., 30°W–180°W) within the six families. Second, we excluded sites in which the Acroporidae did not co-occur with the other families. In each of these data subsets, we calculated diversity trajectories and used them as predictors in a separate MBD model. These models had a merged dataset of Ts and Te of all species included in each case (Indo-Pacific and co-occurrences) as a response variable.Finally, we followed the same modelling procedures described above to investigate the diversity-dependent effects in family pairwise analyses. We applied the MBD model to assess the effects of all other families in each individual family at a time, while also estimating correlations with the key environmental predictors. From the rjMCMC model results for individual families, we extracted the fifty replicates of estimated Ts and Te. Each replicate was then used as input for an MBD run using the relative diversity trajectories of each other individual family as predictors, along with the environmental variables. Once again, we ran 25 million generations of the MBD, with a sampling frequency of 25 thousand, using both linear and exponential correlation models in each age replicate. For all families, we found that the exponential model had a better fit. We then summarized the correlation parameters and the shrinkage weights (Supplementary Fig. 7) derived from the exponential models per family by calculating the median and 95% confidence intervals across replicates.Reporting summaryFurther information on research design is available in the Nature Research Reporting Summary linked to this article. More
