The impact of paleoclimatic changes on body size evolution in marine fishes.

Body size is an important species trait, correlating with life span, fecundity, and other ecological factors. Over Earth's geological history, climate shifts have occurred, potentially shaping body size evolution in many clades. General rules attempting to summarize body size evolution include Bergmann's rule, which states that species reach larger sizes in cooler environments and smaller sizes in warmer environments, and Cope's rule, which poses that lineages tend to increase in size over evolutionary time. Tetraodontiform fishes (including pufferfishes, boxfishes, and ocean sunfishes) provide an extraordinary clade to test these rules in ectotherms owing to their exemplary fossil record and the great disparity in body size observed among extant and fossil species. We examined Bergmann's and Cope's rules in this group by combining phylogenomic data (1,103 exon loci from 185 extant species) with 210 anatomical characters coded from both fossil and extant species. We aggregated data layers on paleoclimate and body size from the species examined, and inferred a set of time-calibrated phylogenies using tip-dating approaches for downstream comparative analyses of body size evolution by implementing models that incorporate paleoclimatic information. We found strong support for a temperature-driven model in which increasing body size over time is correlated with decreasing oceanic temperatures. On average, extant tetraodontiforms are two to three times larger than their fossil counterparts, which otherwise evolved during periods of warmer ocean temperatures. These results provide strong support for both Bergmann's and Cope's rules, trends that are less studied in marine fishes compared to terrestrial vertebrates and marine invertebrates.

Paleoclimatic changes are recognized as strong factors affecting the macroevolutionary dynamics of clades, including their distribution, ecology, and diversification (1). Throughout the course of Earth’s geological history, several large, dynamic climatic shifts have occurred, such as the end-Permian extinction event (ca. 252 Ma), the Cretaceous–Paleogene extinction event (ca. 66 Ma), and the Paleocene-Eocene Thermal Maximum (ca. 55.6 Ma) (2, 3). These periods are often marked by large changes in temperature, ocean acidification, and anoxia, as well as increases in volcanic activity (3). These environmental shifts have led to mass extinction events in fishes (4) and changes in rates and magnitude of body size evolution in amphibians, birds, and mammals (5, 6), among others. Morphological responses to paleoclimate change can be directly observed from the fossil record. As body size correlates with many aspects of a species’ biology, physiology, and ecology, its evolution should be associated with shifts in climate (5, 7, 8).

Bergmann’s rule attempts to summarize body size responses to climatic changes, stating that species within a clade (or populations within a species) tend to grow to larger sizes in cooler environments and smaller sizes in warmer environments. While Bergmann’s rule can apply at multiple evolutionary scales (9), from an interspecific viewpoint it can be defined as an ecogeographical trend where species’ body size varies as a negative function of temperature. Originally studied in mammals, this trend has now been identified in a range of animals such as crustaceans, amphibians, and ray-finned fishes (8, 10–12). Various explanations have been proposed for Bergmann’s rule, from heat conservation in endotherms to oxygen availability in ectotherms (13, 14). Another broad hypothesis summarizing body size patterns is Cope’s rule, stating that species tend to increase in size over evolutionary time. Explanations for Cope’s rule are thought to be linked to fitness advantages at larger body sizes or an increase in size variance as lineages diversify from a smaller ancestor following a passive trend (15, 16). Cope’s rule could simply be an evolutionary or temporal manifestation of Bergmann’s rule if lineages evolve larger body sizes during periods of climatic cooling (8). This idea, termed the Cope–Bergmann hypothesis by Hunt and Roy (8), has received considerably less attention than studies that examine body size trends relating to Cope’s and Bergmann’s rules separately (but see 8, 17).

Species’ responses to climate will vary, but ectotherms that rely on their environment for temperature regulation are likely more susceptible than endotherms to climatic changes (18). Temperature controls a variety of aspects of ectotherm biology and is strongly linked with an organism’s fitness, affecting growth rates and overall body size (14, 18). Understanding ectotherm morphological responses to global paleoclimate change may benefit greatly from examining clades with a rich fossil record. Ectothermic invertebrates have been studied in great detail, particularly brachiopods and marine arthropods (8, 19–21), due to their exceptional fossil record (22). Among ectothermic vertebrates, Cope’s and Bergmann’s rules have been tested in amphibians and reptiles (5, 23), but comparatively less in teleost fishes (but see 7, 24).

