Evolutionary contingency in lingulid brachiopods across mass extinctions

of taxa within the occupied


In brief
Liang et al. analyze the global morphospace occupation of lingulid brachiopods through the Phanerozoic.Detailed examination of lingulid morphology suggests that the limited morphological and ecological diversity of modern lingulids likely reflects disproportionate winnowing of morphospace occupation and ecological disparity due to extinction.

Global morphological occupation
Morphospace occupation for the order Lingulida (comprising three superfamilies: Linguloidea, Discinoidea, and Acrotheloidea 13 ) through the Phanerozoic was constructed using a global dataset of 454 individual valves, including both hand specimens and small shelly specimens that embody the full spectrum of their generic diversity (Table 1 in Mendeley Data; Figure S1).Results of principal component analysis (PCA) indicate that 98% of the variance can be explained by the first two principal components (PC1 = 91.0%and PC2 = 6.5%) (Figure 1 and Table 2 in Mendeley Data).PC1 captures the overall shell shape, with higher values representing a rounded shell outline, while lower values correspond to a sub-rectangular shell shape, as demonstrated by the transformed thin plate splines (TPSs) (Figure 1).Values for PC2 primarily reflect the relative convexity of the posterior margin, with lower values representing an acute posterior margin and higher values corresponding to a more rounded posterior (Figure 1).With a high percentage of the variance explained by PC1, body size data were collected for each specimen to determine if the shell shape was influenced by body size.Our results show that this relationship is not significant (Figure S2), indicating that PC1 values were not influenced by body size.The morphological gradients indicated by PC1 and PC2 correspond to the general shell shape, with the most taxonomically diverse group, Linguloidea, occupying the majority of available morphospace, while Discinoidea and Acrotheloidea are restricted to a smaller area of morphospace, specifically the area represented by high PC1 and PC2 values (Figure 1).
Occupation of morphospace through the Phanerozoic was visualized as boxplots using PC1 values, as this encompasses the majority of shape variation (91.0%).Both the maximum (represents rounded shell shape) and minimum (represents sub-rectangular shell shape) values for PC1 occur during the Early Ordovician (Tremadocian stage) (Figure 2A).After the Ordovician, PC1 values are relatively unstable (Figure 2A).The overall spread of values for modern taxa parallels those of the Early Ordovician, except that the number of intermediate values in the modern fauna is diminished (Figure 2A).To facilitate the visual examination of morphological variation through time and between groups, morphospaces were plotted in epoch-level time bins (Figures 2B and 2C).The Cambrian is dominated mainly by linguloids and a few acrotheloids that collectively populate almost the entire range of possible morphospace but with only a few representatives with low PC1 and high PC2 values (Figure 2B).In the Ordovician, morphospace occupation becomes much broader, with the new genera of linguloids that appear at this time occupying the area of morphospace represented by low PC1 and high PC2 values, as well as the first appearance of the discinoids, which occupy the area of morphospace represented by both high PC1 and PC2 values (Figure 2B).The end Ordovician mass extinction saw a drastic reduction in the proportion of total morphospace.The acrotheloids go extinct at the end Ordovician mass extinction, while the discinoids survive (Figure 2B).The Permian-Triassic mass extinction further contracted the total occupied morphospace, dividing the survivors into two distinct groups, with the linguloids and discinoids exclusively found in the areas of morphospace represented by low and high PC1 values, respectively (Figure 2B).
The end Ordovician and the Permian-Triassic mass extinctions largely shaped the morphospace occupation of the more diverse linguloids, compared with the discinoids and acrotheloids (Figure 2C).After the end Ordovician mass extinction, the Pseudolingulidae diversified and the Lingulidae first appeared (both are the main subgroups of the superfamily Linguloidea), occupying the morphospace represented by low PC1 values with a sub-rectangular shell shape.Linguloids almost entirely disappear during the Permian-Triassic mass extinction, with the exception of the Lingulidae, which function as a ''disaster taxon'' during this time, 15

