Permian–Triassic phylogenetic and morphologic evolution of rhynchonellide brachiopods

Abstract The Rhynchonellida is a major group of brachiopods that survived the “big five” mass extinctions and flourished after the Permian/Triassic (P/Tr) crisis. However, phylogenetic and character evolution in the Rhynchonellida across the P/Tr transition is poorly understood. In view of the widespread homoplasy across this order, we employ a tip-dated Bayesian analysis to reconstruct phylogenetic relationships for late Permian–Triassic rhynchonellides. The same data were also analyzed using three other methods: undated Bayesian, equal-weighting, and implied-weighting parsimony. Compared with trees generated by other methods, those constructed by tip-dating best account for the homoplasy in this group and are closer to previous assumptions on the evolution of this order. Based on the analyses of multiple trees, the major increase in lineage richness occurred in the Early and early Middle Triassic. Also, richness in the Anisian almost reached the highest level seen in the Triassic. According to fossil records, a pronounced reduction in shell size and in the development of ornamentation occurred after the P/Tr extinction, which is largely due to the loss of large and highly sculptured genera and the diversification of small-sized and weakly ornamented genera. Ancestral-state estimation of shell size and development of ornamentation, coupled with comparisons of other characters, indicate that the Early–Middle Triassic mature “small-sized taxa” may have characters displayed by juveniles of their ancestors. This suggests that for these genera, paedomorphosis was possibly a strategy to survive and diversify in the harsh environment after the P/Tr extinction.