Fishes in the order Tetraodontiformes provide a model clade to test patterns of body size evolution in relation to paleoclimate events, owing to their exemplary and well-studied fossil record and extraordinary morphological diversity (4, 25). They constitute a circumglobally distributed taxonomic order of mostly marine, subtropical/tropical dwelling fishes, represented by ca. 450 living species, including the charismatic pufferfishes, triggerfishes, and ocean sunfishes. Tetraodontiforms exhibit a diverse array of body shapes, from nearly square (boxfishes) to globose (pufferfishes) and laterally compressed (filefishes). Species in this order also feature remarkable variation in adult body size, ranging from just 25 mm total length (TL) (e.g., Rudarius excelsus, Carinotetraodon salivator) to 3.4 m TL (e.g., Mola, Masturus lanceolatus). The tetraodontiform fossil record extends to the Late Cretaceous with representatives from 12 exclusively fossil families; all 10 extant families are also present in the paleontological record, and on average, body size is smaller among fossil taxa (25). Their morphological diversity, coupled with a robust fossil record, provides a unique system to test the Cope–Bergmann hypothesis in ectotherms.

This study aimed to investigate patterns of body size evolution in relation to paleoclimate events in ectotherms using tetraodontiform fishes as a model clade. We addressed the following questions: 1) Are paleoclimatic changes correlated with changes in tetraodontiform body size? and 2) If a correlation between paleoclimate and body size is observed, do tetraodontiforms follow the Cope–Bergmann hypothesis? That is, does their body size evolution correlate with a paleotemperature curve where tetraodontiforms are evolving toward a larger body size during periods of climate cooling? To address these questions, we estimated a time-calibrated phylogeny for tetraodontiforms using total-evidence dating approaches that combine genome-wide data from extant species with a morphological matrix coded from both fossil and extant species. We also incorporated body length data and paleotemperature records spanning the past 100 Ma into a series of evolutionary model fitting analyses. We hypothesize that tetraodontiform evolution has been driven by past temperature changes and that body size is strongly linked to past climate.

Results

Phylogenomic Inference, Total-Evidence Dating, and Tree Uncertainty.

Under a total-evidence framework that combines a phylogenomic dataset based on 1,103 nuclear markers and 210 morphological characters (4), we used Bayesian inference and divergence time analyses to generate the most complete tetraodontiform phylogeny to date (Fig. 1). Our approach used extensive taxon sampling that included both newly generated sequences (Dataset S1) and previously published sequences (Dataset S2) (26) for a total of 185 extant (out of ca. 450) and 52 fossil taxa. In addition, we tested other phylogenetic inference methods such as concatenation-based maximum likelihood (ML) and summary multispecies coalescent approaches (). Phylogenetic placements are remarkably stable and congruent with past molecular studies for the group (27, 28) conducted before the advent of genomic datasets, supporting the monophyly of all families and the seven main suborders. Results show some differences between the major approaches conducted, with placements for some lineages (e.g., Triacanthidae + Triacanthodidae) along the backbone varying due to the short internodes at the base of the trees ().

Fig. 1.

Tip-dating tree inferred for Tetraodontiformes and evolutionary model fitting results, including superfamily Plectocretacicoidea. (A) The MCC tree derived from a total-evidence dating analysis using the FBD model in MrBayes including plectocretacicoids (see for a tree excluding plectocretacicoids). The MCC tree is derived from 10,000 trees evenly sampled from the posterior distribution of five independent genomic subsets. provides an expanded version of this tree. (B) Raincloud plots (half-violin plots and boxplots) for each model of body size evolution tested, representing the distribution of likelihood scores from 500 trees evenly selected from the posterior distribution of five independent subsets in the Bayesian analysis (see also Table 1). Dots represent the raw likelihood score for each of the 500 trees analyzed for each model. Evolutionary models include EB, BM, rate trend, OU, mean trend, a climate OU model using tropical ocean temperatures, and a climate OU model using GAT.

Table 1.

Evolutionary model fitting results, including the superfamily Plectocretacicoidea