Figure 1. Morphospace visualization for lingulid brachiopods (order Lingulida)
A bivariate plot of PC1 and PC2 with associated density curves along each axis of variation is shown.Thin plate splines represent hypothetical projections of the maximal and minimal loadings of shape coordinates along both axes.See also Figure S1.
occupying the morphospace represented by low PC1 and high PC2 vales and persisting until the present (Figure 2C).
Disparity pattern, evolutionary model, and phylogenetic signals Morphometric data were analyzed to investigate covariation and changes in overall morphology through time.Two primary indices, the sum of ranges (SORs) and the sum of variances (SOVs), were used to examine different aspects of disparity. 5,9here are two minor instances of decoupling between our selected measures of disparity (Figures 3A and 3B).The first occurs during the transition from the Ordovician to the Silurian, when the morphological SOR dramatically decreases while the SOV only exhibits a relatively small decrease.This indicates that although the total morphospace occupation was heavily reduced during this time, the relative distribution of taxa within morphospace only decreased slightly.The second instance of decoupling occurs in the modern fauna, where the morphological SOV reached its maximum value while SOR values remained low, reflecting a significant increase in the dispersion of taxa within morphospace (PERMANOVA test of SORs and SOVs, Df = 9, F = 5.9556, r 2 = 0.10772, p = 0.001).Other instances of decoupling possibly occur during the Phanerozoic, but analysis of disparity for the Silurian to the Neogene is challenging due to low generic diversity through this interval.
Based on maximum likelihood analysis, stasis is the bestsupported evolutionary model for the lingulid shell shape, regardless of whether shape data are binned by stage or by period (Table S1).As the end Ordovician and the Permian-Triassic mass extinctions had the most marked impact on total morphospace, we also examined if the evolutionary model had shifted during these time intervals.These results show that stasis remains the best-fitting model (Table S1).Body size data were also analyzed using the same method, and results show that stasis is the best supported model for any time interval (Table S2).Importantly, in all cases, there is no support for a directional trend in the chosen traits or any subset of the study period.
Phylogenetic analyses of the order Lingulida (using both maximum parsimony and Bayesian approaches) were undertaken in order to explore if there is a relationship between shell shape and phylogeny (Figures 3C, 3D, and S3).Modern Lingulidae genera (Lingula and Lingularia), together with the Ordovician pseudolingulides, obolides, and Lingulasma, form a monophyletic group, all with low PC1 values (Figures 3C and S3).Acrotheloids, discinoids, and other linguloids, however, form a paraphyletic group characterized by high PC1 values.For PC2 values, only two members of the Linguloidea, Cambrian Eoobolus and Lingulellotreta, both of which possess a large pseudointerarea, exhibit low PC2 values (Figure S3).The phylomorphospace (principal component morphospace into which a phylogenetic tree is projected) that is defined by the first two axes is shown in Figures 3D and S3F.Linguloids with a sub-rectangular shell shape occupy the left part of morphospace, while linguloids with a rounded shell shape are clustered with acrotheloids and discinoids, resulting in multiple crisscrossing branches.Consequently, the phylogenetic signal is very low and not significant for the Lingulida shell shape (K mult = 0.2757, p = 0.373 for the maximum parsimony tree; PC2 values: K mult = 0.1889, p = 0.618 for the Bayesian tree).Moreover, convergence is not identified in Lingulida groups (p = 0.291 for maximum parsimony tree; p = 0.569 for the Bayesian tree.p < 0.05 suggests a significant level of convergence), indicating that multiple lineages of the Lingulida do not independently converge on the same area of morphospace.

