Drivers of beta diversity in modern and ancient reef-associated soft-bottom environments

Triassic data

The Cassian Formation in the Dolomites, Southern Alps, northern Italy, preserves Middle to Late Triassic (Ladinian–Carnian) tropical reef to basin environments with exceptional fossil preservation, allowing comparisons with recent assemblages due to a low taphonomic bias (Roden et al., 2020). The Cassian Formation comprises deposits with considerable differences in depth, from back-reef and lagoonal settings as well as shallow and deeper water deposits from the reef basin (e.g., Bosellini, 1998; Fürsich & Wendt, 1977). The platform-to-basin relief has been proposed to have exceeded 100 m, possibly even reaching bathyal depths (Keim et al., 2001; Urlichs, 2012). The predominantly argillaceous basin sediments were deposited between prograding carbonate platforms that now form the Cassian and Schlern Dolomite (Bosellini, 1998; Hausmann & Nützel, 2015; Keim et al., 2001) and can reach a thickness of over 300 m (Fürsich & Wendt, 1977). The Cassian Formation sensu lato includes all clay-rich Ladinian–Carnian sediments deposited in the interplatform basins of the Dolomites (Part SI). Deposition took place in the Western Tethys in a setting comparable to recent tropical environments, with warm water temperatures, seasonality, and fresh water influx (Nützel, Joachimski & Correa, 2010). Diagenetic alteration and lithification were low and fossil extraction is easy, yielding many well-preserved fossils with original skeletal microstructures and aragonite preservation (Roden et al., 2020). Many localities from this reef-basin assemblage have yielded abundant fossils, ideal for the assessment of diversity patterns and paleoecological studies (Fürsich & Wendt, 1977; Nützel, Joachimski & Correa, 2010; Roden et al., 2020). In a large quantitative study of the Cassian fauna, Fürsich & Wendt (1977) described autochthonous assemblages that are thought to represent communities. These contain recurring sets of species, which favors the idea of community-assembly processes and dynamics determining faunal heterogeneity (Drake, 1991).

Surface and bulk samples from the Cassian Formation were gathered in field campaigns conducted in 2015 and 2016 (Fig. 1). Studied samples belong to the aon, aonoides, and austriacum zones, the time span that covers the vast majority of diverse benthic assemblages of the Cassian Formation sensu lato (Fig. S1), equivalent to approximately 5 myr. The entire formation is distributed over an area of c. 500 km².

Bulk samples were disaggregated in a 7% H₂O₂ solution and wet-sieved with mesh sizes of 2, 1, 0.5, 0.25, and 0.125 mm. Only size fractions >0.5 mm from bulk samples are included in this study. Samples were picked and sorted using a light microscope. Roden et al. (2018) found that beta diversity, calculated as mean proportional dissimilarity, is depicted accurately when only the most abundant taxa in each sample are counted. Whereas the five most abundant species are usually sufficient, we identified and counted the ten most abundant species to avoid issues with samples of high evenness (Roden et al., 2018). All animal taxa were included, providing comparability with previously published data. In disarticulated bivalves, the few incomplete specimens were each counted as single specimens, since valves were of differing shapes and sizes. Only mollusks are among the top ten animal species in all four studied samples.

We inferred paleo-water depth following the criteria of Fürsich & Wendt (1977). Inference is based on the ratio of suspension and deposit feeders, the ratio of carnivores and grazers, the proportion of articulated bivalves, the abundance and diversity of gastropods, the encrustation of specimens, and the presence of coral, sponge, and echinoderm fragments (Fig. 2). Weighting of factors is detailed in Part SII. Age, locality, and diversity of Cassian samples are provided in Table 1.

Field collection was authorized by E. Kustatscher, Museum of Nature South Tyrol. All material collected is accessioned into the collections of the Museum of Nature South Tyrol (collection numbers NMS PZO12554 to NMS PZO12598).

Modern data