Evolutionary model	Parameters	AICc		lnL	AICw
OU climate GAT	5		98.25	−43.99	0.999
OU climate tropical	5		112.42	−51.08	8.35e-04
Mean trend	3		126.91	−60.40	3.60e-06
OU	3		137.88	−65.89	2.47e-09
Rate trend	3		146.78	−70.34	2.89e-11
BM	2		148.01	−71.98	1.56e-11
EB	3		148.03	−70.96	1.54e-11

Model fitting results for the seven macroevolutionary models tested on 500 trees selected from the combined posterior distributions of five genomic subsets. Number of model parameters, mean values for the corrected Akaike information criterion (AICc), mean log likelihood (lnL), and weighted AICw are reported. The strongest support went to the climate-driven OU model using the global average sea surface temperature curve.

To account for phylogenetic and divergence time uncertainty on downstream comparative analyses, we analyzed five independent and largely nonoverlapping genomic subsets, each with a sufficient number of genes to reduce sampling error (29, 30). This approach contrasts with the traditional way of analyzing trees obtained from a Bayesian posterior distribution, which is typically estimated using a concatenated alignment with a scant number of genes. We also used two separate fossil schemes by either including or excluding the superfamily Plectocretacicoidea, whose placement within Tetraodontiformes has been controversial (29, 31). Given the Mesozoic origin of plectocretacicoid fossils (70 to 96.9 Ma), their inclusion/exclusion in the dataset has important implications for tetraodontiform ages (29, 32–34). Our divergence time estimates, including plectocretacicoids, place crown Tetraodontiformes within the Late Cretaceous (92.21 Ma, 95% highest posterior density [HPD]: 86.78 to 113.16 Ma) and the stem age at 98.62 Ma (95% HPD: 96.07 to 114.42 Ma; Fig. 1 and ). Excluding the plectocretacicoids from the analysis shifts the age of crown tetraodontiforms forward to 62.45 Ma (95% HPD: 60.52 to 87.30 Ma; ), which is consistent with other estimations that exclude this extinct superfamily (32). We also assessed the sensitivity of divergence time estimates to root prior choice, finding no strong effects ().

Evolutionary model fitting.

To test the Cope–Bergmann hypothesis in tetraodontiforms, we fit a series of models of continuous trait evolution using mean maximum standard length as a proxy for body size (see ). To account for uncertainty in tree topology and divergence time estimates, each model was tested on 500 trees evenly selected from the posterior distribution of five independent subsets that reached convergence in the Bayesian Markov Chain Monte Carlo runs using the two alternative fossil schemes (i.e., including and excluding Plectocretacicoidea). Models tested included the simple Brownian motion (BM; random walk) model, an early-burst (EB) model, an Ornstein-Uhlenbeck (OU) model, two variants of the trend model (mean trend and rate trend), and a pair of climate-dependent OU models, each fitted using a separate temperature curve as input. While the trend models tested are assumed to explicitly model Cope’s rule, where time is the sole factor responsible for an increase in clade’s body size, the OU climate model allows for tests of both Cope’s and Bergmann’s rules, where both time and temperature can influence the underlying trait (see ). Paleotemperature curves used for the OU climate model include sea surface temperatures at tropical latitudes (15°N to 15°S) and the global average temperature (GAT) for the past 100 Ma (2). Given the broad circumtropical/subtropical distributions and habitat preferences of tetraodontiforms (mainly marine shallow-water dwellers), these two contrasting temperature curves most accurately capture the spectrum of past environmental affinities in this group. Nonetheless, we also tested an additional curve based on deep-sea temperature data (35), which yielded similar results ().

Highest support went to the climate-driven model of evolution using a global average sea surface temperature curve (Fig. 1), with an Akaike weight (AICw) of 0.999 (Table 1). All other models received substantially less support. These results are robust to the inclusion or exclusion of the superfamily Plectocretacicoidea (Fig. 1 and Table 1 and ). The climate OU GAT model resulted in the best fitting model for 329/500 (65.8%) trees analyzed (). To visualize the scale and directionality of body size over time, we reconstructed and plotted ancestral body size as a traitgram. On average, extinct species and families are two to three times smaller in size than extant species and families (Fig. 2). When plectocretacicoids were excluded, similar patterns were observed ().