DISCUSSION
7][18] The maximum morphospace occupation for the Lingulida was reached in the Early Ordovician, suggesting they align with this early burst model (Figure 2A).Both the end Ordovician and Permian-Triassic mass extinction events then drastically reduced total morphospace occupation, ultimately resulting in unpopulated morphospace between the two surviving superfamilies, the Linguloidea and the Discinoidea (Figure 2B).The total morphospace occupied by the Discinoidea remained relatively stable through the Phanerozoic; however, the more diverse Linguloidea lineage experienced marked changes (Figure 2C).Although a large number of Cambrian and Ordovician linguloids possess a rounded shell shape (Figure 2C), species of Ectenoglossa (Figure S4C), which first appeared in the Early Ordovician, 13 already possessed the sub-rectangular shell shape that characterizes modern genera such as Lingula and Glottidia.
In addition to this single species, Pseudolingulidae with a sub-rectangular shell shape also first emerged in the Ordovician (Figure 2C) 13 and subsequently increased in diversity during the Silurian, suggesting that the sub-rectangular shell shape that represents the dominant morphology for modern linguloids had already developed during the Ordovician.Linguloids with low PC2 values gradually went extinct after the end Ordovician, and no surviving linguloids possessed an acute posterior margin after the Permian-Triassic mass extinction (Figure 2C).An acute posterior margin is a result of the elongation of the pseudointerarea, as exemplified by the Cambrian epibenthic genera Eoobolus and Lingulellotreta (Figures 3D and S4B), 13 and is interpreted to be associated with a supportive pedicle that attaches onto substrates for increased stability. 19The disappearance of low PC2 values thus possibly indicates the absence of a truly supportive pedicle that is relatively rigid and employed for permanent attachment to the substrate.This is possibly due to the adoption of an infaunal lifestyle, where a supportive pedicle is unnecessary.For example, the modern taxon Lingula, which lives infaunally, possesses a pedicle that trails posteriorly, providing little or no support. 20,21ological and anatomical implications The relationship between morphological and ecological disparity in the fossil record has received comparatively less attention than that of generic disparity and taxonomic diversity. 5,9However, there is a general expectation that the overall morphological form will correspond to an ecological niche and serve as an effective proxy for understanding the functional ecology of a taxon (Figure 4). 4,5,22Considering the relatively ''simple'' lifestyle of brachiopods, shell shape has been and is herein considered to be directly related to the ecology of brachiopods and is thus reflective of their functional ecology. 23,24For example, modern discinoids are epibenthic (Figure 4G), 13 and given the stable morphospace occupation we identify for the discinoids during the entirety of the Phanerozoic (Figure 2B), fossil discinoids taxa were also likely epibenthic.Examples of clustered epibenthic discinoids attached to various substrates from the Ordovician further support this assertion. 25,26By contrast, changes in morphospace we identify for linguloids likely reflect changes in ecology over time, consistent with previous work. 20,21][29][30][31] Modern linguloids consist exclusively of infaunal forms. 20,21,27he exact timing of when infaunal linguloids first appeared has proved difficult to assess, 27 but detailed examination of their internal anatomical changes, combined with shell shape changes identified in our analysis during the Early Ordovician, may provide some clues.][34] Cambrian linguloids have one pair of dorsal vascula media (Figure 4H) that was subsequently lost in Ordovician and modern forms (Figures 4I and 4J). 35The vascula lateralia on both valves become anteriorly sub-parallel in Ordovician and modern forms, which were likely more easily accommodated within their subrectangular shell shape (Figures 4I and 4J). 20,21These modifications likely reflected how water currents circulated within the brachiopod mantle cavity, potentially as a result of the repositioning of the inhalant water current (anterior laterally) and exhalant current (anterior centrally), which would be more suitable for an infaunal lifestyle.The development of the threesiphon marginal setae seen in modern linguloids would also have been integral to a transition to an infaunal lifestyle (Figures 4F and 4J). 21,36Our morphometric results indicate a gradual reduction of the pseudointerarea from the Cambrian to the Ordovician (Figure S4), as indicated by the elimination of individuals with low PC2 values (Figure 2C).The reduction of the pseudointerarea may positively extend the posterior wall, providing a larger area for the visceral cavity and helping with water circulation efficiency, which could be features needed for the transition to an infaunal lifestyle.It thus appears, based on these changes in anatomy and combined with changes in the shell shape, that all the features necessary for an infaunal lifestyle and that characterize infaunal modern linguloids were likely developed by the Ordovician.