The Recent samples are from a shallow-water area in the Northern Bay of Safaga in the Red Sea, Egypt (Zuschin & Hohenegger, 1998), representing a coral-dominated, subtropical setting with warm water temperatures and seasonality, high salinity, and a highly structured bottom topography reaching down to more than 50 m water depth (Piller & Pervesler, 1989; Titschack et al., 2010). Water temperature and salinity are without any obvious depth gradient due to complete water mixing (Piller & Pervesler, 1989). Terrigenous input along with nutrients occurs mainly along the coast and is due to fluvial transport during flash floods, local erosion of impure carbonate rocks, and aeolian transport by the prevailing northerly winds (Piller & Mansour, 1994). Water energy is relatively weak, but a complex current pattern influences facies development (Piller & Pervesler, 1989). The samples were collected and processed within the scope of previous studies (Zuschin & Hohenegger, 1998; Zuschin & Oliver, 2005).

Standardized bulk sampling in soft substrates was conducted at 13 sites, from shallow subtidal down to 40 m water depth (Table 2), covering an area of approximately 75 km². The dataset consists only of mollusks. Whole shells >one mm were considered and disarticulated valves were counted as individuals (Zuschin & Hohenegger, 1998; Zuschin & Oliver, 2005). For comparability, the Safaga dataset was also reduced to the ten most abundant species per sample. To test whether results are robust, the complete dataset was analyzed and results are provided in Part SIII. Samples taken from the same environment, site, and depth (only several meters apart) were combined to create a by-site dataset comparable to that from the Cassian Formation. Results from the by-sample dataset are provided in Part SIV.

Diversity estimates

We assess beta diversity patterns in the Triassic Cassian Formation based on 8 bulk samples and compare the results with those from the Red Sea samples. For both datasets, we use pairwise proportional dissimilarity (relative Bray–Curtis) to compare samples. Abundance-based indices, such as proportional dissimilarity, have been shown to be relatively insensitive to uneven sample sizes, as opposed to incidence-based indices, which are biased when species richness and/or sampling completeness vary (Krebs, 1999; Magurran, 2004; Wolda, 1981). Proportional dissimilarity is calculated as d_jk = 1 –∑min(x_ij, x_ik), with x_ij and x _ik being the proportions of species abundance in each sample. Ordination plots using non-metric multidimensional scaling are also based on pairwise proportional dissimilarity. The dispersion of homogeneity is based on pairwise proportional dissimilarity as well as the modified Gower measure (Anderson, 2006; Anderson, Ellingsen & McArdle, 2006). For both beta diversity indices, the centroid for the entire dataset is figured for both the Cassian and the Safaga dataset; in addition, the sites were grouped into shallow and deeper water groups (Part SV). Community composition data from samples within a depth group are combined and a Kruskal–Wallis rank sum test is performed to compare the rank distribution of the different depth groups. Non-metric multidimensional scaling is used to visualize differences in community composition. Permutational multivariate analysis of variance (PERMANOVA) using distance matrices is applied to determine significance of depth and age groups. Calculations and visualization were implemented using R Version 3.5.0 (R Core Team, 2016) and the vegan (Oksanen et al., 2016), sads (Prado, Miranda & Chalom, 2014), and visreg (Breheny & Burchett, 2013) packages.

A null model for each dataset is created by first pooling all taxa with their observed abundances, yielding a single vector of total abundances in the dataset. From this, we derived a vector of proportional abundances. The species pool is then randomly resampled, using the proportional abundance vector as probability weights, until the number of specimens and number of sampling sites of the original datasets are obtained (e.g., until there are 8 sampling sites each containing the number of specimens as the original samples). The mean proportional dissimilarity is calculated for each simulated dataset across 1,000 iterations.

