Selection towards larger body size in both herbivorous and carnivorous synapsids during the Carboniferous

Body mass was estimated using a method suggested by Romer and Price (1940). The cross-sectional area of the dorsal vertebrae (estimated from the radius of the dorsal vertebral centrum squared) was used as a proxy for mass. Mass estimates of Sphenacomorpha were taken from the original compendium, and were calculated for those taxa not included in these publications but that possessed dorsal vertebrae. Body masses were log-transformed prior to analysis. The data are presented in Supplementary Material 1.

It should be noted that a more recent method of estimating body mass is available: using the shaft circumference of the femur and the humerus (Campione and Evans 2012). However, the use of this method would result in a smaller sample size, as there are fewer sphenacomorph specimens with both a femur and humerus preserved than with dorsal vertebrae preserved, and the largest specimens of several taxa do not preserve these bones.

The phylogeny used in this study is a composite tree formed from the phylogenies of Mazierski and Reisz (2010) and an expanded version of Brink et al. (2015). The former is the most comprehensive phylogenetic analysis of Edaphosauridae, and the latter is the most comprehensive phylogenetic analysis of Sphenacodontidae. To expand the analysis of Brink et al. (2015), the taxa Ctenorhachis jacksoni and Cryptovenator hirschbergeri were added, along with non-sphenacodontid outgroups (Varanops, Ophiacodon, and Edaphosaurus). The matrix used in this analysis is presented in Supplementary Material 2. In the edaphosaurid phylogeny, a polytomy is present due to the lack of overlap between Glaucosaurus and Lupeosaurus. This polytomy was resolved according to the Bayesian analysis of Brocklehurst et al. (2016), where Lupeosaurus was found to be more closely related to Edaphosaurus than Glaucosaurus. The two resulting phylogenies were combined in such a way that Sphenacodontia was the sister to Edaphosauridae, and non-sphenacomorph outgroups were deleted. The resulting tree is shown in Fig. 1.

The phylogeny was time-calibrated using the method of Brusatte et al. (2008) in the R 3.14 package (R Core Team 2014) paleotree (Bapst 2012). Under this method, zero-length branches resulting from inconsistencies between the order of branching observed in the phylogeny and the order of appearances of the taxa in the fossil record were treated by sharing the length of the branch immediately ancestral to them equally between ancestor and descendants. Uncertainty surrounding the ages of taxa was accounted for using the method of Pol and Norrel (2006). For each taxon, 100 first appearances and last appearances were drawn at random from a uniform probability distribution covering the full possible range of ages for that taxon. One hundred time-calibrated trees were produced from the 100 sets of ages, and all phylogenetic comparative analyses were performed on all 100 of these trees. The age ranges used are presented in Supplementary Material 1.

Taxa for which body mass estimates could not be calculated due to the lack of dorsal vertebrae were retained in the phylogeny during time calibration, as they provide important information on branch lengths and the ages of nodes; however, they were dropped after time calibration. These include Cryptovenator, Dimetrodon borealis, Tetraceratops, and Glaucosaurus. The therapsid lineage was also pruned, as no information is available on the body size evolution of this clade during the Carboniferous and the early Permian. The 100 trees subjected to analysis after pruning are presented in Supplementary Material 3.

Another possible source of uncertainty is the existence of taxa not included in the phylogeny. The analysis of Mazierski and Reisz (2010) included all valid edaphosaurid taxa with the exception of Edaphosaurus mirabilis. However, because this taxon is currently the smallest known member of the genus Edaphosaurus, its absence could potentially have a large impact on the results. The analysis of Brink et al. (2015), though comprehensive in terms of genera, only included four out of the 11 species of Dimetrodon considered valid in the most recent summary of early synapsids (Brocklehurst et al. 2013). There are also species of Haptodus and Sphenacodon not included. To test the impact of these exclusions, another set of 100 time-calibrated trees were produced, in which all missing species for which body masses could be estimated were added at random to the tree, within the genera to which they had been assigned. These 100 trees are presented in Supplementary Material 4. The analyses of phenotypic selection (see below) were performed again on these 100 trees.