Fig. 2.

Evolution of tetraodontiform body size over time, including superfamily Plectocretacicoidea. Ancestral reconstruction of body size in tetraodontiforms, as estimated using the contMap function in the R package phytools (36). The log-transformed mean maximum SL for each species is plotted as a traitgram on the y axis, with time on the x axis. Fish silhouettes are scaled to represent proportional log body size and are colored by family, with extinct families in gray. The estimated ancestral body size of tetraodontiforms is two to three times smaller than the mean of present-day taxa.

In addition, we examined the potential effect of tree age on model selection in our subsets inclusive of Plectocretacicoidea. Our 500 trees ranged in age from 96.25 to 130.86 Ma, with an average root age of 103.18 Ma (). It appears that slightly younger trees were favored by the Mean Trend model (mean age = 99.44 Ma; ) compared with the two climate OU models (GAT mean age = 103.47 Ma, Tropical mean age = 104.49 Ma). However, this could be an effect of the fewer number of trees being favored by the Mean Trend model overall (n = 64) compared with the GAT (n = 329) and Tropical (n = 107) curve climate OU models. While model selection on a small sample of trees can produce biased estimates, using a large number of trees obtained from the posterior distributions of independently assembled gene subsets provides a powerful approach to account for tree uncertainty in macroevolutionary inferences (e.g., 29, 30).

Because there is a global trend toward declining temperatures over the evolutionary history of Tetraodontiformes (i.e., from the Late Cretaceous to the present day), it can be difficult to decouple the effects of temperature (Bergmann’s rule) from other processes that may be correlated with increased body sizes (Cope’s rule), as the patterns generated by either rule can be indistinguishable from one another. To further assess the role of the overall trend, we decomposed the temperature curve into two distinct components: the smoothed overall trend and the fluctuations around this trend (). To identify the model parameters with the strongest weights, we ran three separate analyses on both temperature curves (GAT and tropical latitudes). The first analysis modeled the two independent parameters (overall trend + fluctuations) together, the second modeled only the overall trend, and the third modeled only the fluctuations. If the model with only the overall trend showed the best fit for the data, this would imply that the trend is more important than the temperature fluctuations (i.e., some evidence for Cope’s rule, but inconclusive for Bergmann’s). In our analyses using the GAT curve, we found that most support went to the overall trend + fluctuations model (AICw = 0.769), followed by overall trend (AICw = 0.152) and fluctuations (AICw = 0.077; ) models. Although we observed a different pattern in the tropical latitude curve, with most support going to the fluctuations model (AICw = 0.852 versus AICw = 0.147 for trend + fluctuations model, and AICw = 2.77e-06 for trend model; ), given that the original model fitting results indicated stronger support for a GAT curve over a tropical latitude curve (Table 1), we place more weight onto the decomposed GAT curve model results. All in all, our decomposed model analyses suggest that both the overall trend and the fluctuations around this trend are important to the OU climate model fit, providing support for the Cope–Bergmann rule (8) as an explanation of body-size evolution in tetraodontiforms.

Ecomorphological correlations.

To test whether past ocean temperatures are correlated with tetraodontiform body size, we performed phylogenetic generalized least squares (PGLS) regression analyses under a best-fit model between reconstructed ancestral node body sizes and paleo-ocean temperatures at the age of each node using the maximum clade credibility (MCC) tree ( provides details). While phylogenetically informed statistical tests (such as PGLS) are traditionally used to compare two continuous species traits, we argue that for the purposes of our analyses, ocean temperature can be categorized as a species trait, following the approach of Garland et al. (37), where environmental traits can be used as long as these traits can be passed on from ancestor to descendent species. While many studies use a temperature curve based on deep-sea data for these types of analyses (6, 38), we accounted for temperature variation relating to the actual habitat and ecology of tetraodontiforms by testing two temperature curves based on sea surface temperatures, as described above (2). From the Late Cretaceous (∼100 Ma) to the present day, sea surface temperatures (both at tropical latitudes and the global average) have been gradually decreasing, a trend that correlates negatively with tetraodontiform body size (Fig. 3). These PGLS regressions under a best-fit OU model are statistically significant for the GAT (P = 5.571e-03), but not the tropical latitude temperature (P = 0.0653) curves (). These results reflect the model fitting analyses, where the GAT curve received substantially more model support than the tropical temperature curve (Table 1). Even when the superfamily Plectocretacicoidea was excluded from the analysis, shifting tetraodontiform divergence time estimations forward ∼30 Ma, the PGLS analysis remained statistically significant for the GAT curve (P = 4.84e-05) but not the tropical curve (P = 0.447; ). Examining trends among only fossil species revealed a pattern similar to the analyses where fossil and extant species were combined (). Fossils-only analyses showed that body size was strongly correlated with ocean temperatures for the global average curve (P = 0.0311) but not the tropical curve (P = 0.1804; ).