Alpha-level (local) community diversity is deduced from rank-abundance distributions of the 10 most abundant species per sample as well as the Berger-Parker dominance index (the proportion of the most abundant species) and Pielou’s evenness (J = H/log(s), where H is the Shannon index and s is the number of species) (Berger & Parker, 1970; Krebs, 1999; Pielou, 1966). Hypothesis testing was done using Spearman’s rank correlation, as pairwise dissimilarity as well as alpha diversity values are not normally distributed. Distance decay, the decrease in biological similarity with spatial distance, is measured as the correlation between pairwise proportional dissimilarity and spatial distance.

Environments and alpha diversity

Quantitative faunal data of samples from the localities Costalaresc, Picolbach, Rumerlo cliff and Rumerlo ski slope are first reported herein (Table 3 and Figs. S14–S16). This dataset was supplemented with data from Settsass (Nützel & Kaim, 2014), Stuores (Hausmann & Nützel, 2015), and the two localities Lago Antorno and Misurina landslide (Hausmann et al., 2019, unpublished data), covering various environments (Table 1). The Stuores sample is one of the most diverse assemblages known from the Mesozoic (Hausmann & Nützel, 2015). The assemblage from Settsass is of moderate diversity and differs greatly in taxonomic composition from previously studied Cassian samples (Nützel & Kaim, 2014). The samples from Lago Antorno and Misurina landslide are relatively similar in faunal composition and are both moderately diverse (Hausmann et al., 2019, unpublished data). All samples stem from soft-bottom habitats.

The samples are arranged along a shallow reefal to basinal transect (Fig. 3) based on inferred water depth (Table S1). The bathymetric partitioning among the Cassian localities is only relative, rendering a comparison with recent settings difficult. There were no apparent differences in the sediment matrix among samples. Due to the high proportion of grazers and high gastropod abundance and diversity, we describe Costalaresc as originating from a relatively shallow environment, but—with a high proportion of deposit feeders and all specimens of bivalves being articulated—we interpret a slightly deeper setting than for the Rumerlo localities. Rumerlo ski slope is interpreted as a very proximal back reef setting, due to the large number of fragments of corals, sponges, and echinoids, among other factors (see Part SII). Picolbach probably stems from a deeper setting, as interpreted from mode of life of reported specimens and articulated bivalves. Evenness values, Berger-Parker dominance index (Table 1), and rank-abundance distributions (Fig. S12, Part SVI) demonstrate relatively diverse assemblages, with the exception of Costalaresc. Low alpha diversity in Costalaresc is due to the dominance of the gastropod Helenostylina convexa. There is no significant correlation between inferred water depth and alpha diversity measured as dominance (Spearman’s rho = 0.02, p = 0.98) or evenness (rho = 0.29, p = 0.50) in the Cassian samples.

Environments and depth from which the Safaga samples were taken are recorded and alpha diversity calculated (Table 2). Samples from the reef slope and from sand between coral patches show lowest dominance and highest evenness. There is no significant correlation between alpha diversity and depth in the Safaga samples (dominance: rho = −0.05, p = 0.93, evenness: rho = −0.12, p = 0.78). The range of evenness is slightly higher (J = 0.39 to 0.55) than in the Cassian samples (J = 0.29 to 0.58).

Beta diversity

As shown previously, overall beta diversity is high in the Cassian Formation (mean PPD: 0.91 ± 0.02, ± is standard error; range: 0.52–1; Table 4). In the by-site dataset from Safaga, we measure a beta diversity of 0.89 ±  0.04 (range: 0.35–1.00; Table 5). Null models created for each dataset from the gamma species pool yield much lower beta diversity, with a mean of 0.24 ± 0.0004 for the Cassian dataset (Fig. 4A) and 0.10 ± 0.0002 for the by-site Safaga dataset (Fig. 4B). Histograms show the distribution of pairwise dissimilarity values (Fig. S13). There is no significant correlation between spatial distance and dissimilarity neither in the Cassian samples (Spearman’s rho = 0.14, p = 0.48; Fig. 5A) nor the by-site Safaga dataset (rho = 0.05, p = 0.81; Fig. 5B).

