Periodic Environmental Disturbance Drives Repeated Ecomorphological Diversification in an Adaptive Radiation of Antarctic Fishes

The ecological theory of adaptive radiation has profoundly shaped our conceptualization of the rules that govern diversification. However, while many radiations follow classic early-burst patterns of diversification as they fill ecological space, the longer-term fates of these radiations depend on many factors, such as climatic stability. In systems with periodic disturbances, species-rich clades can contain nested adaptive radiations of subclades with their own distinct diversification histories, and how adaptive radiation theory applies in these cases is less clear. Here, we investigated patterns of ecological and phenotypic diversification within two iterative adaptive radiations of cryonotothenioid fishes in Antarctica’s Southern Ocean: crocodile icefishes and notoperches. For both clades, we observe evidence of repeated diversification into disparate regions of trait space between closely related taxa and into overlapping regions of trait space between distantly related taxa. We additionally find little evidence that patterns of ecological divergence are correlated with evolution of morphological disparity, suggesting that these axes of divergence may not be tightly linked. Finally, we reveal evidence of repeated convergence in sympatry that suggests niche complementarity. These findings reflect the dynamic history of Antarctic marine habitats and may guide hypotheses of diversification dynamics in environments characterized by periodic disturbance.


Introduction
For more than a century, adaptive radiations have provided foundational empirical systems for investigating the factors that promote the diversification of life (e.g., Osborn 1902;Simpson 1944Simpson , 1953Grant 1986;Schluter 2000;Losos 2009). Evolutionary biologists have been particularly interested in understanding the circumstances underlying temporal dynamics of speciation and phenotypic diversification over the course of a radiation. A common expectation holds that adaptive radiations undergo an initial increase in the rates of speciation and phenotypic diversification as new adaptive zones are explored, followed by a decline in diversification rates with increased saturation of adaptive zones and increased specialization of lineages (Simpson 1944(Simpson , 1953, and this expectation has been supported by numerous empirical studies (e.g., Slater et al. 2010;Arbour and Fernández 2013;Price et al. 2016;García-Navas et al. 2018;Borko et al. 2021). However, some radiating clades have endured complex geophysical and climatic contexts that may have promoted the generation of iterative adaptive radiations of subclades with their own unique diversification histories. For example, many regions of the world have experienced climatic oscillations that significantly affected the biotic and abiotic environments. As environmental changes repeatedly drove some lineages extinct, the freeing of ecological space subsequently provided opportunities for new phases of speciation and diversification (Near et al. 2012;Ivory et al. 2016;Pouchon et al. 2018). While it is natural to expect that such cases could produce evolutionary patterns distinct from a single burstlike adaptive radiation, it is less clear what we should expect in such scenarios. A more nuanced understanding of these iterative adaptive radiations is therefore critical for understanding how extinction and changing environmental conditions shape major components of adaptive radiations following their early evolutionary origins.
The iconic cryonotothenioid fishes of Antarctica's Southern Ocean represent an exemplar clade in which to investigate the factors modulating the shape and pace of an adaptive radiation. Cryonotothenioids comprise nearly 80 phenotypically and ecologically diverse species that dominate the abundance, biomass, and species diversity of the region's nearshore teleost fauna (Eastman 1993(Eastman , 2005. The remarkable diversification of these fishes has been attributed to periodic cycles of warming and cooling, which have profoundly shaped the physical environment and biodiversity of the Southern Ocean (O'Loughlin et al. 2011;Near et al. 2012;Wilson et al. 2013;Dornburg et al. 2016Dornburg et al. , 2017. The cyclical advances and retreats of grounded ice across submerged portions of the Antarctic continental shelf (i.e., ice scour) have periodically obliterated benthic nearshore habitat across thousands of kilometers of the Antarctic continental shelf and likely annihilated entire shelf communities (Thatje et al. 2005(Thatje et al. , 2008Smale et al. 2008;Tripati et al. 2009;Allcock and Strugnell 2012;Near et al. 2012;Dornburg et al. 2017). Although laying waste to the incumbent fauna, these ice scour events generate ecological opportunities via the emptying of previously occupied niches. This "resetting of the ecological stage" is thought to have produced a pattern of iterative adaptive radiations nested within the larger cryonotothenioid radiation (Near et al. 2012). The diversification of these subclades likely accounts for the exceptionally rapid speciation rates exhibited by cryonotothenioids relative to other marine teleosts (Rabosky et al. 2018). However, it remains unclear whether and to what extent patterns of phenotypic and ecological divergence within these nested radiations have been impacted by cyclical glacial dynamics.
Most investigations of cryonotothenioid diversification include limited taxon sampling and focus on characterizing patterns across the early radiation. This taxonomically restricted view hinders a detailed understanding of factors that have contributed to iterative adaptive radiations within cryonotothenioids after their early diversification. Recent studies that increased taxonomic sampling revealed that diversification of cryonotothenioids into water column niches spanning the benthic-pelagic axis occurred multiple times in parallel within different cryonotothenioid subclades, suggesting that periodic ice scour disturbances have generated multiple bouts of ecological opportunity for diversification (Rutschmann et al. 2011;Near et al. 2012;Dornburg et al. 2017). However, to our knowledge, no studies have explicitly examined the relationships among phenotypic and ecological parameters within cryonotothenioid subclades, precluding an understanding of whether aspects of phenotypic and ecological evolution are coupled or whether different lineages have evolved distinct phenotypic strategies to occupy similar niches. Moreover, the extent to which the presence of closely related species that co-occur in a given community influences patterns of trait evolution remains unclear. Investigating the influence of both evolutionary history and potential competitive interactions is critical if we want to understand the diversification history of this adaptive radiation and gain new insights into how biological communities recover in the face of periodic disturbance.
Here, we investigate the patterns of phenotypic and ecological diversification within two cryonotothenioid subclades that have been identified as nested parallel radiations (Near et al. 2012): the crocodile icefishes (Channichthyidae) and the notoperches (Trematominae sensu Near et al. [2018], exclusive of the nested non-Antarctic clade Patagonotothen). We compile data on buoyancy, depth occupancy, feeding ecology, body shape and size, and phylogeny to evaluate hypotheses regarding the evolutionary dynamics of ecological and phenotypic trait disparity. We first use phylogenetic least squares regression to evaluate the extent to which the evolution of key phenotypic and ecological traits is correlated. Next, we visualize phylogenetic patterns of trait disparity and use disparity-through-time analyses to evaluate support for alternative expectations of phenotypic and ecological trait diversification within each of our focal radiations. Finally, we test the prediction that interactions among co-occurring species drive elevated trait divergences in adaptive radiation (Dayan and Simberloff 2005;Pfennig and Pfennig 2009). Collectively, our results provide a novel perspective on the patterns of ecomorphological diversification within the cryonotothenioid adaptive radiation.

Trait Data Collection
We compiled data on a series of traits that are hypothesized as central to the adaptive radiation of cryonotothenioids: mean percentage buoyancy, which functions as a proxy for water column niche (Near et al. 2012;Eastman 2017); depths of occurrence (table S1); frequency of occurrence of prey items (figs. S1, S2; table S2); geometric morphometric data on body shape (fig. S3; table S3); and  maximum body size (table S4), represented using the maximum total length (TL) recorded for each species Near et al. 2012;Colombo et al. 2015;Dornburg et al. 2017;Eastman 2017Eastman , 2019Eastman , 2020Daane et al. 2019). Before conducting downstream analyses of trait diversification patterns, we identified major axes of variation in diet using a phylogenetic principal components analysis (PCA) of the prey frequency occurrence data using the R package phytools (Revell 2011). Additionally, the geometric morphometric data on body shape were first transformed using a generalized Procrustes analysis (GPA; Rohlf and Slice 1990;Zelditch et al. 2012) and then subjected to a phylogenetic PCA in the R package geomorph (Adams et al. 2021;Baken et al. 2021). Full details of the trait data collection are provided in "Supplementary Methods" in the supplemental PDF. All data have been deposited in the Dryad Digital Repository (https://doi.org/10.5061 /dryad.bk3j9kdbg; Parker et al. 2022).

Divergence Time Estimation
To estimate a time-calibrated phylogeny of cryonotothenioids for comparative analyses, DNA alignments were taken from Near et al. (2018), who used restriction siteassociated DNA sequencing to capture 104,709 loci across 80 species of notothenioids. Divergence times were estimated using BEAST (ver. 2.4.7;Drummond et al. 2012;Drummond and Bouckaert 2015) assuming an uncorrelated model of molecular rates with a lognormal distribution and a birth-death prior on branching times. We placed four prior age constraints used in previous studies of notothenioid divergence times (see "Supplementary Methods" in the supplemental PDF; Near et al. 2012;Dornburg et al. 2017). Adequate sampling of the posterior distribution for each parameter of each run was assessed via computation of effective sample size (ESS) values, with ESS values 1200 indicating adequate sampling. Marginal posterior probabilities versus generation state were assessed for each parameter in the program Tracer (ver. 1.6; http://beast.bio.ed.ac .uk/Tracer). A maximum clade credibility (MCC) tree was generated using TreeAnnotator (ver. 2.4.7; http://beast.bio .ed.ac.uk/TreeAnnotator) to summarize the posterior probability density of topological and branch-length estimates. For downstream phylogenetic comparative analyses, we extracted from the MCC tree one subtree including all 14 species representing the crocodile icefishes and one subtree including all 16 species representing the notoperches ( fig. S4). Full details of our phylogenetic analyses can be found in "Supplementary Methods" in the supplemental PDF.

Testing for Evidence of Correlated Evolution
To investigate potential evolutionary correlations among body shape and major axes of ecological disparity in ice-fishes and notoperches, we performed multivariate phylogenetic generalized least squares (PGLS) regression using the procD.pgls function in the R package geomorph (Adams et al. 2021;Baken et al. 2021). This approach enabled us to test the relationships between morphology and ecology while accounting for the expected covariance of traits among species due to shared ancestry. In each of our tests, body shape variation was the response variable and was represented using the full set of species-averaged Procrustes shape coordinates from the GPA. We first tested whether changes in body size (represented by maximum TL) are associated with predictable changes in body shape. We then assessed relationships between changes in body shape and changes in habitat and prey usage by conducting a series of regressions of body shape on mean percentage buoyancy, mean depth occupancy, and each of the first three PC axes of diet variation from PCA of the prey frequency occurrence data.
To test whether changes in ecology predict changes in only the major axes of body shape variation, we fitted a series of linear regression models in which each of the first three PC axes of body shape variation from a PCA of the Procrustes shape coordinates was used as the response variable and each of the ecological variables (including mean percentage buoyancy, mean depth, and each of the first three PC axes of diet variation) was used as the predictor variable. We additionally tested whether changes in body size predict changes along each of the first three PC axes of body shape variation as well as whether changes in our ecological variables predict changes in body size. Furthermore, we assessed relationships among each of the ecological variables used to represent habitat and prey usage for icefishes and notoperches. We first tested whether mean percentage buoyancy and mean depth are correlated, given that both traits represent major axes of variation in habitat utilization in cryonotothenioids (Eastman 2017(Eastman , 2020. We then evaluated whether changes in habitat utilization (represented using either mean percentage buoyancy or mean depth) are associated with predictable changes in prey usage (represented by the first three PC axes of diet variation). For all of these tests, we used the gls function in the nlme R package (Pinheiro et al. 2017) to fit linear models of trait correlation in which the covariance structure was defined according to either a Brownian motion model or a Ornstein-Uhlenbeck model of trait evolution on each clade's phylogeny (specified using either the cor-Brownian function or the corMartins function, respectively, in the ape R package; Paradis et al. 2004). Across all tests of correlated evolution, tests were performed independently for icefishes and for notoperches, and r and P values were adjusted using the false discovery rate correction method of Benjamini and Hochberg (1995) implemented using the p.adjust R function.

Visualizing Evolutionary Dynamics of Trait Diversification
We used a series of visualizations to qualitatively evaluate the diversification history of each trait in this study. First, patterns of divergence in bathymetric distribution were visualized using ridgeline plots (Fitzpatrick et al. 2019) in proximity to each clade's time tree. This approach allowed us to directly contrast the range of the water depths occupied by each lineage. We then assessed divergences along major axes of variation in diet by generating phenograms, which project a phylogeny into trait space defined by each of the first three PC axes of variation in prey usage. Phenograms were also used to visually assess variation in body size as well as patterns of morphospace occupancy, where each of the first three axes of body shape variation from the PCA were plotted along with the phylogeny of each clade. All phenograms were generated using the R package phytools (Revell 2011).
In addition to visualizing phylogenetic patterns of trait disparity, we used the R package geiger (ver. 2.0; Pennell et al. 2014) to calculate relative subclade disparity through time (DTT; Harmon et al. 2003) for mean percentage buoyancy, mean depth, and maximum body size. This approach allowed us to test the competing expectations that historical patterns of adaptive radiation have left a signature of among-subclade disparity in the early history of these clades versus an expectation of high within-clade disparity driven by convergence. We restricted our analyses to these traits, as ordination analyses have the potential to bias major axes of variation toward patterns consistent with early-burst models (Uyeda et al. 2015). We compared the empirical DTT results with the expectations of a null model of Brownian motion generated from 10,000 simulations (Slater et al. 2010).

Testing the Relationships between Range Overlap, Clade Age, and Ecomorphological Diversification
Interactions between species living in sympatry are thought to promote increased divergence in ecological traits (Pfennig and Pfennig 2009;Tobias et al. 2014). Alternatively, elevated trait divergences may reflect greater evolutionary time since two species shared a common ancestor. To test whether levels of ecomorphological divergence among icefishes and notoperches can be predicted by range overlap (a proxy for interspecific interactions) or time since common ancestry, we fitted a series of phylogenetic generalized linear mixed models (PGLMMs) using the R package MCMCglmm (Hadfield 2010). In each model, our response variable corresponded to divergence in a given ecological or phenotypic trait for all possible pairwise comparisons of species within each of the icefish and notoperch clades.
We used the dist function in R to calculate pairwise divergence as the Euclidean distance between species in mean percentage buoyancy, mean depth, maximum TL, and PC scores for each of the first three PC axes of variation in diet and in body shape. We tested the effects of three fixed predictors on ecomorphological divergence: (1) the range overlap for a given species pair (modeled as a two-level factor: sympatric, allopatric), (2) the time since common ancestry (millions of years) for a given species pair, and (3) the interaction between range overlap and time since common ancestry. Our models additionally accounted for two random effects: (1) the nonindependence of trait values arising from shared ancestry of lineages and (2) the nonindependence of data points arising from repeated measurements of each lineage (Tobias et al. 2014). We accounted for phylogenetic nonindependence by fitting as a random effect the variancecovariance (VCV) matrix calculated from each clade's MCC tree and connecting the VCV matrix to data points via the most recent common ancestral node shared by a given species pair. Nonindependence arising from repeated measurements was accounted for by fitting species identifiers (as either the focal lineage or the comparison lineage in each pairwise comparison) as random effects.
Model fitting was conducted independently for notoperches and for icefishes, and we fitted separate linear mixed models for each trait. Inverse gamma priors were applied to the variance structure of the random effects. For all models, the posterior distributions of model parameters were simulated using 20 million Markov chain Monte Carlo generations, with sampling every 10,000 generations after discarding the first 10% of generations as burn-in. Model convergence was assessed by visualizing trace plots for each of the estimated parameters in each model and by ensuring that ESS values exceeded 200 for each parameter estimate. Full details of these analyses can be found in "Supplementary Methods" in the supplemental PDF.

Correlated Evolution
Results of our PGLS regressions revealed no evidence that body size or any of the major axes of diet variation strongly influence variation in the Procrustes body shape coordinates across icefishes (P 1 :05 for all regressions; table S5) or notoperches (P 1 :05 for all regressions; table S6). Similarly, for both icefishes and notoperches, our PGLS regressions reveal no evidence that divergences along any of the major PC axes of body shape variation are correlated with divergences in diet, depth occupancy, or percentage buoyancy (P 1 :05 for all regressions; tables S7-S10). We additionally find no evidence that changes in diet, depth occupancy, or buoyancy are associated with predictable changes in maximum body size (P 1 :05; tables S7-S10). Finally, our results do not uncover significant relationships among any of the major ecological axes of divergence in the icefish and notoperch radiations. Specifically, variation in prey usage is not predicted by either depth occupancy or habitat use along the benthic-pelagic axis (P 1 :05), and divergence along the benthic-pelagic axis does not predict divergence in depth occupancy (P 1 :05; tables S7-S10).

Diversification of Feeding Ecology
The first three PC axes explained approximately 95% of the variation in icefish diet ( fig. S5). PC1 (66%) captured variation largely between the consumption of cephalopods and krill ( fig. S5). PC2 (24%) and PC3 (5%) similarly reflected changes in the consumption of krill, cephalopods, and fishes, with PC3 further including differences in the consumption of munids ( fig. S5). Projecting the phylogeny into trait spaces defined by each of these three PC axes reveals striking patterns of convergence in feeding ecology between distantly related lineages ( fig. 1A-1C). For example, Dacodraco hunteri, Chionodraco hamatus, Chionodraco myersi, and Channichthys rhinoceratus all converge on a nearly identical point in the trait space defined by the first PC axis of diet variation, reflecting a convergence on piscivory ( fig. 1A). Similarly, Pagetopsis maculatus and Chaenodraco wilsoni converge on a point in the space defined by the second PC axis of diet, reflecting increased krill consumption ( fig. 1B). In contrast, closely related species and sister species pairs often appear highly divergent along any of the first three PC axes. Examples include the divergences between Champsocephalus esox and Champsocephalus gunnari, C. wilsoni and C. myersi, and Chionobathyscus dewitti and Cryodraco antarcticus ( fig. 1A-1C).
In notoperches, the first three PC axes explain 70% of the variation. PC1 (38%) captured variation largely in consumption of Hyperiella, copepods, and polychaetes ( fig. S6). PC2 (19%) and PC3 (13%) similarly show variation in marine invertebrate prey, including amphipods ( fig. S6). Similarly to icefishes, projections of the phylogeny into trait spaces defined by each of the first three PC axes of diet variation reveal patterns of convergence in feeding ecology between distantly related species ( fig. 1D-1F). For example, in the diet space generated by PC1, distant lineages Trematomus loennbergii and Trematomus scotti converge on polychaete consumption ( fig. 1D), and in the diet space generated by PC2, Nototheniops cf. nudifrons, Trematomus eulepidotus, and Trematomus newnesi all converge on a space between copepods and euphausiids ( fig. 1E). Furthermore, the first three PC axes show recently diverged notoperch lineages diverging in diet space, including between the sister pair T. loennbergii and T. lepidorhinus and between their common ancestor and T. eulepidotus ( fig. 1D-1F).

Major Axes of Body Shape Evolution
The first three PC axes explained approximately 68% of body shape variation in icefish ( fig. S7A). PC1 (accounting for 37% of variance) is primarily associated with variation in body depth, PC2 (19% of variance) corresponds to variation in mouth position, and PC3 (12% of variance) captures variation in snout length (figs. 4A-4C, S7A). Our phenograms reveal several instances of convergence of distantly related icefish species in aspects of body shape. For instance, C. antarcticus, C. esox, and C. gunnari converge on a relatively slender, elongate body shape, while Pseudochaenichthys georgianus and Pagetopsis macropterus converge on deeper-bodied forms ( fig. 4A-4C). However, examination of the first three PC axes also reveals striking divergences in aspects of body shape among closely related species, including sister species pairs P. macropterus and P. maculatus, C. antarcticus and C. dewitti, and C. hamatus and C. rastrospinosus ( fig. 4A-4C).
In notoperches, the first three PC axes explain 62% of the variance in body shape ( fig. S7B). PC1 (39% of variance) is primarily associated with variation in mouth position, PC2 (16%) corresponds to variation in head depth and body depth at anal fin origin, and PC3 (7%) describes variation in caudal peduncle length (figs. 4D-4F, S7B). Similar to our observations for icefishes, we observe multiple instances of trait convergence among distantly related notoperch species. For instance, Nototheniops larseni and T. lepidorhinus converge on a relatively upturned mouth position, and T. newnesi, T. nicolai, and T. bernacchii all converge on a relatively slender, elongate body shape ( fig. 4D-4F). Meanwhile, we find evidence of high divergence in body shape among closely related species, as exemplified by the pairs T. lepidorhinus and T. loennbergii, T. pennellii and T. nicolai, and T. bernacchii and T. hansoni ( fig. 4D-4F).

Body Size Variation
Similar to patterns of disparity in diet, habitat utilization, and body shape, we observe complex patterns of evolutionary  Figure 1: Evolutionary histories of the top three principal component (PC) axes of diet variation reveal complex patterns of repeated convergence in prey usage among distantly related species and high divergence among closely related species. A-C represent projections of the crocodile icefish phylogeny, and D-F represent projections of the notoperch phylogeny. In all panels, the x-axis represents time in millions of years, with branch lengths reflecting estimated divergence times. The y-axis reflects variation along each of the first three PC axes of diet, with tip placement along the y-axis reflecting the PC score for that species along each PC axis.
convergence and divergence in body size among icefishes and notoperches. Within icefishes, distantly related species D. hunteri and the P. macropterus-P. maculatus species pair converge on relatively small maximum body sizes (25-33 cm), while C. gunnari (68 cm) and Chaenocephalus aceratus (76 cm) converge on maximum body sizes approximately two times greater than those recorded for the smallest icefish species ( fig. 5B). High divergence in body size is also obvious between closely related species, such as between C. gunnari (68 cm) and C. esox (35 cm) and between P. georgianus (60 cm) and P. maculatus (25 cm; fig. 5B). In notoperches, we find several closely related species to exhibit high divergence in body sizes, including T. hansoni (45.5 cm) and T. tokarevi

Patterns of Disparity through Time
Analyses of the average relative subclade disparity through time (DTT) provide little general support for an early partitioning of the overall disparity among subclades ( fig. 6). High variance across Brownian motion simulations probably results from small clade sizes and limits the power to statistically reject one hypothesis over the other ( fig. 6). However, the trends revealed by DTT plots mirror the ones depicted in trait phenograms (figs. 1, 4, 5). Analysis of buoyancy disparity through time reveals levels of relative subclade disparity that are higher than the mean trend of DTT expected for icefishes (morphological disparity index [MDI] p 0:30, P p :90; fig. 6A), running counter to the trend for depth (MDI p 20:18, P p :083; fig. 6B). The latter represents the only case of early low relative subclade disparity and contrasts sharply with body size DTT, which depicts a signature of within-clade disparity consistent with high levels of between-clade trait convergence (MDI p 0:36, P p :92; fig. 6C). Notoperches similarly reveal trends more consistent with convergent evolution in buoyancy (MDI p 20:037, P p :32; fig. 6D) and depth (MDI p 0:25, P p :79; fig. 6E) during their evolutionary history. Likewise, the evolution of body size disparity reflects a signal of high levels of convergence in body size between clades rather than a partitioning of body size diversity between clades (MDI p 20:018, P p :41; fig. 6F).

The Relationships between Range Overlap, Clade Age, and Ecomorphological Diversification
For all pairwise comparisons of lineages within icefishes and notoperches, we generally find that range overlap, clade age, and the interaction between range overlap and clade age do not significantly predict observed levels of ecological or phenotypic divergence (table S11, pts. a-g; table S12, pts. a-i). However, within notoperches, patterns of divergence in several traits exhibit exceptions to this general finding. First, our PGLMMs reveal that pairwise divergences in the third PC axis of body shape variation (table S12, pt. e) and in body size (table S12, pt. i) are higher across sympatric species pairs relative to allopatric species pairs. Second, we find evidence for significant positive relationships between clade age and divergence in buoyancy (table S12, pt. a), depth occupancy (table S12, pt. b), and body size (table S12, pt. i). Finally, we find a significant negative relationship between divergence in depth and the interaction between range overlap and clade age (table S12, pt. b).

Discussion
Here, we integrated a strongly supported species-level phylogenetic hypothesis with a comprehensive data set on body shape and size, diet, depth, and buoyancy to elucidate patterns of ecomorphological diversification within two cryonotothenioid subclades that have been identified as iterative adaptive radiations: the crocodile icefishes and the notoperches. Our results reveal several unique patterns of adaptive radiation that have likely been influenced by the high-frequency disturbance regime characteristic of Antarctica's continental shelf. First, we find limited evidence that changes in body shape and size are correlated with divergences along axes of habitat and prey usage. Second, we primarily observe high dissimilarity between closely  related species and repeated convergences among distantly related species in trait space, suggesting repeated trait diversification in both the crocodile icefish and notoperch radiations. Finally, we find only limited support for elevated diversification between sympatric species pairs. Instead, across each of our focal traits, we observe repeated convergence of sympatric species into overlapping regions of trait space, a pattern that is likely maintained by niche complementarity. Our findings are consistent with a scenario in which periodic environmental disturbances provided recurrent generation of ecological opportunities for diversification, thereby guiding the development of hypotheses of trait diversification in other environments characterized by a history of climate instability or heterogeneity.  : Evolutionary histories of the top three principal component (PC) axes of body shape variation demonstrate several instances of convergence of distantly related species in aspects of body shape as well as striking divergences among closely related species. A-C represent projections of the crocodile icefish phylogeny, and D-F represent projections of the notoperch phylogeny. In all panels, the x-axis represents time in millions of years, with branch lengths reflecting estimated divergence times. The y-axis reflects variation along each of the first three PC axes of body shape, with tip placement along the y-axis reflecting the PC score for that species along each PC axis.

Evidence for Correlated Evolution among Traits
Diversification along the benthic-pelagic axis of the water column is commonly associated with parallel divergences in fish body shape, where shifts from benthic to pelagic habitats are often accompanied by evolutionary transitions from deep-bodied forms to more slender fusiform body shapes (Rundle et al. 2000;Hulsey et al. 2013;Tavera et al. 2018). However, we find no evidence that repeated diversification along the benthic-pelagic axis of the water column has significantly influenced changes in body depth in either of our focal radiations (tables S5-S10). These counterintuitive findings about icefish and notoperch ecomorphological diversification are likely explained by the fact that many cryonotothenioid species are known to periodically forage outside of their position along the buoyancybased biotope axis (Casaux and Barrera-Oro 2013;Eastman 2020). This potential ecological plasticity could consequently reduce the expected tight correlations between body shape and niche use in icefishes and notoperches (table S5). Future work investigating the relationship between other axes of phenotypic diversification, such as fin shape or body width, and patterns of diversification in prey and habitat utilization in these radiations are warranted. Body size represents another important axis of phenotypic diversification in many evolutionary radiations, with a hypothesized positive relationship between body size and mean depth occupancy (Merrett and Haedrich 1997). However, similar to patterns observed for body shape, we find no correlations between notothenioid body size variation and disparity in either diet or depth (tables S7-S10). Instead, two species of icefish that lie at either extreme of variation in body size, Dacodraco hunteri and Channichthys rhinoceratus ( fig. 5B), occupy a nearly identical point in diet space, reflecting convergence on piscivory ( fig. 1A). Similarly, one of the smallest notoperch species (Trematomus scotti) and one of the largest notoperch species (Trematomus hansoni) each rely primarily on polychaetes as a prey resource ( fig. S2). One explanation for this lack of correlation is that many cryonotothenioid species appear to be opportunistic feeders, often pursuing the prey resources that are most abundant in a given season or locality (Hopkins 1987;Casaux and Barrera-Oro 2013). Moreover, the lack of correlation between body size and depth, diet, or buoyancy provides no evidence to suggest that patterns of habitat and resource utilization are constrained by body size. Taken together, our results suggest that the ecological plasticity that characterizes many cryonotothenioid species spanning a range of body sizes and shapes has likely been integral to the persistence and diversification of this clade as availability of prey resources and different depth habitats fluctuated during repeated glacial disturbances of the Antarctic continental shelf.

Repeated Ecomorphological Diversification in Icefishes and Notoperches
Trait disparity among closely related species is expected to decline as a result of niche filling (Simpson 1953  icefishes, sister species Champsocephalus esox and Champsocephalus gunnari similarly sit at opposite ends of the third PC axis of diet variation, with C. esox consuming mostly fish while C. gunnari feeds primarily on krill (figs. 1C, S1). These findings of high trait divergence among closely related species suggest that ecological and phenotypic diversification have occurred at recent timescales within both icefishes and notoperches.
Whether or not a lineage radiates in response to emergent ecological opportunities depends in large part on the ability of that lineage to readily adapt to the available niches (Stroud and Losos 2016). Across all of the studied traits, divergences among closely related species are occurring at very short timescales, suggesting a high level of trait lability in cryonotothenioids. For the key ecological trait of buoyancy, high trait lability is almost certainly the case. Recent work using CRISPR-Cas9 has experimentally validated that the modulation of a suite of genes associated with human degenerative bone diseases dramatically alters skeletal ossification and yields a range of skeletal phenotypes in notothenioids (Daane et al. 2019). This genetic mechanism may explain the complex multimodal depth distributions we observe for several species, suggesting a range of potential phenotypes for selection to act on during times of environmental change that allow lineages to take advantage of emergent ecological opportunities. Genomic evidence from island lineages has revealed suites of genes responsible for rapid evolution of body size and skeletal shape (Gray et al. 2015;Parmenter et al. 2016). It is possible that divergence along similar genetic pathways has enabled notothenioids to alter their buoyancy to take advantage of resources in newly assembling communities. Our results highlight that future genomic work focused on the evolution of the notothenioid bauplan, such as genome-wide association studies or notothenioid-guided genome editing in model organisms, is an incredibly promising and exciting research frontier.

Impacts of Range Overlap and Clade Age on Trait Divergence
Coexistence of closely related species in adaptive radiations is frequently attributed to character displacement (Tobias et al. 2014;Gillespie et al. 2020), which is expected to promote increased trait divergences among species pairs that occur in sympatry (Dayan and Simberloff 2005;Pfennig and Pfennig 2009). Alternatively, trait divergence could be better predicted by time since common ancestry than by range overlap, where higher trait divergences may be expected given a greater amount of evolutionary time afforded for accumulation of trait differences. However, we find limited evidence to suggest that either range overlap or clade age significantly predicts patterns of trait divergences within either crocodile icefishes or notoperches (tables S11, S12). We propose that the lack of strong evidence for increased trait divergence among sympatric species pairs is likely driven by repeated convergent evolution among distantly related species along axes of habitat and feeding niche use. Although this finding of rampant convergence in sympatry seems paradoxical, it could be explained by a pattern of niche complementarity. Species that have converged along one dimension of niche space may be sufficiently divergent along another niche dimension that the competitive interactions among them are minimized, thereby facilitating their coexistence in a given area (Schoener 1974;Werner 1977;Losos et al. 2003). For example, the sympatric icefish species Chionobathyscus dewitti and D. hunteri converge in buoyancy yet diverge in prey usage. Chionobathyscus dewitti consumes primarily cephalopods, while the diet of D. hunteri consists primarily of fish (figs. 1A-1C, S1). Additionally, the sympatrically distributed notoperch species Trematomus eulepidotus and Trematomus newnesi exhibit considerable convergence in prey usage ( fig. 1D-1F) yet appear to partition the water column: T. newnesi occurs primarily in shallow-water habitats (mean depth p 160 m), while T. eulepidotus typically occupies deeper-water habitats (mean depth p 389 m; fig. 3; table S1). These and other potential cases of niche complementarity suggest a potential role of interspecific interactions in generating ecomorphological disparity and may explain why convergences among species along the benthic-pelagic axis of the water column are not associated with parallel convergences in bathymetric distribution or prey use (tables S5-S10). Given the general lack of correlation between body shape and buoyancy, diet, or depth, future studies investigating how cryonotothenioids utilize distinct strategies to occupy the same feeding niche and achieve high niche lability will be critical to characterizing the mechanisms of diversification that have shaped this unusual adaptive radiation.

Conclusions
The results of this study reflect a dynamic history of nearshore Antarctic marine habitats, in which benthic communities are subjected to repeated catastrophic events at large geographic scales that decimate their fauna (Huybrechts 2002;Thatje et al. 2005Thatje et al. , 2008Near et al. 2012;Dornburg et al. 2017). The wake of these events has likely facilitated the repeated generation of ecological opportunities for diversification and may be considered analogous to an experiment in which the early stages of adaptive radiation are iteratively repeated (Near et al. 2012). In many ways, this background of widespread environmental devastation and recovery is not unique to the Antarctic. On the contrary, pulses of environmental change have occurred throughout the Cenozoic, not only promoting local extirpation or the extinction of species at large spatial scales but also generating novel ecological opportunities that are correlated with the evolution of phenotypic disparity (Lucek et al. 2018;Folk et al. 2019). As such, our finding that periodic environmental disturbance has promoted the rapid ecomorphological divergence among close relatives that, in turn, results in high levels of convergence among distant relatives may represent a general feature of evolutionary diversification in highly dynamic environments.