Fig. 3.

Tetraodontiform body size and temperature over time. Sea surface temperature for tropical latitudes (15°N to 15°S; orange line) and a global average sea temperature (yellow line) are plotted for the past 100 Ma. The reconstructed ancestral node body size (log mean maximum SL in centimeters) for tetraodontiforms is also plotted against time. Sea surface temperatures have been slowly cooling since the Late Cretaceous, while tetraodontiform body size has gradually increased. provides a version of this plot colored by family.

Discussion

By integrating morphological and genome-wide sequence data in a total-evidence dating framework, we inferred the most complete phylogeny for both living and extinct tetraodontiform fishes. Using this robust phylogenetic framework, which accounts for uncertainty in topology and divergence time estimates, we conducted a suite of comparative approaches to test for Cope’s and Bergmann’s rules. Our model-fitting analyses that incorporate paleoclimatic information provided strong support for the Cope–Bergmann rule (8), where tetraodontiform body size is strongly correlated with sea surface temperatures over the past 100 Ma. Furthermore, by decomposing the temperature curve into an overall trend (Cope’s rule), the fluctuations around that trend (Bergmann’s rule), and a combination of both (Cope–Bergmann rule) (8), we found strong support for the latter. Our results align with previous studies that examined paleoclimatic effects on species evolutionary trends, providing support for a correlation between temperature and body size. Studies of birds and mammals have found support for both rules (13, 39, 40). However, among ectotherms, findings have been mixed, with studies identifying support for Cope’s rule in arthropods (19, 41) and reptiles (42); support for Bergmann’s rule in fishes (10, 43), arthropods (8, 11, 19), reptiles (44), and amphibians (5, 12); an inverse or no support for Cope’s rule in arthropods (45) and reptiles (23, 46); and an inverse or no support for Bergmann’s rule in arthropods (47), reptiles (44), and amphibians (48, 49). Bergmann’s rule, originally hypothesized for endothermic vertebrates (50–52), describes how larger endothermic species might conserve heat better in higher latitudes with cooler temperatures due to an increased surface area to volume ratio. The heat conservation hypothesis behind this rule would likely not apply to ectotherms at all [(53, 54); but see (55)]; thus, other explanations or hypotheses for this pattern are necessary (see below).

Studies examining paleoclimatic effects on evolutionary trends benefit from a robust fossil record. However, variables such as habitat composition, sampling effort, and specimen morphology can influence fossil preservation. Hard-bodied organisms inhabiting shallow marine environments are more likely to fossilize compared to soft-bodied, deep-sea species (56), and there are additional taphonomic biases related to body size. Among fishes, larger fossils are more likely to become disarticulated and then scattered by various hydrological processes, resulting in a lower probability of discovery (57). In contrast, smaller specimens have fragile bones and thus a higher potential of being destroyed compared to larger specimens; thus, they may not be as common in the fossil record (58, 59). All analyses that incorporate fossil data will have taphonomic biases, but a goal should be to minimize potential biases whenever possible. Tetraodontiform fishes are well represented in the fossil record, likely owing to their hardened external anatomy and habitat preferences for shallow marine waters. Their fossil record is rich, with extinct representatives in all 10 living families, as well as 12 exclusively fossil families (4). Among the tetraodontiform fossils in our dataset, most (51 specimens, representing 14 species in 10 families) come from the Monte Bolca Lagerstätten (50.5 to 48.5 Ma). Within this Eocene locality, tetraodontiform fossils are exceptionally preserved (60) and body sizes range from 8 mm SL (Eolactoria sorbinii) to 521 mm SL (Protobalistum imperialis) (Dataset S3). Given the large range in sizes that are found in Monte Bolca and their overall completeness, it is unlikely that preservation and size-related taphonomic biases have major effects on our analyses.