Morphological evolution is not homogeneous. Different selective pressures acting on different lineages will lead to deviations from a “Brownian motion” (random walk) model of evolution. Recently, Baker et al. (2015) proposed a method of identifying lineages along which “positive phenotypic selection” occurred. The concept of positive selection was defined as change beyond what would be expected under Brownian motion. The difference between the morphological change occurring along a lineage and the change that might be expected was calculated using the variable rates model of Venditti et al. (2011). This model has Brownian motion as its basis, but scales branch lengths within the phylogeny to represent deviations from this model. A Markov Chain Monte Carlo reversible-jump framework was used to generate posterior distributions of scalars for each branch and to deduce what pattern of rate variation best fits the morphological data observed. Baker et al. (2015) proposed (somewhat arbitrarily) that positive phenotypic selection could be deemed to have occurred along a branch if the amount of morphological change occurring along that branch was more than double that which would be expected under Brownian motion.

All the time-calibrated trees described above were analysed using the method proposed by Baker et al. (2015). The analyses were performed in BayesTraits V2.0. Markov chains were run for 10 000 000 iterations, with 100 000 discarded as burn-in.

To ascertain the direction in which the morphology was being selected, ancestral body masses were calculated using likelihood with the ace function in the package ape in R (Paradis et al. 2004). This function assumes a Brownian motion model of evolution, but here the branch lengths of the phylogeny were rescaled as described above to represent the variation in rate.

A likelihood approach was used to identify points in the sphenacomorph phylogeny where heritable shifts in the rate of body size evolution occurred. The 100 time-calibrated phylogenies were analysed using the trait medusa algorithm (TM1) in the R package models of trait macroevolution on trees (MOTMOT) (Thomas and Freckleton 2012). This algorithm uses a stepwise Akaike information criterion (AIC) approach to identify the combination of rate shifts that best fit the observed data. First, the likelihood of a single-rate Brownian motion model was calculated. Then, a two-rate model was evaluated, with the likelihood of a rate shift at each node calculated in turn. Models with more than two rate shifts were evaluated with the location of one of the shifts constrained at the node identified in the previous step. This process continued until a user-specified number of shifts had been identified (in this case five). AIC and Akaike weights were used to assess which pattern of rate shifts best fits the body mass data and phylogeny.

While phylogenetic comparative methods allow one to incorporate the relationships of taxa in analyses of evolution, they do have some limitations. First, one is limited by the taxa that have been included in the phylogeny. Also, at present, the model-fitting analyses do not allow one to incorporate information on changes in body size within a single tip taxon. To examine whether incorporating such factors might not alter the results, body mass data were analysed as a time series rather than within a phylogeny. This allowed one to use all species for which a body mass could be calculated. It also allowed the inclusion of different specimens of the same taxon; the largest specimen of each taxon in each time bin was used, thus allowing the same taxon to have different body masses at different points in time. The specimens included in this dataset, their ages, and estimated masses are presented in Supplementary Material 5. In each time bin, the sample size, trait median, and standard error were calculated. Four evolutionary models were fit to the time series. The first, an unbiased random walk (URW) model, represents evolution by Brownian motion; one expects random trait variation around a constant mean, with the trait variance increasing through time (Hunt 2006). The generalised random walk (GRW) model also represents a Brownian motion process, but the trait mean shows a trend in a particular direction. The stasis model represents the trait being held to an adaptive optimum through time by a “rubber band” parameter; the further the trait deviates from the optimum, the more strongly it is drawn to it (Hunt 2006). Thus, the trait mean will remain constant, as will the variance. Finally, a model of punctuated evolution was tested: a stasis model with a shift in the adaptive optimum between two time bins. Model-fitting was performed in the R package paleoTS (Hunt et al. 2008), and models were selected based on the Akaike weights.

Two lineages exhibited positive phenotypic selection in all 100 of the time-calibrated phylogenies (Fig. 2). The first of these lineages was within Edaphosauridae: the lineage leading to the clade containing Lupeosaurus and Edaphosaurus. The second was the lineage leading to Sphenacodontidae. When one examines the inferred ancestral body masses (Fig. 3), it can be observed that both of these lineages demonstrate a rapid increase in body size. All other lineages were found to experience phenotypic selection in either none or an extreme minority of the 100 trees.