Introduction
The rhynchonellide brachiopods, originating in the Early Ordovician, are geologically the oldest and phylogenetically most basal of the extant rhynchonelliforms (Carlson 2016). They survived all "big five" mass extinctions and still inhabit modern oceans Curry and Brunton 2007;Schreiber et al. 2013Schreiber et al. , 2014Carlson 2016). One of the most important evolutionary changes occurred across the Permian/Triassic (P/Tr) boundary. Unlike some orders (e.g., Orthida, Productida, and Spiriferida) that went extinct during the crisis, the Rhynchonellida survived the P/Tr mass extinction and became the second most diverse order (after the Terebratulida) in the Mesozoic and Cenozoic (Lee 2008;Carlson 2016); they, however, were generally a more minor component of brachiopod communities during the Paleozoic (Carlson 2016). When compared with the relatively better-known, cladistically based phylogenies of other orders within the Brachiopoda (Carlson 1991a(Carlson ,b, 1993(Carlson , 1995Alvarez et al. 1998;Carlson and Leighton 2001;Jaecks and Carlson 2001;Carlson and Fitzgerald 2008;Harnik et al. 2014;Congreve et al. 2015;Lee and Shi 2016;Guo et al. 2020b), the phylogenetic evolution of this successful group across the P/Tr transition still remains poorly understood. Xu (1990) published the first cladistic analysis of some Triassic genera, based on relatively few characters, and subsequent research confirmed that morphology-based phylogeny is a powerful tool revealing the evolution of this and other groups (Carlson 1991a(Carlson ,b, 1993(Carlson , 1995Carlson and Fitzgerald 2008;Congreve et al. 2015;Sclafani et al. 2018; Guo et al. 2020b). Accordingly, phylogenetic analysis of Permian and Triassic rhynchonellides provides new insights into the successful evolutionary strategies of the Rhynchonellida across the P/Tr transition.
This paper aims to provide phylogenetic analyses of the late Permian to Triassic rhynchonellides at the genus level, coupled with morphologic analyses to investigate successful strategies and morphologic selectivity during their diversification following the P/Tr mass extinction. Shell size and ornamentation are two basic traits of brachiopods and other invertebrate shells that play important roles in the evolution of these animals (Payne 2005;Vörös 2010;Zhang et al. 2015;Schaal et al. 2016;Wu et al. 2019). Combining these characters and the trees generated, it is possible to investigate the role of heterochrony in the rhynchonellides during this critical period, which may provide useful information on the diversification dynamics of this clade (Gould 1977;McNamara 2012;Schreiber et al. 2013).
Traditionally, to reconstruct phylogenetic trees, morphologic characters are analyzed using maximum parsimony and all characters are weighted equally. Considering the homoplasy in data, some parsimony approaches that rescale characters in relation to their homoplasies have been developed (e.g., implied weighting; Goloboff 1993;Goloboff et al. 2008). Recently, model-based methods, such as Bayesian analysis using the Mk model, have become increasingly popular for phylogenetic analysis and are also important approaches for evaluating biotic parallel or convergent evolution (Wagner and Marcot 2010;Wagner 2012;Wright 2017Wright , 2019. Some simulated studies have indicated that Bayesian inference provides more accurate trees than maximum parsimony (Wright and Hillis 2014;O'Reilly et al. 2016Puttick et al. 2017Puttick et al. , 2019, implying that this inference could be the default method for phylogenetic estimation from phenotype datasets (Puttick et al. 2017). Despite many studies, the best method for phylogenetic analysis among those approaches is still not universally agreed upon (Wright and Hillis 2014;Congreve and Lamsdell 2016;O'Reilly et al. 2016Puttick et al. 2017Puttick et al. , 2019Goloboff et al. 2018;Smith 2019;Barido-Sottani et al. 2020;Keating et al. 2020;Mongiardino Koch et al. 2021).
Phylogenetic relationships and divergence times are critical for macroevolutionary studies of organisms. In the past, approaches to build time-calibrated trees involved estimating the topology and branch lengths in separate and sequential analyses (Bapst 2014;Bapst and Hopkins 2017). However, the temporal data of taxa include important information on the evolutionary dynamics of the group, and the traditional "node-dating" method may discard key information (Ronquist et al. 2012a). By incorporating various sources of information from fossil records, the tip-dating method can deal with this problem and infer both the topology and divergence times simultaneously while accounting for their uncertainties in a coherent Bayesian statistical framework (Ronquist et al. 2012a). This method has been successfully applied to many fossil groups (Bapst et al. 2016;Matzke and Wright 2016;Wright 2017; Zhang and Wang 2019; King and Beck 2020) but has not been used for phylogenetic analysis of rhynchonellides. In addition, for analyzing morphologic and stratigraphic data simultaneously, tip-dating is an ideal tool to study the data that include the convergent evolution of taxa with temporal gaps (Lee and Yates 2018;King and Beck 2020). In the Rhynchonellida, homoplasy is common and the same external appearance can be repeated several times in the evolutionary history of the group (Ager et al. 1972;Cooper 1972;Manceñido and Owen 2001;Savage et al. 2002). Therefore, in this study, tip-dated Bayesian (TB) analysis is employed to reconstruct phylogenetic trees for the late Permian and Triassic rhynchonellides. To investigate the effect of tipdating on tree topology, the same data are also analyzed using other methods: equalweighting parsimony (EW), implied-weighting parsimony (IW), and undated Bayesian (UB).

Occurrence and Duration Data
The occurrence data and stratigraphic ranges of genera were downloaded from the Paleobiology Database (paleobiodb.org). Some doubtful records were amended or excluded. The generic durations based on original descriptive papers and recently published monographs were revised (Shen et al. 2017;Sun et al. 2017). Some newly established genera, not included in Treatise on Invertebrate Paleontology were also added. In the following analysis, the stage bin was used for time intervals. When a genus occurred before and after a stage in the actual fossil record, then it was also assumed to be present within that stage.

Taxon Selection and Character Coding
Taxon Selection.-Unlike higher-level phylogenetic analyses that usually involve a restricted number of taxa in each family/subfamily, the genus-level analysis includes more lower-rank morphologic variation and also helps test the reality and integrity of familylevel groups. For late Permian rhynchonellides, the superfamily Lambdarinoidea is unique in both its external and internal characters when compared with other more "normal" rhynchonellides. Only one species of this superfamily was reported from the upper Permian (Grant 1988;Baliński and Sun 2008); this superfamily therefore is excluded from the present analysis.
Previous studies have shown that a small number of characters selected from a large number of taxa may generate poorly resolved consensus trees (Puttick et al. 2017;Schrago et al. 2018;Barido-Sottani et al. 2020). In addition, computation time may be greatly increased if too many taxa are included. Therefore, the genera whose internal characters (especially the crura, which is a very important structure in traditional classification of this group) are not well known were excluded from the present analyses. Coledium was also deleted, because the species (Coledium erugatum) that perfectly confirms the diagnosis of this genus occurs in the lower Carboniferous and has a large time gap between itself and other coded taxa. The Griesbachian (early Induan) genus Meishanorhynchia is one of the earliest representatives of Mesozoic-type rhynchonellides and is crucial in understanding the initial recovery of this clade after the P/Tr mass extinction, but the anterior parts of its crura remain poorly defined (Chen et al. 2002). However, some newly obtained material from the type locality (in Meishan, eastern China) of the type species of Meishanorhynchia reveals that its crura are laterally compressed and almost straight anteriorly, typical of spinuliforms (Supplementary Fig. S1; see also Supplementary Text). A total of 71 genera (including Meishanorhynchia) were coded for analysis. The type species of most genera were coded based on original and updated taxonomic descriptions. Where some important characters are unknown in the type species, a well-described species that has the same generic affiliation was selected. The list of genera coded and uncoded is provided in the Supplementary Material.
Characters.-The morphologic terms follow Williams et al. (1997) and Savage et al. (2002). Fifty-seven discrete characters were employed in the analysis (Supplementary Data). Of these, 24 characters describe the external appearance of shells, including outline, sulcus and fold, valve convexity, ornamentation, ventral umbo, pedicle opening, and stolidium. The remaining 33 are associated with internal structures and shell perforations; seven of them are related to crura. Characters could be binary or multistate and were treated as unordered in the analysis. Most characters adopted here are comparable with those used by Schreiber et al. (2013). However, some characters used in Schreiber et al. (2013) are not adopted in the present analysis. This is because the latter study is focused on living species, and some characters (e.g., sizes of muscle scars, shapes of teeth, development of socket ridges) can be observed directly from living shells, but these characters may not be observed in fossilized brachiopods whose internal structures are usually reconstructed by a means of serial sections.

Tip-dated Bayesian Inference Analyses
Ages of tips are very important information in these analyses, and the stratigraphic interval of each tip must therefore be determined before the analyses. Here we employed two approaches: "species-dating" and "genusdating" to calibrate the tips (see Supplementary Data for tip dates used in this study), and the morphologic character data were analyzed using both approaches. In the first approach, a tip was dated based on stratigraphic occurrences of selected species from which the morphologic data were derived. For instance, the characters of Abrekia were summarized based on those of its type species Abrekia sulcata described by Dagys (1974) from the Induan in the Russian Far East region, and the tip date of Abrekia therefore is 251.9 to 251.2 Ma, the occurrence range of A. sulcata. This method has the highest levels of accuracy and precision in estimating divergence times (Püschel et al. 2020), although this is not the main goal of our study. In the second tip-calibration approach, the species of a genus were assumed to have the same character states. The age of a tip was given based on the first geologic stage in the range of a genus, and the Wuchiapingian was treated as the tip date of late Permian taxa. For example, Abrekia has a range from the Induan to the Anisian, and it first appeared in the Induan; its tip age is therefore given as 251.9 to 251.2 Ma. Tip ages of a taxon are sometimes identical for these two methods (e.g., Abrekia), but they can be different as well. In both approaches, tip ages were assigned uniform priors over the range of uncertainty .
TB analysis was performed using BEAST v. 2.6.2 (Bouckaert et al. 2019). The Mk model (Lewis 2001) was used, with a gamma distribution to account for rate variation across sites. Characters were partitioned according to the number of character states. The clock model was an uncorrelated lognormal clock (Drummond et al. 2006), and the tree prior was a sampled-ancestor fossilized birth-death model (Gavryushkina et al. 2014). Because the evolutionary history of the rhynchonellide brachiopods is poorly studied, the parameters thus varied within a relatively wide range; the distributions of parameters are listed in Supplementary Text. The analysis was run for 800 million generations, sampling every 200,000, saving 4001 trees. The first 50% of trees (2000 trees) were discarded as burn-in to ensure that the chains had reached a stationary condition. Convergence of four independent runs was assessed based on effective sample size values >200 and confirmed in Tracer v. 1.7 (Rambaut et al. 2018) and the R package RWTY (Warren et al. 2017). The default maximum clade credibility (MCC) consensus tree method was applied to summarize the posterior sample of trees (there are four runs and 2001 trees are retained for every run, so the posterior sample has 8004 trees). Additionally, another consensus tree method, the 50% majority-rule consensus (MRC) tree was also deployed, using PAUP* v. 4.0a166 (Swofford 2003). Previous studies showed that the MCC tree is probably not a good method to summarize the results of Bayesian phylogenetic analyses of morphologic data, because it often includes more incorrect clades than the MRC tree (O'Reilly and Donoghue 2018). However, although the MRC tree has higher accuracy, it is usually poorly resolved (especially for small datasets like ours; Puttick et al. 2017;Schrago et al. 2018), and it is almost impossible to capture much useful information if there are substantial polytomies that cannot be used in the character analysis. The fully resolved MCC tree has more incorrect nodes than the MRC tree, but at the same time, it possibly contains more correct nodes (O'Reilly and Donoghue 2018). In this study, many nodes recovered by the MCC tree are not well supported (see "Results"). To avoid utilizing the MCC tree as the only input in the following analysis, multiple posterior trees were also analyzed to account for the uncertainty of the tree topology and divergence times (Wright et al. 2015Bapst et al. 2016;Soul and Wright 2021).

Lineage Richness
Although the trees generated by TB include information on branch lengths and divergence times, every taxon is considered as a point occurrence in time, and therefore, the trees do not display the stratigraphic range of the taxa (Bapst et al. 2016). Also, during the analyses, the presence of mass extinction was not considered, and many Triassic lineages were dated back to the Permian ( Supplementary Fig. S5). Previous studies on empirical datasets have shown that the tip-calibrated analyses tend to inaccurately recover old divergence-time estimates, and the dating approaches can affect the results of downstream analyses (Ronquist et al. 2012a(Ronquist et al. , 2016Bapst 2014;O'Reilly et al. 2015;Bapst et al. 2016;Matzke and Wright 2016;Bapst and Hopkins 2017;Püschel et al. 2020;Simões et al. 2020). Therefore, to calculate lineage richness over time, we adopted the topology of the tip-dated trees and recalibrated them using the cal3 method (Bapst 2013) in the R package paleotree (Bapst 2012). This approach uses the probabilistic model of branching, extinction, and sampling processes to date samples of trees, and it stochastically draws divergence dates given a set of rates for those processes (Bapst 2013). Simulation studies imply that the cal3 method performs better than other posterior time-scaling methods (Bapst 2014;Bapst and Hopkins 2017). The three rates (branching, extinction, and sampling rates) needed in cal3 were calculated based on the BEAST2 posterior estimates, and their median values were selected. The time ranges of genera are in discrete stratigraphic intervals; therefore, the function bin_cal3Time-PaleoPhy was used, and the first and last appearance dates of each taxon were placed within their first and last intervals under a uniform distribution. This method was applied to both species-dated and genus-dated trees, although the values input for the three rates for species-dated trees may be less accurate, as tip ages in the species-dated analysis are not identical to the first appearance of genera. One thousand dating replicates were performed for each of the two MCC trees, respectively, and one median curve for each tree was calculated. In addition, we also calculated lineage richness for all the posterior samples of trees (8004 trees) generated under both the tipcalibration methods (species-dating and genusdating). Every tree was calibrated once, and all 8004 curves were plotted on the same figure.

Character Analysis
Shell Size.-The shape of the rhynchonellide shell was approximated as an ellipsoid, so the formula 4π × 3 −1 × (0.5L) × (0.5W ) × (0.5T ) was applied to calculate the shell volume of each genus, where L is the shell length, W the width, and T the thickness. For each coded species, the largest specimen in the original or emended descriptions was measured. Then, the shell volume data were normalized by log 10 transformation (Supplementary Data).
Ornamentation.-The development of ornamentation was evaluated by the ornamentation index (OI). The product of "ornamentation length," "ornamentation coarseness," and "ornamentation strength" was used to describe the OI. The ornamentation length describes the distribution of ribs. A rib having a full length was given a value of 1, while a smooth shell was given a value of 0. A rib arising slightly anterior to the umbo, beginning at midlength, and confined to the anterior margin was given values of 0.8, 0.5, and 0.2, respectively. Ribs confined to the sulcus were coded as 0.5. The ornamentation coarseness is measured and coded as follows: coarse and fine plicae (number of ribs <20) possess a value of 1, and very fine ribs (number of ribs >20) are coded as 0.5. Antidichotomous ribs (fine costae that merge anteriorly to form coarse plicae) were treated as coarse plicae, and only the anteriorly coarse part was considered in the OI calculation. The ornamentation strength has three values: 0, 0.5, and 1, which represent absent, faint (e.g., Cyrolexis), and distinct ornamentation (e.g., Rhynchonelloidea), respectively. The OI values therefore range from 0 to 1 (Supplementary Data). To sum up, the longer, coarser, and stronger ribs have the greater OI values, and when ribs are absent, the OI value is 0. Although the OI has discrete values, it was treated as a continuous character in this analysis.
Evolutionary Trend of Shell Size and OI in Fossil Record.-To investigate the evolutionary trends of size and OI, we calculated shell size and OI values for all genera from the upper Permian and Triassic based on their type species. Variations in these values were calculated for each stage from the late Permian to Triassic based on their durations from fossil records; the structures of trees are not considered in this analysis. For shell size, the Mann-Whitney test was applied to evaluate whether there is a significant change in shell size between two adjacent stages. The variation of ornamentation development was represented by the mean of the OI.
Ancestral-State Estimation.-To investigate the distribution and evolution of size and OI on phylogenetic trees, ancestral-state reconstruction was performed using the fastAnc function in the R package phytools (Revell 2012), which generates maximum-likelihood ancestral states assuming a Brownian motion model. The reconstructed ancestral states were plotted on the MCC trees using the contMap function. Moreover, we paid special attention to the character evolution of the Early and early Middle Triassic "small-sized taxa" (see "Results"). We compared the character states of the tips (i.e., the small-sized taxa) and those of their ancestral nodes. The most recent ancestral node of a tip was labeled as "AN1," and the most recent ancestral node of AN1 was termed "AN2," and so on. In total, character states of four groups of ancestral nodes (AN1 to AN4) were considered. As mentioned earlier, the branch length or the tree topology may affect the comparative analysis, and we conducted this analysis for both BEAST2-generated (not recalibrated) trees and cal3-recalibrated trees. For both cases, a thousand (1000) trees were randomly sampled from the posterior sample of trees, and the recovered character states of AN1-AN4 of the 1000 trees were recorded and displayed in box plots. In a very small number of trees, some genera do not have AN4 and/or AN3, and the states of those nodes were not considered. All calculations were performed using R (R Core Team 2020).

Phylogenetic Analyses Using Other Methods
Out-Group.-In addition to tip-dated analyses, the same data were also analyzed using EW, IW, and UB. For these analyses, a camerelloid genus, Camerella was chosen as the outgroup (Schreiber et al. 2013). The Camerelloidea belongs to the Pentamerida, but that superfamily was closely allied to the Rhynchonellida Schreiber et al. 2013). Instead of the type species, another species, Camerella bella, was selected to be coded to represent Camerella following Schreiber et al. (2013).
IW and EW Analyses.-IW can down-weight characters according to their homoplasy. The extent of down-weighting under IW is controlled by a concavity constant (k) (Farris 1969;Goloboff 1993). A very low k-value strongly penalizes the homoplastic characters, and IW begins to exhibit the undesirable properties of clique analysis. Nevertheless, when a k-value is very high, IW performs like EW (Goloboff et al. 2008(Goloboff et al. , 2018Smith 2019). In the present analysis, we applied the TNT script setk.run, written by Salvador Arias, to search for the appropriate value of k. The script returned a value of k = 12.63 of our data; that value was adopted in IW.
Parsimony analyses were performed using TNT v. 1.5 (Goloboff and Catalano 2016). For both EW and IW, heuristic searches were performed using traditional search algorithms. Ten thousand replicates were executed, saving 10 trees per replicate. The strict consensus trees were calculated for the most parsimonious trees (MPTs). Bootstrap values were calculated to examine the results of EW analysis, and symmetric resampling support values were calculated for the IW analysis.
UB Analysis.-UB analysis was performed in MrBayes v. 3.2.6 (Huelsenbeck and Ronquist 2001;Ronquist and Huelsenbeck 2003;Ronquist et al. 2012b). The Mk model (Lewis 2001) with gamma distribution priors for site rate variation was used for analysis. The prior on the gamma shape parameter was a uniform distribution between 0 and 10. The analysis was run for 80 million generations with four runs of four chains that sampled every 20,000. Convergence of the four runs was confirmed in Tracer v. 1.7 (Rambaut et al. 2018) and standard deviation of split frequencies <0.01. The first 50% of trees were discarded as burn-in. The 50% MRC tree was constructed for the posterior sample of trees using the sumt function in MrBayes. Moreover, the MCC tree was also calculated using TreeAnnotator in BEAST v. 2.6.2 (Bouckaert et al. 2019). All consensus trees and the character-taxon matrix in Nexus format are provided in the Supplementary Material.

Tree Topology
Tip-dated Bayesian Inference Analyses.-For both calibration approaches, the 50% MRC trees are poorly resolved, as expected (Supplementary Fig. S2). Therefore, only the MCC trees are displayed here (Fig. 1). When tips are calibrated using the species-dating method, the members of the Paleozoic superfamily Stenoscismatoidea form a clade in the basal part of the tree (Fig. 1A). Other Paleozoic taxa and the Triassic genus Lissorhynchia are included in a monophyletic group, which is a sister group of the clade consisting of other Mesozoic taxa. In the larger Mesozoic group (from Homoeorhynchia to Caucasorhynchella), taxa are classified as two poorly supported clades.
The smaller clade includes halorellids, wellerelloids, pugnacoids, some norelloids, and two raduliform genera: Caucasorhynchella and Aorhynchia. Although their crural types and superfamilial affiliations are variable (see Supplementary Data), they all lack developed septal plates and a dorsal septalium. Except halorellids and the two raduliform genera, all other members in this clade are grouped into three subclades: the Batongorhynchia-Laevirhynchia clade, which has weak ornamentation and lacks dental plates; the Austriellula-Trigonirhynchella clade, which features weak dental plates and weak ornamentation; and the Calcirhynchia-Veghirhychia clade, which is characterized by fully and densely costate shells. The other large monophyletic group in this tree consists of the raduliform (or variations thereof) genera and the Dimerellidae. Most of the genera in this group have developed dental plates, septal plates, and a dorsal septalium. Some monophyletic groups are also recognized in the tree, but many of them are not strongly supported.
The MCC tree generated using the genusdating method displays two large clades (Fig. 1B). The lower one consists of wellerelloids, pugnacoids, some norelloids, and members of Halorellidae. Some Paleozoic members grouped in this clade are located in its basal parts. Some members of the Norelloidea (Norella, Austriellula, Batangothynchia, and Laevirhynchia) are grouped in a monophyletic group, together with three wellerelloid genera: Sinorhynchia, Apertirhynchella, and Moisseievia. The other large clade includes members of the Stenoscismatoidea, Dimerellidae, Rhynchonelloidea, Hemithiridoidea, and Rhynchotetradoidea, and some norelloids also appear in this group. Except for the Stenoscismatoidea, which is located in basal positions, three smaller clades are classified: the Carapezzia-Aorhynchia clade, which has variable internal structure and ornamentation; the Fissirhynchia-Timorhynchia clade, which is characterized by an entirely costate shell; and the Herangirhynchia-Paranorellina clade, which has a partially costate shell.
Other Analytical Methods.-Consensus trees generated using the EW, IW, and UB are provided in the Supplementary Material. In all these trees (Supplementary Figs. S3, S4), the Stenoscismatoidea forms a clade in the basal part of the tree. Nevertheless, in contrast to the tip-dated MCC trees discussed earlier, many Paleozoic members are located in relatively derived positions. Raduliform (or variations thereof) genera are often reconstructed as paraphyletic associations (they are included in a monophyletic group in the MCC tree generated by UB; Supplementary  Fig. S4), and the septifal and arcuiform members are included in a monophyletic group. Although topologies are different among these trees and many nodes are not strongly supported, the general pattern of taxonomic and character distributions is highly comparable: raduliform (or variations thereof) genera in the lower part, members of the Dimerelloidea (sensu Savage et al. 2002) in the middle part, and septifal or arcuiform elements in the upper part.

Lineage Richness Variation
Generic richness of fossil data underwent a dramatic drop across the P/Tr boundary ( Fig. 2A) and remained rather low during the Early Triassic (Induan-Olenekian). The Anisian (early Middle Triassic) saw the first increase in generic richness of Triassic rhynchonellides, with the number of genera exceeding pre-extinction levels. Then, generic richness experienced a gentle decline in the Ladinian, followed by a stepwise increase through the Ladinian to the Norian, with a pronounced peak in the Norian.
The lineage diversities calculated based on the two MCC trees are shown in Supplementary Figure S5. Both curves show a gradual increase in lineage richness from the late middle Permian to the early Middle Triassic, although there are some minor fluctuations, and a small drop in lineage richness occurred in the late Permian on the curve based on the genus-dated MCC tree. From the Middle Triassic to the end of the Triassic, the lineage richness varied strongly, and there are several obvious losses and increases in lineage diversity.
The lineage richness calculated using raw tip-dated MCC trees does not consider the last appearances of genera. By contrast, the stratigraphic ranges were added by recalibrating trees using the cal3 method. The lineage richnesses calculated by cal3-calibrated trees are displayed in Figure 3. Although the trees generated using these two calibration approaches (genus-dated and species-dated) are not mutually consistent in topology, the lineage richness variations calculated based on these two methods are very close to one another (Fig. 3). In addition, the general trends of lineage richness outlined by posterior tree samples (8004 trees) and by MCC trees are comparable. Lineage richness derived from those two approaches shows no evident loss across the P/Tr boundary, which is expected, as many Changhsingian genera that disappeared in the P/Tr extinction, were not included in phylogenetic analyses. Instead, a pronounced increase in lineage richness occurred from the Early Triassic to the Anisian. Also, richness in the Anisian almost reached the highest level seen in the Triassic. After the Anisian, the richness in both lineages fluctuated slightly and maintained a high level until the middle Norian. The lineage richness, coupled with generic richness, declined markedly in the Rhaetian (Fig. 3).

Shell Size Variation
Fossil records show that rhynchonellide shell sizes suffered a conspicuous reduction across the P/Tr boundary and reached their lowest values in the Induan (Fig. 2B). The rather low levels of body sizes through the Induan to Olenekian within the Early Triassic interval indicate the miniaturization of rhynchonellide brachiopods in the aftermath of the P/Tr extinction. Then, body sizes increased significantly in the Anisian, followed by a gentle, stepwise increase through the Middle-Late Triassic, but never returned to pre-extinction levels until the Norian (Fig. 2B).
After an estimation of ancestral states, shell sizes and the OI values are plotted on the MCC trees generated by TB analysis (Figs. 4,5,. In both speciesdated and genus-dated trees, there are some small-sized elements in the aftermath of the P/Tr extinction, and most are in the upper parts of the figures, such as Meishanorhynchia, Paranorellina, Nucleusorhynchia, Abrekia, Parabrekia, and Lichuanorelloides. They are not within a monophyletic group, but all appear to have close relationships with some members of the Rhynchonelloidea (sensu Savage et al. 2002). In the following discussion, they are termed "small-sized taxa." Their ancestral nodes usually have a darker color (i.e., larger size) than the small-sized taxa (see Fig. 4), suggesting that the small-sized shell is a derived state. Moreover, there are some taxa of distinctly larger size: Sulcirostra, Rhynchonellina, Carapezzia, Halorella, and Halorelloidea ("largesized taxa" in the following discussion). These taxa belong to the Dimerelloidea (Dimerellidae and Halorellidae; sensu Savage et al. 2002), but the latter is not a monophyletic group according to our results. In some parts of the trees, the descendant lineages seem to have larger shell sizes than their parent lineages. Although there are several exceptions, both the Homoeorhynchia-Paranorellina clade in the species-dated MCC tree and the Herangirhynchia-Paranorellina clade in the genus-dated MCC tree display an increasing trend in shell size along the lineages (Fig. 4, Supplementary Figs. S6, S7).
The character states of the ancestral nodes of the small-sized taxa were estimated across 1000 trees randomly selected from the posterior samples. The analytical results for genus-dated trees and species-dated trees (Fig. 6, Supplementary Fig. S10A,B) are almost identical to one another, so only the former is shown here. For all six genera, there is a downward trend in shell size from the AN4s to the tips in both tip-dated trees and recalibrated trees, and it is more obvious in the three Early Triassic genera Meishanorhynchia, Lichuanorelloides, and Paranorellina. The only exceptions are that in some cases, the sizes of Paranorellina and Abrekia are slightly larger than the median sizes of their respective AN1 and/or AN2.

Ornamentation Variation
The OI measures the development of ornamentation on shells. Clearly, the OI values also experienced a dramatic decline across the P/Tr boundary (Fig. 2C). This decline continued from the Induan to the Olenekian (Early Triassic). The same proxy increased rapidly in the Anisian, followed by a gentle, stepwise increase in the Middle-Late Triassic. The postextinction OI values eventually reached preextinction levels in the Late Triassic. In our dataset, two large groups of taxa have the most prominent ornamentation: fully costate and raduliform elements (mostly hemithiridoids) and fully costate and septifal taxa (mostly wellerelloids and pungnacoids). In addition, there are two weakly ornamented groups: the semicostate and raduliform taxa (roughly equal to the small-sized taxa) and the smooth or slightly ribbed arcuiform and septifal genera (e.g., Norella, Austriellula, Batangorhynchia, Laevirhynchia, Sinorhynchia, and Trigonirhynchella). In the two MCC trees, neither increase nor decrease in OI is obvious in the "arcuiform and septifal groups," but a trend of increasing OI is possibly present in the "semicostate and raduliform group" (Fig. 5,Supplementary Figs. S8,S9): the Early to Middle Triassic small-sized taxa have the weakest ornamentation within the semicostate group, and the later diversified clades have relatively higher OI values. As for shell size, the ancestral nodes of the small-sized taxa also have a darker color, indicating their ancestors may have had more developed costae (Fig. 5, Supplementary  Figs. S8, S9). Similarly, the estimates of the states of their ancestral nodes across multiple trees exhibit a decreasing trend in OI from the AN4s to the tips as well (Fig. 6, Supplementary  Fig. S10A,B), and again, this trend is rather prominent in Meishanorhynchia, Lichuanorelloides, and Paranorellina.

Analytical Methods and Tree Topology
Traditionally, morphologic data are analyzed using parsimony methods. To diminish the effect of homoplasy, the IW approach was introduced (Goloboff 1993; Goloboff et al. 2008). For our dataset, the external characters exhibit substantial homoplasy, and therefore, they were down-weighted more under IW. Several nodes in the consensus trees clearly indicate the effect of IW ( Supplementary Fig. S3). For example, in the EW strict consensus tree, Laevirhynchia is grouped with Uniplicatorhynchia, because they are very similar in external appearance such as having smooth shells and a uniplicate commissure. Nevertheless, their internal characters are very different. Unlike Uniplicatorhynchia, Laevirhynchia lacks dental plates and a septalium and shows a FIGURE 5. Ancestral-state reconstruction of ornamentation index (OI), plotted on the "species-dated" maximum clade credibility (MCC) tree (A) and the "genus-dated" MCC tree (B). Darker color means higher OI and more pronounced ornamentation (online version in color). See Supplementary Material for the same figure with all taxa labeled. closer relationship with Batangorhynchia and Sinorhynchia in the IW tree. The latter two genera also lack dental plates and a septalium. The 50% MRC tree generated by UB is poorly resolved (Supplementary Fig. S4). Even so, general patterns of taxa and character distribution in this tree are also comparable with those demonstrated in the IW and EW consensus trees, as stated earlier. Also, many clades in the MCC tree generated by UB are identical with those of the EW and IW consensus trees. Therefore, although topologies are not completely identical, trees generated using the three methods: EW, IW, and UB are broadly similar (King 2020).
By incorporating stratigraphic information into phylogenetic analysis, TB generates quite different trees from the other three methods, and tip-dated MCC trees obviously have better stratigraphic congruence (King 2020). In tipdated MCC trees (Fig. 1), all the Paleozoic taxa are positioned in the basal part; nevertheless, many of them (mostly wellerelloids) are reconstructed in derived positions in EW, IW, and UB consensus trees (Supplementary Figs.  S2-S4). This discordance of topology and stratigraphic position can be caused by a variety of factors (Carlson and Fitzgerald 2008). Intuitively, it is possible that the choice of out-group determines the position of those wellerelloids. The out-group taxon Camerella has external ornamentation and crura similar to those of the stenoscismatoids and Triassic semicostate rhynchonelloids. Consequently, these taxa are grouped in the basal part of the tree, and the wellerelloids, having very different characters from Camerella, migrated to relatively derived positions in the trees. Also, the earliest members of this group are not included in this study, and character states of the Triassic rhynchonellides may have appeared many times in the long history of this order. In this case, the choice of out-group will be challenging. TB, which does not need an out-group taxon, is more appropriate for the analysis of our data.
Although the IW method can lessen the effect of homoplasy, it does not consider the temporal disparity between taxa. If two species with a large stratigraphic gap develop very similar characters independently (convergent evolution), they will be grouped together in FIGURE 6. Estimates of character states of the "small-sized taxa" (tips) and their respective ancestral nodes (AN1-AN4) based on 1000 trees randomly selected from the posterior samples. A, Results calculated from non-recalibrated "genusdated" trees. B, Results calculated from cal3-recalibrated "genus-dated" trees. Outliers are omitted. IW trees (e.g., the Late Triassic Himalairhynchia and the Paleozoic Rhynchopora; Supplementary  Fig. S3). In contrast, taking account of the temporal gap between these two species, they are positioned in two mutually remote clades in the tip-dated trees (Fig. 1). Thus, for such groups as rhynchonellides, which have a limited number of characters and sometimes display homoplasy, our study indicates that TB maybe a better method than EW, IW, and UB to reconstruct the phylogenetic relationships. In this study, the evolutionary relationships of the rhynchonellides displayed in tip-dated MCC trees are much closer to their actual stratigraphic occurrences in the fossil record, and thus the relationships are more reliable, if we assume that the fossil records of this group are relatively complete and dependable.
However, this advantage of TB may also generate questionable relationships if the fossil records are incomplete. For instance, if two closely related species are discovered from two horizons with a large temporal gap and no comparable fossils are reported in this gap, the tree generated using TB may mistakenly place them in two different clades, especially when the character signal is weak (King 2020). This can be reflected by trees generated by two different calibration approaches. In the genus-dated MCC tree, Nudirostralina and Homoeorhynchia are closely located (tip dates: Nudirostralina, Olenekian; Homoeorhynchia, Carnian), as well as in the EW, IW, and UB consensus trees (Fig. 1B, Supplementary Figs. S3,  S4). However, in the species-dated MCC tree (Fig. 1A), they are far apart from one another, which is possibly caused by the large temporal gap between the coded species of the two genera (tip dates: Nudirostralina, Anisian; Homoeorhynchia, Pliensbachian). Another example is the different position of Norella. It has a long range from the Olenekian to the Rhaetian. In the species-dated tree (Fig. 1A), however, it was positioned in a relatively derived position, because its tip age was given as "Norian" in that analysis. The differences in topology between the genus-dated and species-dated MCC trees indicate that for some taxa, character signals can be overturned by temporal signals (King 2020). Because the temporal gaps between genera are generally smaller in genus-dated analysis when compared with those used in species-dated analysis, character signals played relatively more important roles in the former. Therefore, the MCC tree generated by the genus-dated method is closer to traditional classifications that are proposed according to characters.
It is noticeable that in all the consensus trees (Fig. 1, Supplementary Figs. S3, S4), many nodes are not strongly supported. This is caused by the limitations of this morphologic dataset, which includes many genera with limited numbers of characters and character states, thus, it is difficult to generate consistent topologies for Bayesian or bootstrap analyses. As stated in "Data and Methods," the MCC tree is not a recommended method to summarize a posterior sample of trees (O'Reilly and Donoghue 2018). However, the MCC tree is a tree with the greatest production of clade probabilities, and at least, it represents one of the most credible tree structures; moreover, it provides more valuable information than the MRC tree if the latter recovers rare nodes and clades. The low posterior support of nodes (i.e., highly variable topology) is not uncommon in empirical datasets (O'Reilly and Donoghue 2018; Barido-Sottani et al. 2020), and the tree topology and divergence times affect the results of phylogenetic comparative analysis (Bapst 2014;Bapst and Hopkins 2017). However, the downstream comparative analyses are not impossible if multiple trees, rather than a single point estimate of phylogeny, are analyzed (Wright et al. 2015;Bapst et al. 2016;Soul and Wright 2021).

Lineage Evolution of the P-Tr Rhynchonellides
The stratigraphic ranges of genera were treated differently in the raw tip-dated MCC tree and recalibrated trees (the last appearances were not considered in the former analysis). Thus, the lineage richness in the late Middle to Late Triassic cannot be compared directly, because many taxa that originated in the early Middle Triassic and persisted to the Late Triassic were not counted in the curves based on raw MCC trees. However, an outstanding discrepancy occurs between the "raw" and "recalibrated" lineage richness on the left parts of the curves. In the raw MCC trees, many nodes were dated back to the middle and late Permian, resulting in an increase in lineage richness at that time. By contrast, the nodes recalibrated by the cal3 method seem to have a relatively younger age, thus, the lineage diversity increased frequently in the Early and Middle Triassic. Considering that the P/Tr extinction is the most severe extinction event in Earth's history (Chen and Benton 2012), many Paleozoic lineages were likely truncated across the P/Tr boundary, and most Triassic lineages diversified after the extinction. If this is true, the richness calculated by the cal3 method is more reliable. The tendency of tipcalibrated analyses to recover inaccurately old divergence-time estimates (or "deep root attraction") has been observed in many datasets (Ronquist et al. 2012a(Ronquist et al. , 2016O'Reilly et al. 2015;Bapst et al. 2016;Matzke and Wright 2016;Püschel et al. 2020;Simões et al. 2020). Many factors may have caused the deep root attraction, such as an inadequate model of morphologic evolution and the prior distributions of the parameters (Ronquist et al. 2012a(Ronquist et al. , 2016Püschel et al. 2020). Maybe more sophisticated models or carefully constrained values of parameters are needed to deal with it (Ronquist et al. 2016;Simões et al. 2020).
All the calculated lineage richnesses show an evident increase in the Early and early Middle Triassic and reach the highest level in Anisian. In contrast, generic diversity peaked in the Norian. The lineage diversification of the rhynchonellides in the Early and early Middle Triassic is partially consistent with taxonomic diversification in fossil records. The small-sized taxa appeared early in the ocean after the P/Tr mass extinction (Chen et al. 2002(Chen et al. , 2005a(Chen et al. , b, 2007Wang et al. 2017;Wu et al. 2021), and they represent the earliest Mesozoic-type rhynchonellides. In the Olenekian, in addition to the raduliform (or variations thereof) genera, there are other newly originated genera such as the arcuiform Norella and Lissorhynchia. The Anisian also witnessed a great biodiversification of rhynchonellides (Yang and Xu 1966;Dagys 1974;Pálfy 2003;Chen et al. 2005aChen et al. , 2018Guo et al. 2020a). The Anisian taxa developed a variety of external and internal structures, for example, smooth or semicostate or fully costate shells, possession or lack of dental plates, a septalium and dorsal median septum, and occurrences of various types of crura. Many higher-level classification units in the traditional classification schemes therefore have already appeared in the Anisian. All these data indicate that the diversity and morphologic complexity of rhynchonellides recovered in the Anisian and even exceeded the pre-extinction levels.
Moreover, lineage richness of cal3-calibrated trees decreased gradually from the late Norian to the Rhaetian probably due to the "Signor-Lipps effect" (Signor and Lipps 1982;Wagner 2019). Because the last appearances of genera are uniformly distributed through their last stratigraphic intervals, generic extinctions occur before the end-Rhaetian and gradually increase in frequency toward the Triassic/Jurassic boundary, resulting in a gradual loss in lineage richness. Similarly, the diversification in the Early and early Middle Triassic also may be biased by the "Jaanusson effect" (i.e., the first appearances of genera may be younger than their true origin time, making the pattern of very rapid diversification appear to be gradual; Jaanusson 1976) to some extent, which usually delays the recovery of diversity. Regardless, these results indicate that the lineage richness reached almost the highest level before the end-Anisian.
It is noteworthy that we used the relatively simple fossilized birth-death model in this study. However, it is highly possible that the evolutionary dynamics of the rhynchonellides differ during the extinction, radiation, and background intervals. Therefore, further efforts based on more complicated models (such as the skyline model; Stadler et al. 2013;Gavryushkina et al. 2014) and informative priors should be made in the future (Simões et al. 2020;Wright et al. 2021). Stratigraphic range data of taxa may also be retained in analyses with the uncertainties of ages considered (Stadler et al. 2018;Barido-Sottani et al. 2019. If the divergence times are properly estimated, there will be no need to recalibrate trees using posterior approaches, which is the original goal of this method.

Shell Size and Ornamentation Evolutions
Shell size was not included in the phylogenetic analysis, and its distribution on trees therefore is not predicted by the tree topology. As stated earlier, some genera are conspicuous in all consensus trees in terms of shell sizes (Fig. 4, Supplementary Figs. S6, S7). Of these, the large-sized brachiopods belong to the Dimerelloidea (sensu Savage et al. 2002); they are not a monophyletic group on our trees. These taxa are unusual and were thought to have inhabited hydrocarbon seeps and hydrothermal vents (Campbell and Bottjer 1995;Peckmann et al. 2007Peckmann et al. , 2011Peckmann et al. , 2013Sandy 2010;Kiel et al. 2014). The specialized hydrothermal habitats offered sufficient nutrients for dimerelloids to grow large (Kiel et al. 2014). The apparent miniaturization of shell size in Early Triassic was accentuated by the disappearance of the large Permian genera and the diversification of small-sized individuals (Fig. 4). Some parts of the MCC trees also show the increasing shell sizes along lineages, as noted in the "Results" (Fig. 4, Supplementary Figs. S6, S7), but the mechanism of enlargement is poorly understood.
With respect to the variation of ornamentation, a sharp drop in OI coincided with the P/ Tr extinction (Fig. 2C). Within the Early Triassic, most of the newly originated elements (e.g., Abrekia, Laevorhynchia, Lichuanorelloides, and Meishanorhynchia) have smooth or semicostate shells that are lower in OI compared with fully costate genera, which resulted in the low OI value in the Early Triassic. The reason why the Induan saw a slightly higher OI than that in the Olenekian is that a Permian fully ribbed taxon (Terebratuloidea) temporarily survived the P/Tr extinction (Chen et al. 2009) and increased the mean OI of the Induan. For some completely smooth or costate groups, closely related taxa may have very different ornamentation (e.g., ribbed Halorella vs. smooth Halorelloidea, capillate Rhynchonellina vs. ribbed Sulcirosta; Fig. 5); the development of ornamentation does not correlate significantly with shell size. For semicostate groups, however, it is suggested that the development of ornamentation varies during ontogeny, with larger individuals often having longer and more prominent plications than smaller ones (Cooper and Grant 1976;Wang et al. 2017; Fig. 7). This is also observed on our trees: the variations of shell size and OI of the superfamily Rhynchonelloidea and its affinities are highly mutually consistent (Figs. 4,5). Brachiopod ornamentation is thought to have prevented predation (Leighton 1999;Vörös 2010). Strong external sculpture is an important protective mechanism for Mesozoic brachiopods, and the escalation of predators was also reflected by the increased ornamentation of brachiopod clades during the Mesozoic marine revolution (Vörös 2010). However, considering the different patterns of ornamentation development in various clades, more detailed studies are needed to investigate the relationships between OI and predation pressure.

Paedomorphosis in the Small-sized Taxa
Paedomorphosis is a pattern of heterochrony wherein growth is retarded during ontogeny in descendants, compared with their ancestors, and is common among extinct and extant animals, including brachiopods (Gould 1977;MacKinnon 2001;McNamara 2012;Bitner et al. 2013). The occurrence of the small and smooth Meishanorhynchia in the aftermath of the P/Tr extinction was assumed to be an example of paedomorphosis because of its small size and spinuliform crura (Manceñido and Motchurova-Dekova 2010). In fact, in addition to Meishanorhynchia, other small-sized taxa are also possibly paedomorphic taxa. These taxa were traditionally assigned to two superfamilies: the Rhynchonelloidea and Norelloidea. Our analysis does not support a monophyletic group, but they appear to be closely related to some members of the Rhynchonelloidea (sensu Savage et al. 2002). As members or close relatives of the Rhynchonelloidea, the small-sized taxa display a smaller size and weaker ornamentation compared with larger relatives. The estimated ancestral states show that in general, the ancestors of the small-sized taxa have a relatively larger size and more pronounced ornamentation than their descendants, with these features being more prominent in the Induan Meishanorhynchia and Lichuanorelloides (Figs. 4-6). As mentioned earlier, Paranorellina may have slightly larger size than its AN1. This is because it is always recovered as a sister taxon of Meishanorhynchia, which has a smaller size than Paranorellina, but appeared earlier than the latter. The estimated size of their AN1 (they have the same AN1) therefore is closer to that of Meishanorhynchia and smaller than that of Paranorellina. In addition to a small size and weak ornamentation, all these taxa have a depressed profile, a low fold and shallow sulcus, low convexity, and spinuliform to incipiently raduliform crura, which are often regarded as features of the juveniles of some strongly convex rhynchonellide genera (Ager 1962;Dagys 1968;Wang et al. 2017 ; Fig. 7). These characters, together with the reconstructed ancestral states for shell size and OI, imply paedomorphic origins for the small-sized taxa among the Triassic rhynchonellides.
Paedomorphosis is induced by progenesis, postdisplacement, and neoteny (McKinney and McNamara 1991;McNamara 2012). Any one of these processes or their combinations can result in paedomorphosis. Usually, only when life histories of taxa are well known and compared, can a specific cause be detected (McKinney and McNamara 1991;Jaecks and Carlson 2001). Nevertheless, for many fossil taxa such as rhynchonellide brachiopods, it is not easy to judge which process has played the most important role in the formation of the heterochronic pattern due to the incompleteness of fossil preservation. Paedomorphosis can be stimulated by many extrinsic factors, for example, temperature, predation pressure, habitats, and nutrient supply (McNamara 1983;McKinney and McNamara 1991). Meishanorhynchia appeared rapidly in the aftermath of the P/Tr extinction (Chen et al. 2002(Chen et al. , 2007, and other small-sized taxa were already diverse in the Early and Middle Triassic. Accumulating evidence shows that some harmful factors, such as high seawater temperature, limited nutrition, low oxygen content, and microbial bloom associated with the P/Tr extinction, recurred in Early Triassic oceans (Erwin 2006;Knoll et al. 2007;Chen and Benton 2012;Sun et al. 2012;Song et al. 2014;Chen et al. 2015Chen et al. , 2019Huang et al. 2017Huang et al. , 2019. For rhynchonellide brachiopods, lack of food, low oxygen, and low carbonate saturation in seawater may limit their ability to produce calcium carbonate to form shell substance, resulting in small-sized shells with weak ornamentation. To survive and develop in the inhospitable, unstable, and unpredictable environment of the Early Triassic, it is possible that these animals had to mature rapidly to reproduce and shorten their life spans, essentially in "r-selection" mode (Reznick et al. 2002). These processes also indicate that progenesis is probably the primary trigger of paedomorphosis in this group. More confident conclusions are possible if their "ancestors" are more fully FIGURE 7. Possible ontogeny of members or close relatives of the Rhynchonelloidea. Some species are listed here, and the rectangular shape indicates ontogenetic stages observed in the fossil record (Ager 1962;Dagys 1968Dagys , 1974Chen et al. 2002;Wang et al. 2017;Guo et al. 2020a). Note that not every species undergoes all these developmental stages. M., Meishanorhynchia.
sampled, and their life history can be more completely studied.

Conclusion
Phylogenetic analysis of the late Permian to Triassic rhynchonellide brachiopods was performed with four methods: EW, IW, UB, and TB. The results indicate that analytical methods have a significant effect on the topology of trees and any following analysis based on these trees. Compared with trees generated by EW, IW, and UB in which Permian and Triassic taxa are mixed together, tip-dated trees appear to be more reasonable and the revealed evolutionary relationships more consistent with the fossil record. Although the topology and branch lengths vary greatly among trees generated by different analytical approaches (species-dating vs. genus-dating, tip-dating vs. posterior rescaling), the downstream analyses based on multiple posterior trees generate some comparable results. According to the cal3calibrated trees, the major increase in lineage richness occurred in the Early and early Middle Triassic and reached its highest level in the Anisian. The Anisian taxa evolved complex and diverse internal and external characters, implying the full recovery of this order. Both shell size and the strength of ornamentation diminished rapidly after the P/Tr extinction. The decline in these two measures was likely caused by the disappearance of larger and sculptured Permian genera and the radiation of minute and weakly ornamented taxa. Ancestral-state estimations of shell size and the development of ornamentation, along with comparisons of other characters, show that the Early to Middle Triassic small-sized taxa probably have characters displayed by juveniles of their ancestors, implying that paedomorphosis was a likely survival strategy that developed in the adverse environments after the P/Tr extinction.
In this study, we provide a preliminary tipdated analysis for the late Permian and Triassic rhynchonellides using the relatively simple fossilized birth-death model. However, further efforts based on more complicated models and informative priors will be a key future line of enquiry. In addition, further paleontological studies on the Early Triassic rhynchonellides are needed to develop more persuasive conclusions for the paedomorphosis of the small-sized taxa. Directly sampling their ancestors and studying their ontogenies will provide tests for these hypotheses.

Acknowledgments
We sincerely thank A. M. Popov, B. Radulović, M. Siblík, and X. Guo for their help procuring some publications. We are very grateful to S. Arias for sharing his TNT script. T. L. Stubbs and M. J. Benton are acknowledged for helpful guidance on the analytical approaches. We thank M. Clapham, three anonymous reviewers, and associate editor J. Crampton for their critical comments and constructive suggestions that have improved greatly the quality of this paper. The contributors of the Paleobiology Database are also thanked for their provision of the taxonomic data for the database. This study was supported by three NSFC grants (41821001, 41930322, and 41772007), 111 Project of China (grant no. BP0820004), and the Fundamental Research Funds for National Universities, China University of Geosciences (Wuhan). This is a contribution to the IGCP 630: "Permian-Triassic Climatic and Environmental Extremes and Biotic Responses."