Cope’s rule is often explained as a selective advantage toward larger body sizes. Despite the fact that overall body size increased over time, we found that support for a mean trend model (i.e., Cope’s rule), while higher than that of the OU, rate trend, EB, and BM models, was not the strongest, implying that a selective advantage to large body sizes over time alone does not provide enough explanatory power. Indeed, we found support for a Cope–Bergmann rule (8) in our climate OU models where body size increased over time as sea surface temperatures gradually cooled. Cope’s rule was previously thought to be a manifestation of Bergmann’s rule during periods of climate cooling (Cope–Bergmann rule) (8), and any trend in increased body size during these cooler periods would just be seen as coincidental with Cope’s rule; thus, the two rules might be combined into a single Cope–Bergmann rule. Contrary to one of the three predictions of the Cope–Bergmann rule postulated by Hunt and Roy (8), which suggests that increased body sizes are attributable to evolutionary changes within lineages/species, the evolutionary trajectories observed in tetraodontiforms may be indicative of lineage turnover resulting from species selection. While, unlike Hunt and Roy (8), we lack a high-resolution temporal dataset to assess intraspecific evolution, a clear general pattern emerged from our dataset: as smaller and older tetraodontiforms species became extinct, they were replaced by larger and younger ones. Ultimately, these patterns may be indicative of Bergmann’s rule operating at deeper macroevolutionary scales (9, 61).

In other marine ectotherms, temperature is a driving factor for increased body size. Hunt et al. (19) examined trends in 19 species of deep-sea ostracods and found that body size increased in 84% of examined species from the Eocene (40 Ma) to the present, during which global deep-sea temperatures gradually cooled. They also identified strongest support for a temperature-tracking model (a simple model where expected changes in body size correspond to changes in temperature), finding significant body size increases only during periods of sustained cooling. Taken together, these results highlight the importance of incorporating environmental variables into macroevolutionary studies that examine trends over time.

It is unlikely that a global explanation for body size evolution exists for all ectotherms in response to paleoclimate change and, by extension, Bergmann’s rule. However, temperature seems to be an important driver of their evolutionary patterns, as ectothermic species have a reduced capacity for heat conservation compared to endotherms. An increase in body size in colder environments may be related to the temperature-size rule, which states that ectotherms reared in colder temperatures in experimental settings tend to grow more slowly and to mature with larger body sizes relative to ectotherms reared under warm temperatures (62). Paleotemperature changes over large timescales are correlated with changes in ectotherm body size (8, 38). Additionally, because many environmental variables are temperature dependent, temperature may play an indirect role in determining ectotherm body size, with other confounding factors coming into play. Among aquatic ectotherms, oxygen may be a more limiting factor for body size evolution. Reduced dissolved oxygen in warmer waters is limiting for ectotherms dependent on aquatic respiration. This temperature-dependent oxygen limitation has been proposed to relate to Bergmann’s rule in aquatic ectotherms (14), and these reductions in body size are greater in aquatic taxa than in terrestrial taxa (63). Other studies have linked warming-induced anoxia as a driver of decreased body size in Early Jurassic marine invertebrates (20). Temperature-dependent oxygen limitation may also explain why past studies of terrestrial ectotherms found either an inverse Bergmann’s rule (48, 49) or a nonexistent trend (44, 47), meaning that increased body size in relation to temperature cooling may be much stronger in aquatic-respiring ectotherms (14).