Evolutionary implications
By the Ordovician, maximum morphospace occupancy for the Lingulida was essentially filled, with both sub-rectangular shapes and rounded shapes prevalent.Given that this increase in morphospace occupation does not coincide with an increase in the breadth of their ecological niche, 21,24 the observed changes may possibly be due to neutral processes rather than adaptive change.Our inability to find evidence for directional evolution for both shell shape and body size through the Phanerozoic (Tables S1 and S2) suggests that the observed changes in shape or body size over time are unlikely to be due to selection pressures. 37Moreover, phylomorphospace and convergence analyses provide no evidence of significant convergent evolution, indicating that the surviving linguloid lineages did not adapt to an infaunal niche over time.Considered together, these results suggest that the decrease in morphospace occupation through the Phanerozoic for the Lingulida is not representative of an evolutionary trend or a shift in the ecological niche 38,39 but can instead likely be attributed to disproportional winnowing of morphospace occupation and ecological disparity through extinction.The end Ordovician mass extinction had a differential effect on linguloids, disproportionally eliminating epibenthic forms with rounded shell shapes while linguloids with sub-rectangular shells and an infaunal lifestyle survived both the end Ordovician and Permian-Triassic mass extinctions, persisting until the present day.Both morphospace occupation and the life strategies of the epibenthic discinoids remain unchanged across these two mass extinction events. 25,264][25][26][27][28][29][30][31] This is not to say that other evolutionary processes, such as natural selection, did not shape the evolutionary trajectory of the Lingulida.However, subsequent evolution of the group following the end Ordovician and Permian-Triassic mass extinctions was likely constrained by the much more limited morphological diversity of the survivors.
This series of events is consistent with the concept of contingency (as defined by Gould 6 ), where idiosyncratic events, in this case the differential effect of mass extinction, drive an evolutionary outcome rather than deterministic processes such as natural selection. 40Based on these results, modern linguloids are predominately infaunal not because of natural selection or an evolutionary inevitability but because of somewhat random pruning of the evolutionary tree (given mass extinction events represent rare, chance episodes where adaptive change is not possible). 41,42While this pruning process could have potentially removed the necessary genetic variation needed to reoccupy a larger area of morphospace, the lingulid genome has been shown to be actively evolving, duplicating at a rate approximately two to four times faster than that of other lophotrochozoans, 43 and substantial cryptic diversity exists within modern forms. 44As such, the failure of linguloids to re-radiate into the morphospace they occupied prior to the end Ordovician mass extinction may possibly be due to other constraints, such as a lack of ecological opportunities due to incumbency. 18

STAR+METHODS
Detailed methods are provided in the online version of this paper and include the following:

Geometric Morphometric Analyses
Landmark-based geometric morphometrics was used to collect shape data, using one homologous landmark, defined by the maximum curvature at the posterior position, and 28 equidistant semi-landmarks along the ventral borders (Figures S1C-S1E).
Both landmarks and semi-landmarks were generated with the program TpsDig v. 2.31. 47The definition of landmarks and semi-landmarks for analyses of the pseudointerarea followed the method of Zhang et al. 19 (Figures S1F-S1H).Landmarks were then appended to curves, and the sliding direction of the semi-landmarks was determined using TpsUtil v. 1.78. 47A generalized Procrustes analysis (GPA) with a minimized bending algorithm was used to remove the effects of size, location, and orientation of the specimen images, using the two-dimensional landmark data.The shape of each specimen was then plotted using the landmark configuration from the GPA.Body size data (log-transformed shell area) (Table 2 in Mendeley Data) were collected for each specimen to determine if the shape differences between these groups are the result of ontogenetic allometries, and then multivariant regression between body size and the overall shape were conducted using PAST v3.25 (Figure S2). 48