The addition of taxa not yet tested in phylogenetic analysis to the tree did not change the inference of positive phenotypic selection within Edaphosauridae; the same lineage was found to have experienced significant selection in all trees. Edaphosaurus mirabilis, one of the taxa added to random locations within Edaphosaurus and the smallest included member of this genus, was also found to have experienced significant phenotypic selection in most of these trees. However, the inference of positive selection towards larger body size in the lineage leading to Sphenacodontidae was less certain when the new taxa were added; significant selection along this lineage was found in only 42 of the 100 phylogenies.

In 91 of the 100 phylogenies, the TM1 analysis identified a single rate shift as best fitting in the observed data. This shift is located at the node ancestral to the two youngest and largest species of the genus Edaphosaurus: E. cruciger and E. pogonias (Fig. 2). The shift is towards lower rates of body size evolution in this clade. The results of the 100 TM1 analyses are presented in Supplementary Material 6.

According to the Akaike weights, the time series of body masses best fits the punctuated model (Table 1), with an increase in trait optimum occurring at the end of the Gzhelian (Fig. 4). The support for this model was not overwhelming, however, and the URW model received an almost identical Akaike weight. The lack of overwhelming support for any particular model is likely due to the short time series.

An increase in body size in Permo-Carboniferous non-mammalian synapsids has been qualitatively noted since specimens were first placed in a stratigraphic context (Romer and Price 1940; Olson 1962; Gould 1967; Pivorunas 1970; Gould and Littlejohn 1973; Bakker 1986; Tracy et al. 1986; Berman et al. 2001; Tomkins et al. 2010; Reisz and Fröbisch 2014; Brink et al. 2014, 2015). Reisz and Fröbisch (2014) proposed that a trend towards increased body size can be observed in all herbivorous lineages, and also posited a body size increase in the Sphenacodontidae; although they suggested that this increase was less than that observed in herbivores. The results of this study support the hypothesis of Reisz and Fröbisch (2014); a body size increase was observed in edaphosaurids and sphenacodontids. However, this increase was not a gradual trend occurring throughout the Carboniferous and early Permian. Only two sphenacomorph lineages consistently showed phenotypic selection towards larger body size regardless of the uncertainty surrounding the ages of taxa. All other lineages exhibited changes in body mass not significantly different from what would be expected from Brownian motion.

Both of these short, rapid increases in body size occurred during the latest Carboniferous, during the Kazimovian and Gzhelian stage. This observation is supported by the time series of body masses, which showed an increase in mass at this time, followed by a period of relative stability (Fig. 3). The fitting of evolutionary models to the time series suggests that the punctuated model, with stasis during the Carboniferous and early Permian and a shift to a larger body size in between, best fits the data (although one should again emphasise the weakness of the support for this model). The timing of the shift in body mass is interesting, as it contradicts a suggestion made previously on body mass evolution in pelycosaurian-grade synapsids. In an analysis of changes in species richness of early synapsids (Brocklehurst et al. 2013), it was tentatively suggested that a rapid trend towards warming and drying during the early Permian at the end of the Sakmarian (Montanez et al. 2007), which reduced the extent of the equatorial forests, may have promoted the evolution of larger body sizes; a more open environment would allow larger inhabitants (Brocklehurst et al. 2013). The results here show that this is unlikely to have been the case; although many of the largest pelycosaurian-grade synapsids do not appear in the fossil record until after this event, the initial evolution of large body size appears to have occurred much earlier, during the late Carboniferous. The transition to a more open habitat may have allowed larger taxa to become more abundant, and thus more likely to be preserved (Brocklehurst and Fröbisch 2014), but it did not precipitate the increase in body size.

The evolution of the dorsal sail showed differing patterns in edaphosaurids and sphenacodontids. The dorsal sail evolved independently in each clade, and evolved more than once in sphenacodontids (once in Secodontosaurus and once at the base of the Dimetrodon clade). Sphenacodontids increased in body size before the first occurrence of the dorsal sail in Secodontosaurus, and species within Sphenacodontidae that do not have extremely elongate neural spines (Sphenacodon Ctenorhachis and Ctenospondylus) still evolved large body sizes, comparable to contemporaneous species of Dimetrodon. This pattern suggests no correlation between the evolution of large body sizes and the development of a dorsal sail in sphenacodontids. However, edaphosaurids increased in body size after the first occurrence of the dorsal sail in Ianthasaurus, and so it is possible that the two traits are related. Given the differences in patterns between body size and dorsal sail evolution, the differences in neural spine shape and ornamentation, and that thermoregulatory hypotheses for the dorsal sail have been questioned, these results lend support to the hypothesis that elongate neural spines are a secondary sexual character in Sphenacomorpha (Tomkins et al. 2010). The dorsal sail, along with increased body mass, became more prominent after the decline of rainforests at the end of the Carboniferous (Sahney et al. 2010), where a dorsal sail would have been less cumbersome.