Dissimilarity between depth categories does not differ much from dissimilarity within depth groups in either dataset. Grouping Cassian samples by depth yields a dissimilarity of 0.83 between the four shallower and the four deeper localities (see Table 1 and Table S1). Using three depth categories (the four shallowest localities, two intermediate and two deeper localities, see also Fig. 3), we find the following dissimilarities: shallow/intermediate: PPD = 0.83, intermediate/deep: PPD = 0.78, shallow/deep: PPD = 0.94). Dissimilarity between the two localities from deeper environments is lower than among the shallower localities, while the two sites from an intermediate depth are very dissimilar (Fig. 3). Non-metric multidimensional scaling shows no clear distinction between localities based on depth (Fig. 6). Neither two nor three depth groups are statistically supported by PERMANOVA (R² = 0.21, p = 0.15; R² = 0.23, p = 0.06, respectively). The Cassian data show relatively broadly dispersed sites when plotting the dispersion of homogeneity of all assemblages. Differentiating between shallow- and deeper-water assemblages yields two separate groups for PPD but overlapping centroids when the modified Gower measure is applied (Fig. S10 , Part SV). Kruskal–Wallis rank sum tests show that the shallow- and deeper-water assemblages have statistically indistinguishable distributions (p = 0.07). Using three groups, rank sum tests also yield the same rank distributions (shallow/deep: p = 0.07, shallow/intermediate: p = 0.31, shallow/deep: p = 0.06). There is no correlation of paleo-depth with mean PPD of each site (Spearman’s rho = 0.24, p = 0.58). Temporal turnover between the ammonite biozones is only slightly higher than dissimilarity within the ammonite biozones, which strongly varies (Part SI). We find age to explain 43% of the variation in community dissimilarity (p = 0.02).

At Safaga, there is also no clear relationship between beta diversity and water depth (Fig. 7). As in Cassian, there is no distinct association of samples from similar water depths at Safaga (Fig. 8). Grouping samples into two or three depth groups yields non-significant clusters (R² = 0.17, p = 0.28). Dispersion of homogeneity yields similar results as in the Cassian dataset (Fig. S11). Samples taken from deep, muddy settings exhibit the highest mean beta diversity relative to other samples (0.98 ±  0.02). Otherwise, there is no relationship between sedimentary attributes and mean dissimilarity. By grouping the samples into depth ranges, we cover several environments for each. Samples from shallower environments have a more similar community composition than samples from deeper environments (Figs. 7 and 8). There is no correlation of water depth with mean PPD of each site (Spearman’s rho = 0.02, p = 0.98). Dissimilarity between the four shallower and the four deeper samples is 0.77 and therefore lower than other values measured within depth ranges. Using three depth categories (the four shallowest localities, two intermediate and two deeper localities), we find the following dissimilarities: shallow/intermediate: PPD = 0.79, intermediate/deep: PPD = 0.60, shallow/deep: PPD = 0.77). In contrast to the Cassian data, Kruskal–Wallis rank sum tests show that all depth groups have different rank distributions (p <  10⁻¹⁶) in the Safaga dataset.

Environment

We find no distinct trends in alpha diversity with depth in either dataset. Measures of evenness and dominance show relatively high alpha diversity for most samples in both datasets, from all environments and depths. Low dominance is corroborated by another study on soft-bottom taxa in the central Red Sea (Alsaffar et al., 2019). Patterns of a bathymetric diversity gradient in benthic fauna in the literature are conflicting (e.g., Brown & Thatje, 2014; Gray et al., 1997; Holte, Oug & Cochrane, 2004; Martins, Quintino & Rodrigues, 2013).

Beta diversity has not been studied extensively in soft-bottom habitats, with only a few studies in non-reefal soft-bottom settings (e.g., Aldea, Olabarria & Troncoso, 2009; Ellingsen, 2002; Koulouri et al., 2006). Our study shows high beta diversity in two reef-associated soft-bottom assemblages from warm-water settings. We find no significant contributions of water depth, spatial distance, and stochastic processes to beta diversity. There is a moderate influence of age in the Cassian samples.