QUANTIFICATION AND STATISTICAL ANALYSIS Morphospace and Disparity Analyses
The morphometric data matrix was analysed using TpsRelw v. 1.70 47 and PAST v3.25 48 to explore potential changes in morphospace and to visualize shell shape using thin-plate splines.A principal component analysis (PCA) was carried out on the aligned, post-GPA landmarks to show the morphological variation through geological time and between lineages.Morphological data were further analyzed to investigate covariation and changes in shell shape through time.Two primary indices were used to summarize different aspects of disparity.Sum of ranges (SOR) captures the total amount of morphospace occupied and thus reflects patterns of overall expansion and contraction.Sum of variances (SOV) measures the dispersion of taxa around the centroid and reflects the density of taxa within the occupied region of morphospace. 49A bootstrap analysis (1,000 replicates with replacement) was applied for disparity calculations, except the groups without enough data.The disparity analyses were carried out using the R package 'dispRity'. 50Significant shifts in morphospace were tested through nonparametric multivariate analysis of variance (NPMANOVA), also known as PERMANOVA, using 1,000 permutations carried out across all axes of variation. 51All analyses were performed using R. 52 Phylogenetic and Comparative Methods 33 characters were coded for 21 taxa, one outgroup taxon and 20 genera of the Order Lingulida (Table 3 in Mendeley Data).The chosen taxa are considered exemplars of the major Lingulida families and subfamilies.For the selected genera, we predominantly include two genera for each family.For those genera that have a long stratigraphical duration, one species was placed at the base of the clade (like Lingularia for family Lingulidae, and Trematis for family Discinidae), and another species that exhibits fully developed characteristics (like Lingula for family Lingulidae, and Discinisca for family Discinidae) of the family.For other genera, like Paterula, Elkania, and Dysoristus etc., that have shorter stratigraphical durations and few stratigraphic occurrences, only the type genus was selected to include in the phylogenetic analysis.Two recently reported genera, Aulonotreta 53 and Neobolus, 54 were added to the data matrix, modified from Holmer and Popov. 13The phylogenetic data matrix was built in Microsoft Excel 2016.Phylogenetic trees were inferred using both maximum parsimony and Bayesian methods.Parsimony analysis was performed using PAUP* (v.4.0a169). 55A non-parametric bootstrap search based on 1,000 replicates was conducted using a heuristic search algorithm, with starting trees built using stepwise addition and branch swapping undertaken using tree bisection and reconnection (TBR).Results of this bootstrap analysis were summarized as a 50% majority rule consensus tree (Figure S3A).Bayesian analyses were run using MrBayes (v.3.2.7) 56 and the Mkv model, 57 with gamma-distributed rate variation and variable coding.The analysis used a sampling frequency of 1,000, two concurrent runs, four Metropolis-coupled chains, and was run for 10 million generations.A 25% relative burn-in was implemented for all summary statistics.The resulting phylogenetic tree is presented in Figure S3C.Shell shape was not included in the phylogenetic analysis, and its distribution on the tree is therefore not predicted by the tree topology.
To test if closely related taxa tend to have a more similar shell shape to each other than more distantly related taxa, the possibility of a phylogenetic signal in the shape data was calculated using the K mult statistic, a method that measures the similarity of trait values in relation to a Brownian motion model of evolution. 58The K mult value was calculated using physignal function from the R package 'geomorph'. 59A K mult value greater than 1 implies that closely related species are more similar than expected under a null model of Brownian motion (more phylogenetic signal), whereas a value that is less than 1 suggests that relatives resemble each other less than expected (less phylogenetic signal).The significance of the K mult value was then evaluated through a permutation test (1,000 iterations) that randomized the data across the tips of the reference phylogeny.The time-calibrated branch lengths were added to the maximum parsimony and Bayesian tree using bin_timePaleoPhy function with 'basic' model from the R package 'paleotree'. 60In order to visualize the distribution of the shape variables on the phylogenetic tree, the first two PC scores were mapped on the phylogeny using the contMap function from the R package 'phytools'. 61o explore the evolutionary history of shell shape diversification, a phylomorphospace approach was used. 62The phylomorphospace plot was produced by the phylomorphospace function in the R package 'phytools'. 61This function graphs the phylogeny in a two-dimensional space and reconstructs ancestral states along morphospace axes.Phylomorphospace thus not only illustrates current patterns of shell shape variation among taxa, but also reconstructs the evolution pathways that have led to these patterns.To further explore if separate lineages independently converge on a similar area in morphospace (convergence), a result that would indicate that ecological/environmental pressures are a primary driver of patterns of morphological change, we used the R package 'convevol'. 63We applied the function convnumsig, which simulates evolution under a Brownian motion model, to our genus-level tree, using parameters from the first two principal components from our PCA analysis. 63The simulated values were compared with the convnum metric to determine the significance of the convergence.This metric refers to the number of times a lineage has entered a region in a given morphospace.Convnumsig performs 1,000 evolutionary simulations along our Lingulida phylogeny using parameters derived from the observed shell shape data based on the first two PC scores.The convnumsig function requires the number of species selected for comparison within a given group to be larger than the number of variables, in this case, PCs, which is the case for our dataset.
To statistically test for the existence of a trend in shell shape or body size, we fit several evolutionary models of trait evolution, using the R package 'paleoTS'. 64Evolutionary models were analyzed with 'joint' parameterization to examine support across four statistical models of shell shape (represented by PC1 value) and body size change: unbiased random walk, directional trend, stasis, and strict stasis (Tables S1 and S2).A directional trend is the expected model when selective pressures favor shell shape increase or decrease to a constant extent over time, whether due to species sorting or within lineage trends.An unbiased random walk is the null model in the absence of any selective pressure on shell shape. 65A strict stasis model assumes that the mean of a trait is constant through time.Stasis is expected when shell shape is optimized to invariant biological, ecological, and environmental factors. 66Statistical support was compared across these three models using the small-sample corrected version of Akaike's information criterion (AICc) and the associated Akaike model weight. 67Lower scores for the Akaike information criterion and higher Akaike weights indicate stronger support for the respective model.Each interval within a given model include at least six geologic stages long which can constitute a true, long term evolutionary model.