Past climatic changes strongly shaped the trajectory of phenotypic evolution across many clades. On a global scale, large paleoclimate changes are associated with extinction events, such as those that occurred during the end-Permian (ca. 252 Ma) and Cretaceous–Paleogene (ca. 66 Ma) events (64). These events can be size selective and are often thought to favor small-bodied taxa (65, 66). In general, small-bodied species are presumed to be at a lower risk of extinction, potentially owing to their faster generation times and increased fecundity (24). But they may also be at increased risk of extinction due to factors such as geographic range, which is often smaller than that of their large-bodied counterparts (66). Although extinction events favoring small-bodied taxa have been documented (24, 66, 67), the opposite pattern has been observed as well, implying that this may simply be a clade-specific effect (66). When extinction risk is examined at higher taxonomic levels, these events tend to favor large-bodied taxa (65), suggesting that additional biological factors play a role in extinction risk [e.g., lesser predation of large-bodied tetraodontiforms due to increased toxicity of flesh (68)] and scale is an important factor to consider. Among ectotherms, those in marine environments may be most affected by global extinction events related to climate warming compared to terrestrial habitats (69). This was the case for the largest mass extinction event, the end-Permian, where temperature-induced hypoxia drove a majority of marine species to extinction (64). Examining these past patterns of ectotherm evolution may provide insights into how species will respond during the next chapter of global climate change.

In conclusion, we found a strong link between the evolution of body size in tetraodontiforms and past climate and paleotemperatures. Gradual climate cooling over the past 100 Ma (especially pronounced during the past 50 Ma) was associated with increases in average tetraodontiform body length. Our results are robust to a number of factors driving uncertainty in macroevolutionary inferences, including the use of different genomic subsets and root priors for time tree inference, the inclusion or exclusion of the controversial plectocretacicoid fossils, the utilization of different paleotemperature curves, and the implementation of alternative comparative approaches. While the evolution of body size in tetraodontiforms appears to conform to the Cope–Bergmann rule, other factors (e.g., ocean acidification, dissolved oxygen concentrations) could affect this trait and thus deserve further investigation.

Materials and Methods

Taxonomic Sampling and Genomic Data.

Extended materials and methods are reported in the . We generated genomic data for 141 individuals representing 131 species of the order Tetraodontiformes and four species of its sister group, the Lophiiformes (33) (Dataset S1). All tissue samples were associated to voucher specimens deposited in museum collections (). We shipped DNA extractions to Arbor Biosciences for library preparation and target enrichment. Sequencing of libraries was conducted at The University of Chicago Genomics Facility (Illumina HiSeq 4000). Target capture used the Eupercaria probe set of Hughes et al. (26, 70) to enrich 1,105 single-copy nuclear exon markers. We assessed sequence quality and removed two exons due to high levels of missing data, leaving 1,103 exons in total. We excluded one newly sequenced species (Rhinecanthus verrucosus) due to low capture efficiency. After quality control, we aligned all exons by considering their reading frames. We further increased our taxonomic sampling by adding sequences from 55 additional tetraodontiform species and one outgroup species retrieved from the National Center for Biotechnology Information (Dataset S2). Our combined genomic dataset contained 185 tetraodontiforms and five outgroup taxa.

Phylogenomic Inference.

We inferred phylogenetic trees and associated support values in an ML framework in IQTREE v.1.6.12 (71) (Dataset S4). In addition, we conducted a multispecies coalescent analysis in ASTRAL-III (72) based on IQTREE ML gene trees ( and Dataset S5). To account for the effect of missing data in our dataset, we conducted two concatenation-based ML analyses using all 1,103 exon markers: one including all newly sequenced taxa and four previously published transcriptomes (134 tetraodontiform species, 47% missing data overall; and Dataset S6) and a second excluding taxa with more than 65% missing data (102 tetraodontiform species, <33% missing data overall; and Dataset S7). Because the topology and branch lengths were largely in agreement between the two analyses, all downstream phylogenetic analyses used the complete dataset.

Integration of Fossil and Extant Species.

To combine the fossils and extant tetraodontiform species, we used the morphological matrix of Arcila and Tyler (4), which consists of 210 characters coded for 17 extant and 52 fossil tetraodontiform species, plus two additional outgroup taxa. We combined the morphological matrix with our genomic dataset for a total of 237 tetraodontiform species and seven outgroups. Our analyses used the GTRGAMMA and Mk models with four partitions, three for the molecular sequences (one for each codon position) and one for the morphological dataset.

Phylogenetic Uncertainty and Total-Evidence Dating Using the Fossilized Birth Death (FBD) Process.