Time averaging generally lowers beta diversity in an assemblage (Tomašových & Kidwell, 2009). However, we can assume that the individual assemblages in the Cassian are time-averaged to a similar degree as in Safaga (i.e., within the same habitat over hundreds to few thousands of years). But although time averaging does not hamper the comparison, the differing temporal ranges may (5 myr vs. <1,000 years, respectively). High beta diversity within ammonite biozones in the Cassian points to a true ecological signal, but unfortunately, we do not know what temporal range lies between two samples from the same biozone. Temporal differences between some samples are in the range of millions of years.

The high beta diversity in the Cassian Formation is remarkable, considering that all samples stem from soft-bottom habitats adjacent to reefs—most likely a more uniform habitat than the reef structure. The slightly lower beta diversity of modern Safaga compared with ancient Cassian is surprising considering that certain environments such as mangrove and seagrass, which contribute to beta diversity at Safaga, did not yet exist in the Triassic period (Hemminga & Duarte, 2000; Ricklefs, Schwarzbach & Renner, 2006). In addition, the bottom topography of Safaga is highly structured for a soft-bottom habitat (Zuschin & Oliver, 2005).

Results from grouping the Cassian localities by inferred depth are ambiguous. While dissimilarity between the two localities from deeper environments is lower than most other values, data are too limited to allow for robust conclusions on beta diversity vs. depth. However, the Safaga samples—with a larger number of sampling sites and specific recorded depths—yield slightly higher beta diversity in deeper habitats. In addition, dissimilarity between two depth categories (four shallower samples vs. four deeper samples) is lower than within depth categories for both datasets. The high beta diversity at both Cassian and Safaga may partly be driven by the great range of environments included. Specifically, beta diversity has been shown to increase with the variance of depths, due to the larger range of communities included in the study (Harborne et al., 2006).

Change in faunal composition with depth is usually linked to variations in environmental factors, such as temperature, salinity, oxygenation, or sediment characteristics (Durden et al., 2015; Ellingsen & Gray, 2002; Holte, Oug & Cochrane, 2004; Laine, 2003). Due to water mixing in the Bay of Safaga (Zuschin & Oliver, 2005), there is probably no depth gradient related to temperature, salinity, or oxygenation, but light penetration may play a role. Differences in sedimentary attributes are not related to faunal dissimilarity. While samples from a muddy environment have the highest mean dissimilarity with the other samples from Safaga (see also ordination in Zuschin & Hohenegger, 1998: fig. 11a), this may be related to other factors, as there is no trend in beta diversity from sand to muddy sand to mud. Muddy habitat is common in the Red Sea in deeper waters and occurs in a shallower settings in protected depressions. Especially the very high dissimilarity between samples from similar depths at Safaga and in the Cassian contradicts findings by Ellingsen & Gray (2002) that samples from similar depths show higher similarity. While the Safaga samples differ in habitat and grain size, differences in sedimentary attributes were not noted in the Cassian Formation, leading us to assume that other factors may play a larger role in determining faunal composition. While depth is often found to be an important driver of beta diversity (Ellingsen & Gray, 2002; Harborne et al., 2006), we find no changes in beta diversity with depth in our sites. However, depth can determine faunal composition through depth-range restriction of species (Rex & Etter, 2010), which may explain high overall beta diversity when differences in water depth are large.