Figure 2 .
Figure 2. Morphospace occupation from Cambrian to present (A) Boxplot of PC1 scores based on global dataset and stage time binning scheme.The width of the box indicates the sample size.(B and C) Morphospace occupation for the order Lingulida (B) and superfamily Linguloidea (C) through time and across mass extinctions.The first two axes (shown) for morphological datasets represent 91.0% and 6.5% of the variation, respectively.Each plot represents all taxa present within a time bin.End Ordovician and Permian-Triassic mass extinctions are marked by red lines.

FigureFigure 4 .
Figure Disparity and phylogenetic comparative analysis(A and B) Disparity measured as the sum of ranges (SORs) and the sum of variances (SOVs) for the Lingulida based on 12-epoch time binning.(C) Mapping of PC1 scores on the time-calibrated phylogenetic tree; the warmer colored clades indicate a higher PC1 score and vice versa.(D) Phylomorphospace of the shell shape (represented by PC1 and PC2 values) for the Lingulida based on maximum parsimony tree.See also Figure S3 and STAR Methods.

TABLE
d RESOURCE AVAILABILITY B Lead contact B Materials availability B Data and code availability d EXPERIMENTAL MODEL AND SUBJECT DETAILS B Taxonomic Scope and Dataset Assembly d METHOD DETAILS B Geometric Morphometric Analyses METHOD DETAILS