Selection towards larger body size in edaphosaurids appears to coincide with the evolution of herbivory. Although no skull material from Lupeosaurus has been preserved, the ribs show curvature throughout their length (Sumida 1989). This produces the barrel-like chest cavity characteristic of herbivorous pelycosaurian-grade synapsids, thought to be necessitated by the long intestine required for gut fermentation of plant material (Sumida 1989; Reisz and Fröbisch 2014). Edaphosaurus’s status as an herbivore is clear; the leaf-shaped serrated teeth and occluding tooth plates formed from the palatal and dentary teeth are indicative of an herbivorous diet (Romer and Price 1940; Modesto 1995; Sues and Reisz 1998).

The coincidence of the evolution of herbivory and selection towards large body size in edaphosaurids provides support for the observation of Reisz and Fröbisch (2014) that the earliest terrestrial herbivores showed a tendency towards larger body size, although it should be noted that the universality of this “rule” has been challenged; in captorhinids, the evolution of large size was found to be unrelated to diet (Brocklehurst 2016). One must also consider the precise reason behind this relationship. In African ungulates, larger body size is correlated with reduced predation; larger taxa have fewer predators and a lower percentage of their mortalities are caused by predation (Sinclair et al. 2003). However, the increase in body size in the predatory sphenacodontids, the largest and most abundant terrestrial predators of their time, occurred later than the increase in herbivorous edaphosaurids (Fig. 3). Therefore, increased predation seems an unlikely selective pressure towards larger size. Alternative explanations are physiological. The abundance-packet size hypothesis (Olsen 2015) suggests that because larger taxa have greater energy requirements, a lineage that evolves a larger body size will experience a selective pressure towards finding a more abundant food source, such as plant matter (or alternatively find food in larger “packets”, i.e., feed on larger prey). However, this explanation does not fit in the case of edaphosaurids. The abundance-packet size hypothesis suggests that the transition to herbivory should occur in larger taxa; that is, the diet transition should follow the change in body size. This hypothesis provides no selective pressure towards larger body size in herbivores, and cannot explain the phenotypic selection identified in this lineage.

A more fitting model linking herbivory and body size is the Jarman-Bell principle (Geist 1974). This principle posits that although larger animals have larger absolute energy requirements, smaller animals have higher metabolic energy requirements relative to their body size. Therefore, larger herbivores are able to subsist on lower quality plant material such as leaves. Because such low quality plant material is more abundant than higher quality material such as roots and fruits, this model does supply a selective pressure towards larger body size in an herbivorous lineage, and can therefore provide an explanation for the positive phenotypic selection observed in edaphosaurids.

One must consider the small specimen from the Stephanian-aged Rakonitz coal basin (equivalent to Gzhelian) of the Czech Republic, assigned to the species Edaphosaurus mirabilis. If its assignment to Edaphosaurus is correct, this would represent the smallest species of this genus, with an estimated mass of 5 kg (Romer and Price 1940). The assignment of this specimen to Edaphosaurus is unclear due to its incompleteness (only a single dorsal centrum and partial neural spine is preserved), and it has been suggested that it may belong to Ianthasaurus (Modesto and Reisz 1990a), a small insectivorous edaphosaurid from the Carboniferous of the USA (Reisz and Berman 1986; Modesto and Reisz 1990b; Mazierski and Reisz 2010). Nevertheless, even if its assignment to the herbivorous Edaphosaurus is assumed to be correct, this does not alter the inference of phenotypic selection; when E. mirabilis is added 100 times to random positions within Edaphosaurus, positive phenotypic selection towards larger body size is still identified along the same lineage leading to the herbivorous edaphosaurids. Significant phenotypic selection is also identified in E. mirabilis itself in the majority of these 100 phylogenies, indicating that if its assignment to Edaphosaurus is correct and it is not a juvenile, it represents an example of miniaturisation. It does not overturn the inference of selection towards larger body size in the herbivores.