Comparing dissimilarity among different datasets is made more difficult by the fact that spatial resolution influences beta diversity. Small spatial scales yield higher beta diversity than larger scales in coral reefs and rain forests, as the heterogeneity of environments is limited (Pandolfi, 2002). However, studies on the effects of spatial scale on beta diversity reach differing conclusions (Lennon et al., 2001). Specifically, small-scale environmental differences—often created by the biota itself—can increase faunal variability. These small-scale differences can in turn increase beta diversity at a larger scale (e.g., Thrush, Lundquist & Hewitt, 2005). Applying Pandolfi’s (2002) three-phase model of variability to the reef-adjacent soft bottoms of the Cassian Formation and the Bay of Safaga, we find that variance in community composition is expected to be lower at the observed spatial extent (1–100 km distance) of both datasets than at smaller (<1 km) or larger (>1,000 km) scales. We therefore conclude that the high values of beta diversity are genuine and not due to spatial scale.

Stochasticity and distance decay

The lack of correlation between spatial distance and community dissimilarity in the two sites lets us assume that there is no noteworthy distance decay over the limited spatial extent covered by the fossil and modern datasets. Studying an environmentally very homogeneous soft-bottom habitat on the Norwegian continental shelf, Ellingsen (2002) also found a weak relationship between spatial distance and dissimilarity. In addition, beta diversity in both study sites is much higher than predicted from the null models. Therefore, the spatial patterns of community composition among the study sites are not random and probably not constrained by dispersal. Since the two datasets differ in size (the Safaga dataset contains a higher number of samples, specimens, and species over a smaller area), the higher beta diversity in the Cassian null model is probably due to increased randomness by sampling fewer specimens and species from a slightly smaller species pool. With a beta diversity of 0.91 for the Cassian Formation despite the small species pool, stochasticity is clearly not the main driver of beta diversity. With the lack of stochasticity on top of the lack of a spatial pattern related to distance, other factors must have a stronger control on beta diversity.

Community assembly

Although we find no significant change in beta diversity with water depth, the factors driving beta diversity may vary with depth. In shallower settings, physical disturbances, such as storms, may have a large impact on benthic communities. In greater depths, competition may be a stronger driving force due to less prevalent disturbances and predation (Harper & Peck, 2016; Klompmaker & Finnegan, 2018). Klompmaker & Finnegan (2018) hypothesize predation to be an important factor in structuring modern and fossil soft-bottom communities, supporting earlier assessments (Stanley, 2008). However, predation has increased since the Triassic (Huntley & Kowalewski, 2007), therefore the Cassian fauna might have been less affected by predation than the Safaga fauna.

Environmental heterogeneity is often considered the main driver of faunal heterogeneity (e.g., Carlos-Júnior et al., 2019; Ellingsen, 2002; Pandolfi & Jackson, 2001; Stegen et al., 2013), with more uniform environments generally yielding a more homogeneous fauna (e.g., Clarke & Lidgard, 2000; Ellingsen & Gray, 2002; Fagerstrom, 1983; Shin & Ellingsen, 2004). However, it is often ignored that the environment can be structured by the biota itself. Organisms with hard parts that live on the sediment are found to directly contribute to environmental heterogeneity and in turn to faunal heterogeneity (Hewitt et al., 2005). The studied assemblages include both epifauna and infauna, which may directly contribute to high beta diversity through an increase in diversity driven by the combination of higher sediment stability in infaunal assemblages and biotic interactions in epifaunal assemblages (Van der Zee et al., 2015).

Besides structuring of the seafloor by the biota, we presume priority effects may be one of the drivers of beta diversity in the modern and ancient reef-related soft-bottom habitats. While there are many mechanisms that lead to differences in community composition, the order in which organisms settle in a habitat directly affects the community that subsequently assembles (Drake, 1991; Fukami, 2015; Thrush, Lundquist & Hewitt, 2005). Priority effects influence community composition at smaller sites by affecting the regional species pool as well as local population dynamics (Fukami, 2015). In addition, the first species to arrive may even gain an evolutionary advantage, as they adapt to environmental conditions sooner than subsequent arrivals (De Meester et al., 2016). Further research may reveal one or both of these effects as drivers of beta diversity. Our results show that water depth, stratigraphic age, geographic distance, and random dispersal are not the key determinants of beta diversity in the studied reef basin assemblages.