FACILITATING CORALS IN AN EARLY SILURIAN DEEP-WATER ASSEMBLAGE

: Corals are powerful ecosystem engineers and can form reef communities with extraordinary biodiversity through time. Understanding the processes underlying the spatial distribution of corals allows us to identify the key biological and physical processes that structure coral communities and how these processes and interactions have evolved. However, few spatial ecology studies have been conducted on coral assemblages in the fossil record. Here we use spatial point process analysis (SPPA) to investigate the ecological interactions of an in situ tabulate and rugose coral community (n = 199), preserved under volcanic ash in the Silurian of Ire-land. SPPA is able to identify many different sorts of interactions including dispersal limitation and competition within and between taxa. Our SPPA found that the spatial distribution of rugose corals were best modelled by Thomas clusters ( p d = 0.834), indicating a single dispersal episode and that the tabulate corals were best modelled by double Thomas clusters ( p d = 0.820), indicating two dispersal episodes. Further, the bivariate distribution was best modelled by linked double clusters ( p d = 0.970), giving signiﬁcant evidence of facilitation between the tabulate and rugose populations, and identifying the facilitators in this community to be the tabulate corals. This interaction could be an important ecological driver for enabling the aggregation of sessile organisms over long tem-poral periods and facilitation may help to explain trends in reef diversity and abundance during the Ordovician biodiversiﬁcation and in the early Silurian.