The known Carboniferous and early Permian Sphenacodontia form an entirely carnivorous set of taxa. While therapsids are included in the clade Sphenacodontia, and do include herbivorous taxa, they are unknown from the Carboniferous and early Permian fossil record with the exception of a set of sacral vertebrae (Spindler 2014). As such, therapsid body size evolution cannot be compared with the Sphenacodontidae or Edaphosauridae, and is discounted from this discussion.

As was observed in edaphosaurids, Sphenacodontia showed significant selection towards a larger body size occurring along a single lineage leading to the Sphenacodontidae, the dominant clade of large terrestrial marco-carnivores in the early Permian. This rapid increase in body size occurred shortly after it occurred in Edaphosauridae, and the two may have been linked. It is unsurprising that larger carnivores have been shown to be able to feed on a wider range of prey (Sinclair et al. 2003). As the appearance of large Edaphosauridae provided a wider range of available prey sizes, there would have been a selective pressure towards larger carnivores able to exploit the greater variety of prey items available.

It should be noted, however, that this explanation does rely on the inference of a predator–prey relationship in the absence of definite evidence. While Dimetrodon and Edaphosaurus are often assumed to have existed in a predator–prey relationship (Romer and Price 1940), there is little actual data for the diet of sphenacodontids in the fossil record (although one notable example is a skeleton of the temnospondyl Zatrachys in the body cavity of Dimetrodon milleri (Romer and Price 1940)). There is a definite disparity in the abundance of sphenacomorphs, with edaphosaurids being considerably less common than their supposed sphenacodontid predators (Olson 1966). This latter point implies that sphenacodontids must have, at the very least, supplemented their diet with other tetrapods or fish (Romer and Price 1940; Olson 1966; Bakker 2015).

Moreover, the nature of the body size increase in Sphenacodontia is less clear than in Edaphosauridae. While positive phentotypic selection is well supported in all 100 time-calibrated trees containing the taxa subjected to cladistic analysis, when the sphenacodontian taxa not subjected to phylogenetic analysis were added at random to their supposed genera, support for phenotypic selection was only found in 42% of analyses. It appears that the inference of significant selection towards larger body size at the base of Sphenacodontidae is highly dependent on the taxa included in the analysis.

Even if we accept the results of the analyses including the smaller sample of taxa, the magnitude of the increase in size at the base of the Sphenacodontidae is less than that of Edaphosauridae. A further difference between the two can be observed in the period following these two rapid increases. While the carnivorous Sphenacodontidae continued to increase in size throughout the early Permian, the herbivorous members of Edaphosauridae showed little change in body mass following their initial increase (Fig. 3). The results of the TM1 analysis support this observation; a shift towards lower rates of size evolution was identified at the node containing the largest species of Edaphosaurus (Fig. 2), which presumably have converged on this adaptive optimum.

This last contrast allows us to infer differing modes of phenotypic selection acting on the high-fibre herbivores and the macro-carnivores. While the lineage leading to herbivores does exhibit selection towards a mass larger than that of the ancestral taxa, this appears to be selection towards an adaptive optimum rather than a constant trend. Once this optimum is reached, there is a shift towards lower rates of evolution and body sizes remain within a relatively narrow range. On the other hand, the sphenacodontid top-predators continued to evolve ever-larger body sizes throughout the early Permian. By the Kungurian, the lineage leading to Dimetrodon grandis had exceeded the largest known edaphosaurid. The increased diversity and abundance of some of the largest terrestrial animals of the early Permian, the caseids, in the latter stages of the Cisuralian (Olson 1954, 1968; Reisz et al. 2011; Brocklehurst et al. 2013; Romano and Nicosia 2014) may provide a reason for the continued increase in body size: the need to expand the possible diet range to include these even larger prey items. The sphenacodontids also evolved a more sophisticated dental apparatus at this time: “true” ziphodonty, where the serrations of the teeth have denticles with a dentine core (Brink and Reisz 2014).