In addition to the phylogenomic analyses described above, we conducted divergence time estimations under a total-evidence, or tip-dating, framework using the FBD model in MrBayes v 3.2.7a (73). To account for topological uncertainty, we assembled 15 largely independent genomic subsets containing ca. 50 randomly selected loci subsampled from the complete genomic dataset (1,103 loci total) (Dataset S8). All subsets overlapped in only five “anchor” genes to maintain the same set of species for each subset (29, 30). In addition to genomic data, each subset contained the morphological dataset with fossil and extant taxa. We provide a list with the fossil ages in Dataset S3 as well as a list of prior distributions used for node dating from previous studies that included Tetraodontiformes in their analyses (). We ran all 15 subsets in MrBayes. After 6 mo of total runtime, only five (of the 15) subsets reached convergence based on estimated sample size values close to or above 200, and we retained those for downstream comparative analyses. Because there is no consensus yet on whether the superfamily Plectocretacicoidea should be considered as stem tetraodontiforms and the exclusion of this superfamily has the potential to drastically affect the age estimations (4, 31), we used two different calibration schemes, including and excluding plectocretacicoids. We sampled 100 trees from the posterior distribution of each retained subset (500 trees total). In addition, we constructed an MCC tree from 10,000 trees evenly sampled from the posterior of all five subsets using TreeAnnotator v.2.6.0 (74).

Trait Data.

We compiled standard length (SL) data for most fossil and extant tetraodontiform species in our dataset using museum collection databases, public datasets, and published articles (Dataset S3). Given the highly fragmented nature of some fossil specimens, we excluded three extinct species from the body size analysis. Additionally, because of the bias for smaller specimens in museum collections, we omitted any measurements from individuals that were more than 20% smaller than the maximum recorded size and averaged the largest specimens to obtain a mean maximum SL per species. We performed all analyses using log-transformed values. Mean maximum length was chosen as an indicator for how large a species could potentially reach. Although some tetraodontiform clades have a more three-dimensional body shape compared to other fish groups, PGLS analyses between SL and volume (SL versus volume: P = 0.0005) and between SL and surface area (SA; SL versus SA: P = 0.0248) revealed a strong correlation (). These analyses were restricted to a subset of 41 tetraodontiform species (across all 10 extant families) for which computed tomography scans (CT-scans) were generated or obtained from Morphosource (Dataset S9). Measurements other than SL were not included due to homology concerns arising from fragmented fossil specimens.

Paleotemperature Data.

We obtained temperature curves that spanned the nearly 100 Ma evolutionary history of tetraodontiforms from Scotese et al. (2). These authors used oxygen isotope data to reconstruct past global average ocean temperatures and sea surface temperatures between tropical latitudes (15°N to 15°S). These two contrasting temperature curves most accurately captured the spectrum of past environmental affinities in this group (see ).

Evolutionary Model Fitting.

We conducted model fitting analyses in R version 4.0.2 (75). We fitted models of continuous character evolution using the R package “geiger” and the “fitContinuous” function. To more explicitly assess an increase in body size in response to paleoclimatic changes over time (Cope–Bergmann rule), we fitted a climate-dependent OU model, where the parameter of the model that represents the optimum trait value is time-variable and follows a temperature curve (76). We tested this model using the two temperature curves mentioned above from Scotese et al. (2). Additionally, we tested the climate OU model on a deep-sea curve from Cramer et al. (35). The and provide a description of other models tested. To account for tree uncertainty, we tested all models with the 500 trees selected from the posterior distribution. We accounted for interspecific sources of variation by incorporating measurement error into our model fitting analyses ( provides more details).

Ecomorphological Correlations.

To further determine whether tetraodontiforms adhere to the Cope–Bergmann rule, we examined patterns of body size in relation to past ocean temperature changes for the two temperature curves mentioned above. We performed ancestral state reconstructions of body sizes for all nodes and mapped these onto the MCC tree using the “contMap” function in the R package “phytools” (36). Resulting estimated ancestral sizes for each node were incorporated into a PGLS analysis and modeled for each of the two temperature curves from Scotese et al. (2). We tested various evolutionary models for the PGLS analysis (e.g., OU model, BM model, and a nonphylogenetically informed model [i.e., ordinary least squares]) to determine best fit. ( provides more details).