Abstract: Corals are powerful ecosystem engineers and can form reef communities with extraordinary biodiversity through time. Understanding the processes underlying the spatial distribution of corals allows us to identify the key biological and physical processes that structure coral communities and how these processes and interactions have evolved. However, few spatial ecology studies have been conducted on coral assemblages in the fossil record. Here we use spatial point process analysis (SPPA) to investigate the ecological interactions of an in situ tabulate and rugose coral community (n = 199), preserved under volcanic ash in the Silurian of Ireland. SPPA is able to identify many different sorts of interactions including dispersal limitation and competition within and between taxa. Our SPPA found that the spatial distribution of rugose corals were best modelled by Thomas clusters (p d = 0.834), indicating a single dispersal episode and that the tabulate corals were best modelled by double Thomas clusters (p d = 0.820), indicating two dispersal episodes. Further, the bivariate distribution was best modelled by linked double clusters (p d = 0.970), giving significant evidence of facilitation between the tabulate and rugose populations, and identifying the facilitators in this community to be the tabulate corals. This interaction could be an important ecological driver for enabling the aggregation of sessile organisms over long temporal periods and facilitation may help to explain trends in reef diversity and abundance during the Ordovician biodiversification and in the early Silurian. P U T A T I V E corals have been known since the early Cambrian (Sorauf & Savarese 1995;Hicks 2006) but it was not until the Middle Ordovician when two major orders, Tabulata and Rugosa, radiated to ecological prominence (Webby et al. 2004;Lee & Riding 2018). These remained the two major orders of Palaeozoic corals until their demise in the Permo-Triassic mass extinction (Fedorowski 1989), with Silurian corals established as important reef-builders (Webby 2002) as exemplified by Wenlock reef frameworks comprising of branching tabulate and rugose corals as well as stromatoporoid sponges and bryozoans (Scoffin 1971). Notably, reefs are important in generating and exporting diversity to different environments over geological time (Kiessling et al. 2010). Therefore, understanding the ecology of Palaeozoic corals is crucial for understanding Palaeozoic macro-ecology. A crucial aspect of reef ecology is understanding the mechanisms by which sessile reef-building organisms interact with one another and how these interactions influence and structure the ecosystem as a whole (Wood 1999;Karlson et al. 2007).
The spatial distributions of sessile organisms reflect the biological and physical interactions within and between their populations and with their environment (Illian et al. 2008;Wiegand & Moloney 2013). Therefore, analyses of the spatial distribution of in situ coral communities in the rock record has the potential to illuminate the ecological and environmental influences that influenced and structured fossilized communities through deep time. Spatial point pattern analysis (SPPA) is a set of powerful techniques to enable the analysis and modelling of the ecology of sessile organisms (here fossil specimens), by describing fossils as spatial points in an observational window (the exposed bedding surface) (Illian et al. 2008;Wiegand & Moloney 2013). Different properties can be allocated to spatial coordinates such as size, orientation, and/or taxonomic identification. The relative spatial patterns of these properties, known as 'marks', can be analysed (Illian et al. 2008;Wiegand & Moloney 2013). The key advantage of quantifying the spatial patterns of organisms is that it allows the statistical comparison of observed patterns with known ecological models. This enables the inference of the most likely underlying process behind an observed pattern (Illian et al. 2008). SPPA has been developed and used widely in extant plant species distributions to resolve influences such as physical environment (Wiegand et al. 2007a), organism dispersal/ reproduction (Seidler & Plotkin 2006) and competition for resources (Getzin et al. 2006). Application to the fossil record is more limited, but SPPA has been applied to Ediacaran organisms to infer reproductive traits (Mitchell et al. 2015), competitive interactions (Mitchell & Kenchington 2018) and interactions with different environmental settings (Mitchell et al. 2019(Mitchell et al. , 2020. Recently, SPPA has also been applied to investigate drilling predation distribution in bivalves (Rojas et al. 2020).
The application of SPPA to corals, both extant and fossil, has been limited. Extant coral-related SPPA has been used to evaluate the spread of diseases in coral-reef communities (Zvuloni et al. 2009;Easson et al. 2013;Deignan & Pawlik 2015) and more recently to infer the mortality, population and community dynamics of a deep-sea coral and sponge community (Mitchell et al. 2020). SPPA has also been suggested as a method of monitoring changes in coastal benthos over time (Piazza et al. 2020) as well as investigating mortality due to adult proximity (Gibbs & Hay 2015) or density-dependence (Edmunds et al. 2018). However, to our knowledge SPPA has not been used to analyse coral palaeocommunities.
The majority of analyses of non-encrusting fossil coral spatial patterns have used the nearest-neighbour metric (e.g. Rich 1981;Weidlich et al. 1993). However, nearest neighbour analyses are limited by the maximum distance between two points: if the maximum distance between neighbouring points is 50 cm, spatial scales over 50 cm go undescribed by the analysis, which does not allow for resolution of more complex ecological patterns (Mitchell & Butterfield 2018). Complex patterns over larger spatial scales can be analysed using Pair Correlation Functions (PCFs) which describe the density variations of points as a function of radius (Fig. 1). The value of PCF is the mean density of a set of points (in this case the centres of corals) within a given radius r. PCF = 1 describes complete spatial randomness, PCF > 1 describes aggregation (or clustering) and PCF < 1 describes segregation for a given value of r. Hence the higher the value of PCF, the greater the aggregation at that radius, and vice versa. For example, when PCF = 4, the density is four times greater than the mean density, and when PCF = 0.5, the density half that of the mean density. To avoid overfitting to noise, the observed PCF distribution is smoothed by averaging points over a distance (dr) dependent on the sample size; the more points, the less smoothing is required, and for a given distribution the optimal smoothing is calculated (Illian et al. 2008). This smoothing is required because a small number of points may easily create outliers (noise), which can change the result if not taken into account. The calculated PCFs can then be compared to different models using Monte Carlo simulations and goodness of fit tests to infer spatial processes that may have resulted in the given distribution (Wiegand et al. 2007b(Wiegand et al. , 2009Diggle 2013;Wiegand & Moloney 2013).
For sessile organisms there are four different sets of processes that that could underlie a given spatial distribution: (1) dispersal; (2) habitat association; (3) intra or inter-specific competition and/or facilitation; and (4) density-dependent mortality (Seidler & Plotkin 2006;Wiegand et al. 2007b;Getzin et al. 2008;Lingua et al. 2008). For instance, one can imagine in a community of undisturbed sessile organisms that the dominant control on the distribution is that of the dispersal of their propagules. There are two major types of dispersal method that extant corals employ: most Scleractinia are broadcast spawners (external fertilization), a subset are brooders (internal fertilization) where eggs are fertilized in the gastrovascular cavity and then released into the water column (Harrison 2011). Dispersal can be modelled in spatial ecology as cluster processes, where the point density within the cluster is normally distributed about the centre, known as Thomas clusters (TC) (Wiegand et al. 2007b). If there are two dispersal episodes of the same taxon (univariate) this leads to two clusters at different spatial scales. The first dispersal event forms a cluster around the initial colonizers, and the second set of clusters around the first generation; these dispersal events are modelled by double Thomas cluster (DTC) processes (Lin et al. 2011). However, if there is significant environmental patchiness (and therefore disturbance) then this patchiness may be more important than the dispersal in determining the spatial distribution of these organisms and can be modelled on a heterogeneous background with a Poisson distribution (HETP) (Lin et al. 2011). If the pattern of sessile organisms is consistent with complete spatial randomness (CSR), then this data would best be described by a homogenous Poisson distribution (Illian et al. 2008) as the Poisson distribution applies to discrete objects (points) in a space (observational window) that are not dependent on one another. When analysing the relationship between two taxa, their spatial distributions can be described by bivariate PCFs. Bivariate PCFs are calculated using a similar method to univariate (one taxon) PCFs; however, the PCF is calculated using the coordinates of one group in relation to the coordinates of the other (e.g. PCF < 1 implies the groups are segregated in relation to each other at the given value of r).

GEOLOGICAL SETTING
The Kilbride Peninsula (Republic of Ireland) straddles the Galway-Mayo border and houses Llandovery and Wenlock Epoch (Silurian) marine sediments of the Irish Midland Valley platform (Clarkson & Harper 2016). The succession shows a transgressive sequence from crossbedded sandstones of the Lough Mask Formation, interpreted as braided river facies that overlie the sandstones and siltstones of the marine Kilbride Formation, and succeeded by red shales of the Tonalee Formation, which were deposited in an off-shore setting (Fig. 2). The Kilbride Formation is especially fossiliferous, containing corals, brachiopods, trilobites and Skolithos burrows (Gardiner & Reynolds 1912).
One particularly well-preserved bedding plane in the Kilbride Formation shows a diverse fauna dominated by tabulate corals preserved in situ by a volcaniclastic surge (Harper et al. 1995). The volcaniclastic material can be traced along bedding surfaces that have a pale, feldspathic dust ( Fig. 3B, D). Weathering through the exposed surface has revealed coral colonies that have rapidly decalcified to leave moulds of the corals. Most of the corals can be identified to at least to genus level (Harper et al. 1995). Based on calcification growth increments, the most long lived organisms in the community survived for~15-20 years before their demise (Harper et al. 1995).
Sedimentation restricted the growth of the corals in some cases (Harper et al. 1995, fig. 4). Brachiopods of the Clorinda community in the upper part of the formation and an absence of storm-generated sedimentary features suggest a deep-shelf dysphotic environment (>50-60 m), with intermittent sedimentation (Brett 1995;Harper et al. 1995). This assemblage is one of the deepest known Silurian coral fauna (Scrutton 1998). This near-instantaneous, Pompeii-esque mass mortality allows us to investigate the community ecology of this deep-water assemblage, moments before their demise. We use SPPA to determine the controls on the spatial distribution of these extinct corals and its ecological implications.

Data collection
The fossiliferous surface (154°/~45°SE; 53°34 0 29.0″N, 9°26 0 48.9″W) was LiDAR scanned to a resolution of AE 1.68 mm with a FARO Focus X330 and photographed using a Canon 70D over the course of two days in August 2019. A 3D photogrammetric model of the surface was created in AgiSoft Photoscan (v1.4.5) and combined with a high-resolution LiDAR scan to ensure morphometric accuracy of the 3D model. A 2D projection of the model F I G . 1 . Simplified diagram for calculating PCFs on a hypothetical surface. A, the observation window and spatial data; coloured bands represent the radial bins each defined by dr. B, the PCF summary statistic calculated for one point using the radial bins defined in A: the total number of points in each radial bin is plotted against the radial distance from the reference point. In this example, a PCF is only computed for one point (within the red circle); in practice this computation is done for every point in the window and so PCF = 1 represents the mean radial density. A PCF line, joining the crosses and representing all observations, is omitted for clarity.
(perpendicular to bedding plane) was used to identify the corals on the surface. The corals were marked up as ellipses using Inkscape 0.92.3; the two axes of the ellipse were calculated, and their intersection of the axes defined as the coordinate that defines the position of any individual coral on the surface (see Mitchell et al. 2019). Approximating the coral position and dimensions as ellipses is warranted as the corals are in situ, aligned and their (usually round) top surfaces exposed (Fig. 3). The resultant vector map of 34.4 m 2 was used for subsequent analyses (Dhungana & Mitchell 2021, fig. S1). Corals could be identified to genus-level with the same taxa identified by Harper et al. (1995), and correspondingly had a high genus-level coral diversity (n > 20). Species level identification was not made as sometimes subtle species differences reduce the confidence of certain species level identification, and therefore impact the accuracy of the analyses. Identification of specimens below class was not always possible due to erosion of taxonomically relevant features. We chose to analyse at the level of order (Tabulata or Rugosa) as there were two groups with n > 30 (Table 1) and finer taxonomic analyses (e.g. genera) would have had smaller sample sizes, therefore making it harder to detect signals within the community.

Spatial analyses
Four different spatial distributions were investigated in this study: (1) all the specimens within the community (ALL); (2) tabulate corals (TAB); (3) rugose corals (RUG); and (4) rugose and tabulate corals (TAB-RUG). Three of these are univariant distributions (ALL, TAB, RUG) and one is bivariant (TAB-RUG). Preliminary analyses were performed in R (using the package Spatstat F I G . 2 . Simplified geological map and stratigraphic column of part of the Kilbride peninsula) Scale bar on main map represents 1 km. Smaller inset map is of the Republic of Ireland, with a pin denoting the position of the Kilbride Peninsula. Silurian strata are coloured for emphasis. Locality surveyed in this study is denoted on the map and the stratigraphic column with a star. Note the Bencorrag Fm. is included in the stratigraphic column for completeness and does not crop out in this section of the peninsula. Dating information from Leake (2014) and references therein. *Exact dating of the Dalradian schist is not known, therefore the geological Era (Neoproterozoic) is given in the Period column.
To determine the structure underlying the observed spatial distributions, two different methods were used to compare the fossil spatial distributions to different ecological models. The spatial distribution model fit was evaluated using Diggle's goodness-of-fit test (Diggle 2013) and visual comparison to a 999 Monte Carlo simulation envelope with the lowest and highest 49 simulation values removed (see Mitchell et al. 2019). Each process is simulated for the same number of observed points (fossils) for each observational window shape (sampled area) and size. Note, this envelope is not a confidence interval due to the non-independence of points and because the size of the simulation envelope will decrease as sample size increases. Diggle's goodness-of-fit test calculates a p d value which represents the total squared deviation between the observed pair correlation function and the model pair correlation function (Diggle 2013). For the goodness-offit test, p d = 1 describes a perfect fit of the calculated  . For example, if there are significant aggregations or segregations at smaller spatial scales, and CSR for larger spatial scales, models should only be fit to the non-CSR portions of the PCF. To simulate a heterogeneous background to the data (e.g. for HETP), an estimate of the background density distribution is created from the observed data. These maps are defined by a mean of the density of points with a circle of radius R, smoothed over the whole surface. Here we created multiple HETP models in which we changed r from 1 cm to 1 m in 10 cm increments to capture different scales of heterogeneity; the best fit model radius was used. The best fit density radius for all distributions on the surface was 50 cm (Table 1). To aid computation, instead of determine the densities over a continuous scale, the radial distances are discretized into a series of radial bins of ring width (dr).
To find the best-fit models to the data the following procedure was followed: 1. The PCF and L-function were calculated for each distribution (Levin 1995). The L-function is a linearized form (varying with r) of Ripley's K-function which measures spatial association as a cumulative function of radius (and therefore varies with r 2 for homogenously distributed data). 2. For each of the distributions, the regions where PCF > 1 were fitted to Thomas-cluster processes for the PCF. Not fitting to the random fluctuations around PCF = 1 ensures that actual aggregations are being modelled. Programita used the minimal contrast method to find the best-fit model ( 7. Finally, simulations were conducted on the univariate models CSR, HETP, TC, DTC and the bivariate models CSR, SPTC and LDTC. The coral surface areas (from the marked spatial data) of the two univariate resolved groups (RUG and TAB) were analysed using Bayesian information criterion (BIC) to resolve the populations in each distribution performed using the R package mclust v. 5.4.7 (Fraley & Raftery 2006). Gaussian mixed models with different numbers of clusters (in this study populations of corals) and parameters (means and standard deviations) can be compared by using the value of the maximized log likelihood function (which describes the probability of observing the values by the given model) with a penalty on the numbers of different parameters in the model is represented by the BIC value (Schwarz 1978;Fraley & Raftery 2007). Within the total distribution, there are two sets of cohorts that are fitted to the data. First, the best-fit cohorts that have Gaussian distributions with different means (and samples) but the same variance are fitted, then second where the variance is allowed to be different. The number of components are the number of these Gaussian distributions that fit within the total observed distribution. Two models with a difference of BIC scores larger than 10 are strongly significantly different; between 5 and 10 is significant, 3 to 5 weakly significant and less than 3 not significant (Fraley & Raftery 2007).

RESULTS
This community was dominated by the tabulate corals, comprising of 71% of the identified corals (Table 1) with the rugose corals consisting of 29% of the mapped community.
The three univariate distributions (ALL, TAB and RUG) were analysed to find the best fit models and their corresponding parameters. In this community, none of the univariate distributions analysed exhibited significant CSR or were best modelled by the HETP model (see Table 1). The best fit models for the univariate distributions were single or double Thomas-cluster models (for p d values see Table 1). All the corals (ALL) and tabulate corals (TAB) best fit a DTC model. ALL had two best fit (p d = 0.856) cluster sizes of 24.4 cm and 20.9 cm in diameter (Table 2; note cluster radii are 2r). The TAB best-fit clusters (p d = 0.820) had diameters of 10.6 cm and larger clusters of 154.8 cm ( Table 1). The rugose corals (RUG) best-fit (p d = 0.834) a TC model with cluster diameter of 78.73 cm. The smallest Thomas clusters of the tabulate corals have a notably smaller radius than the rugose Thomas clusters (Table 2), which have large, less densely populated dispersal clusters (Tables 1, 2). These results suggest that the distribution of the corals on the surface is strongly determined by two dispersal episodes for the tabulate corals and one for rugose corals.
BIC analyses of the coral surface areas find two cohorts of tabulate corals and two cohorts of the rugose coral distribution (Dhungana & Mitchell 2021, fig. S2c-d); tabulate corals have a relatively even distribution of specimens between two cohorts (41% and 59%) whereas rugose corals show that the density distribution is dominated by one peak comprising of 81% of the fossils. The cohort of tabulate corals which comprise the majority, are larger than most rugose corals (Dhungana & Mitchell 2021, fig.  S2a-b). The two populations of tabulate corals may correspond to the two Thomas clusters of the DTC model ( Table 2). The dominance of one of the two rugose populations is probably reflected in the best fit single cluster model (Table 1).
For the bivariate distribution of the tabulate and rugose corals (TAB-RUG) there was a non-random distribution (poor fit to CSR); instead, this distribution strongly fits an LDTC model with a p d value close the perfect fit value of 1 (p d = 0.970; Table 1) with the linked clusters centred on the tabulate clusters. The smaller clusters in this model have diameters of 13.0 cm, with the larger clusters of 250 cm (Table 2).

DISCUSSION
Spatial patterns can be complex or subtle, and thus very difficult if not impossible to ascertain by visual inspection alone (Illian et al. 2008). In comparison to forests, there are fewer studies of coral spatial ecology. However, from the work that has been done, spatial distributions of extant corals have been found to exhibit a variety of different spatial patterns including regular spacing (Endean et al. 1997), randomly and aggregated distributions (Muko et al. 2013) and to show remarkable differences in spatial distributions over relatively small (metre scale) distances (Mitchell & Harris 2020). These extant coral analyses can be used to help link observed biological processes to the spatial distributions (e.g. Mitchell et al. 2015Mitchell et al. , 2019. Investigation of dispersal clusters may give insight into dispersal preferences of the rugose and tabulate corals. In our study, the dominant control on the spatial distribution of the corals is dispersal process (Table 1). Extant (scleractinian) corals have preferences for settlement on substrates; the type and composition of attachment substrate can affect the larval metamorphosis rate (Ritson-Williams et al. 2010). Small spatial clusters have been linked to short dispersal distance of the planula larvae around the parental colonies (Muko et al. 2013). Whereas for a spawning genus, the larvae tend to disperse away from the parental populations, leading to much larger radii for dispersal clusters (Muko et al. 2013). It is assumed that most Palaeozoic corals had similar prerequisites concerning the attachment sites of their planula larvae, but little else is known concerning Palaeozoic coral preferences in this context (Scrutton 1998). Given that the tabulate corals on this surface are the facilitator taxa, the larvae of the tabulate corals probably had a greater tolerance for attachment (and survival) on the soft silty substrates of this relatively deep-water environment. The tabulate corals are more highly aggregated (PCF~4 at r = 0; Fig. 4A), whereas the rugose coral distributions have a more diffuse aggregation throughout the spatial scales (PCF > 1) but the extent is aggregation near coral centres is less than half that of the tabulate corals (PCF~1.4 at r = 0, Fig. 4C) and has smaller number of individuals on the surface (Table 1). We tentatively suggest that the tabulate corals on this surface could be dominated by brooding processes due to the high aggregation near coral centres (and hence parental colonies). To confirm this hypothesis, lower classification identification is needed (genus, species) which unfortunately is not possible for this community due to erosion of taxonomically relevant features and relatively low sample numbers. The rugose corals on this surface could be dominated by spawners due to the less aggregated distribution, implying high infant mortality. Alternatively, they could be dominated by brooding processes (and fewer larvae per coral) but have planulae that are more selective regarding attachment sites, thus not as aggregated at small spatial scales. Univariate processes were the main drivers for the spatial distributions as indicated by the lower PCF values of the bivariate interactions (Fig. 4). Nonetheless, the bivariate spatial distribution of the rugose and tabulate corals was significantly non-random, and best-fit to the LDTC T A B L E 2 . Parameters of the best fit models for each of the distributions analysed.

Distributions
Best-fit model Parameter r is a measure of cluster size, with 2r representing the cluster radius. q is a measure of the intensity (density) of the distribution. n is the average number of points in each cluster. Largest clusters are given as the first set of data for each distribution (ALL, TAB). For the TAB-RUG distribution, the first set of parameters are for the tabulate corals and the second set for the rugose corals on the surface.
(p d = 0.970), suggesting a facilitative interaction (Getzin et al. 2006;Lingua et al. 2008). Facilitation is defined as the ability of a resident species to increase the survival of another species, for example by ameliorating the local environmental conditions for another taxon, permitting the recruitment (settlement, growth and survival) of the inferior taxa (Brooker et al. 2008). One such possible mechanism by which facilitation can occur is ecosystem engineering, in which an organism induces a physical state change in the habitat causing a change in the resource distribution (Jones et al. 1997;Wright & Jones 2006). Our analysis identified the LDTC as the best fit model between the tabulate and rugose corals, demonstrating that the rugose coral recruitment or survival was significantly enhanced by the tabulate corals. Modern evidence for facilitation is generally found in harsher environments, where the facilitator is more a stress-tolerant taxon and is providing an otherwise limiting resource (Bertness & Callaway 1994). If the conditions become less hostile then it is possible to see shifts from facilitation to competition (Graff et al. 2007;Maestre et al. 2009). Such shifts cannot be assessed using SPPA using an instantaneous snapshot of a community. However, we found no evidence of competition on the distribution of rugose corals on the surface, even after 15-20 years of growth of the community (Harper et al. 1995), indicating that environmental conditions remained difficult for the rugose corals for this length of time.
Our study demonstrates that not only can SPPA demonstrate facilitation from spatial fossil data but can also distinguish the facilitator and facilitated taxa (in the absence of physical evidence such as encrusting relationships). SPPA can be used to analyse any sessile, in situ, non-time-averaged community in the fossil record since all that is required is specimen positions which reflect the specimens life history (Illian et al. 2008), as in the example of Jurassic crinoids (Hunter F I G . 4 . Best fit models for each univariate and bivariate PCF analysis performed; for p d values see Table 1. A, univariate PCF of the tabulate corals (TAB); grey shading indicates the best fit model (999 Monte-Carlo simulated 'envelope' of a Thomas double-cluster model) for the distribution on the surface. B, bivariate PCF of the tabulate-rugose (TAB-RUG) distribution, with a linked double Thomas cluster simulated 'envelope'. C, univariate PCF of the rugose corals (RUG) and best fit (single double-Thomas cluster) model grey 'envelope'. D, univariate PCF of all the corals on the surface, also best fit with a double-Thomas cluster model. The average value of each best-fit model is given by the dotted line. et al. 2020). Therefore, this study demonstrates that facilitating relationships between fossil taxa can be identified using spatial data. The technique could be applied to a range of sessile taxa, such as sponges, tunicates and brachiopods.
Coral growth in this assemblage has been shown to be hindered by sedimentary influxes (Harper et al. 1995); the larger clusters of tabulate corals on the surface (which are the majority of tabulates) have larger coral diameters on average than all the rugose corals analysed (Dhungana & Mitchell 2021, fig. S2a-b) and probably shielded the rugose corals from the local sedimentary and hydrodynamics conditions, thus increasing the survival of the rugose corals. Amelioration of conditions by one taxon can thereby increase habitable area for the facilitated taxon and so could be considered ecological engineering, as it is likely that the tabulate corals modified local physical conditions.
The evolutionary implications of positive, facilitative interactions are more explicit in some environments, such as reefs where a new habitat is formed by the reef-building taxa (Bruno et al. 2003;Idjadi & Edmunds 2006). However, the data for the role of facilitation in aiding survival of corals is more sparse in the fossil record. Our work provides rare, clear evidence of facilitation between specific orders of Palaeozoic corals, and an insight into their ecological interactions. The identification of tabulate corals as facilitators in soft sediment settings that can enhance coral survival, as we have shown here, could have macroevolutionary implications for the initiation of reef building. In the Middle Ordovician, there was an expansion of metazoan reefs dominated by robust skeletal reef-builders (including tabulate corals) from hard to soft bottom sediments. This expansion coincided with the great Ordovician biodiversification (Kr€ oger et al. 2017) and it has been hypothesized that corals, as physical ecosystem engineers, enhanced diversity during the biodiversification event (Erwin & Tweedt 2012). From the Middle Ordovician into the Silurian, there was also an increase in the abundance of coral-dominated reefs (Kiessling 2009). The ecological drivers for these processes remain unclear, but our study demonstrates that tabulate corals are capable of acting as facilitators in the Silurian and indeed this has the potential to play a role in these expansions. Further investigations into the relationship between reef-expansion and facilitation, using data such as that published by Penny et al. (2020), would enable the role of facilitation in evolutionary trends of specific reef-builders into new habitats to be evaluated with SPPA.
In our study, we have demonstrated how the physical presence of sessile organisms can facilitate other sessile taxa (see Hastings et al. 2007). This facilitation may enable the maintenance and perhaps increase the diversity of sessile organisms by allowing the inferior (facilitated) taxa to survive in a greater range of environments through time. By expanding the niches of taxa, more diverse communities can be sustained and the presence of more (facilitating) sessile taxa can therefore increase diversity. This evolutionary effect has been demonstrated in plant lineages where Quaternary genera in the Mediterranean have ameliorated the environment for 'Tertiary' genera, allowing the older plant species, which are less well adapted to the newer environmental conditions, to survive to millions of years more (Valiente-Banuet et al. 2006;Lortie 2007). Evolutionary parallels are likely to have occurred in metazoan evolution. It is conceivable that similar positive interactions with newer evolutionary fauna have allowed older evolutionary faunas to survive for longer. This effect would effectively drag the older fauna through time, contrasting with the widely accepted displacive view (Sepkoski & Miller 1985;Sepkoski 1996). Ecological interactions structure communities in deep time, and these effects are not restricted to competitive ecological theory but should be observed through the lens of both positive and negative interactions. We have demonstrated that these ecological techniques can be used for in situ, non time-averaged metazoan assemblages of sessile organisms in the Phanerozoic, and it is useful to do so in terms of deciphering complex patterns that conventional palaeontological analyses or visual inspection alone cannot resolve. Better understanding of the interactions, reproductive strategies and substrate preferences of extinct coral genera through techniques such as SPPA could help us inform future, long-term predictions of living coral ecology and survivorship as they face anthropogenic change. The evolutionary implications of facilitation as an ecological mechanism remain crucial for understanding diversity through time. Detection of these relationships allows the development of testable predictions about the abundance, occupancy or diversity change in specific taxonomic groups, thereby increasing our understanding of the role of facilitation in ecological and biodiversity change in the fossil record.
Palaeontological Association Sylvester-Bradley grant (PA-SB201802 to AD), Worts Travelling Scholars Fund and the Cowper-Reed Palaeontology fund (Dept. Earth Sciences, Cambridge) to AD and the Natural Environment Research Council (grant numbers 627 NE/P002412/1 and Independent Research Fellowship NE/S014756/1 to EGM). We are grateful for the comments made by the two anonymous referees, who helped improve the manuscript.
Author contributions. AD first conceived the project. AD and EGM collected the field data. AD performed the analyses. Both authors discussed results and prepared the manuscript.

DATA ARCHIVING STATEMENT
Supporting information, data and R code used in the study are available here: https://doi.org/10.6084/m9.figshare.13286621