Molecular early burst associated with the diversification of birds at the K–Pg boundary

Complex patterns of genome and life-history evolution associated with the end-Cretaceous (K– Pg) mass extinction event limit our understanding of the early evolutionary history of crown group birds [1-9]. Here, we assess molecular heterogeneity across living birds using a technique enabling inferred sequence substitution models to transition across the history of a clade [10]. Our approach identifies distinct and contrasting regimes of molecular evolution across exons, introns, untranslated regions, and mitochondrial genomes. Up to fifteen shifts in the mode of avian molecular evolution map to rapidly diversifying clades near the Cretaceous-Palaeogene boundary, demonstrating a burst of genomic disparity early in the evolutionary history of crown birds [11-13]. Using simulation and machine learning techniques, we show that shifts in developmental mode [14] or adult body mass [4] best explain transitions in the mode of nucleotide substitution. These patterns are related, in turn, to macroevolutionary shifts in the allometric scaling relationship between basal metabolic rate and body mass [15, 16]. In agreement with theoretical predictions, this scaling relationship appears to have weakened across the end-Cretaceous transition. Overall, our study provides evidence that the Chicxulub bolide impact [17] triggered integrated patterns of evolution across avian genomes, physiology, and life history that structured the evolutionary potential of modern birds.


Introduction
It has been over forty years since Alvarez et al. [17] provided chemical evidence indicating that the Cretaceous-Paleogene (K-Pg) mass extinction was associated with an extraterrestrial impact. Subsequent research has confirmed and refined our understanding of this abrupt and cataclysmic event (e.g., [7,[18][19][20][21]). Mounting evidence suggests that the K-Pg extinction event triggered convergent patterns of life-history evolution [22]. For example, some lineages experienced a transient "Lilliput effect" in which average body sizes became smaller, likely through faunal sorting, dwarfing, or miniaturization [23][24][25][26][27]. While great effort has been devoted to investigating extinction patterns among various groups across the K-Pg boundary [17,21,28], the impact of the end-Cretaceous mass extinction event on the genomes of surviving lineages has received little attention.
Given that life-history parameters such as body mass, generation length, and metabolic rates have been linked to different aspects of molecular evolution [29,30], convergent patterns of life-history evolution are expected to leave an imprint of the K-Pg extinction event in the genomes of surviving lineages [4,6]. Only a few studies of animal taxa have attempted to directly investigate how the aftermath of the K-Pg mass extinction event shaped genome evolution (e.g., [4,8,[31][32][33][34]). For example, the evolution of polyploidy may be associated with the K-Pg transition in plants [35,36], and avian substitution rates may have increased due to extinction-related size-selectivity [4].
Many studies attempting to connect events in Earth's history to patterns of genome evolution rely on inferences from molecular clock analyses (e.g., [4,35]). These approaches can reveal heterogeneous patterns in the tempo of molecular evolution (e.g., [6,37,38]), but typically assume that the underlying data evolved according to the expectations of a uniform substitution model. If this assumption is violated, time-homogeneous models may obscure important evolutionary patterns. Nevertheless, techniques that enable substitution models to vary across a clade's evolutionary history have not yet seen widespread adoption in the macroevolution literature (e.g., [39][40][41][42]). Detecting model transition points on a phylogeny may provide evidence of evolutionary transitions in the "mode," or process that generated the observed data [43] -a concept that, when combined with hypothesis testing, has been termed "phylogenetic natural history" [15,44]. Thus, investigating patterns of model shifts in the context of genome and life-history evolution may reveal unknown links between Earth's history and evolutionary processes.
Here, we combine approaches from molecular systematics and phylogenetic comparative methods to investigate molecular model heterogeneity across the avian tree of life. We apply a stepwise approach to estimating where molecular substitution model parameters have shifted across a clade's evolutionary history, implemented in Janus [10]. Our analyses reveal wellsupported shifts in estimated equilibrium base frequencies across exons, introns, untranslated regions (UTRs), and mitochondrial genomes. Remarkably, model shifts are mostly constrained to previously hypothesized clade originations associated with the K-Pg boundary. This pattern indicates a temporal interval early in crown bird evolutionary history in which the mode of avian genomic sequence evolution rapidly diversified, consistent with an "early burst" macroevolutionary process (e.g., [12,13,43,45,46]). We examine how this pattern of molecular model shifts is related to macroevolutionary changes in avian life-history variation and focus on aspects of breeding ecology, development, senescence, and metabolism that are thought to have experienced intense selection or relaxation of evolutionary constraints during the K-Pg transition (e.g., [4-6, 47, 48]). Specifically, we investigate the hypothesis that molecular shifts indicate shifts in the adaptive landscape [49] or evolutionary allometries [15]. Our findings demonstrate a roadmap for linking molecular shifts inferred from the comparative analysis of genomic sequences with macroevolutionary patterns detected from the fossil record.

Molecular model shifts
We report details of dataset assembly as supplementary material. Fifteen molecular model shifts are identified on total-clades estimated to have originated near the K-Pg boundary [1, 3, 50, 51] ( Figure 1, Supplementary Figure 1, 7a-d). Considering overlap across data types, Janus detects thirteen phylogenetic regimes (one ancestral + twelve derived) that are required to explain the heterogeneity in modes of sequence evolution across exons, introns, UTRs, and mtDNAs, relative to the [51] MRL3 topology (Figure 1). Most molecular shifts are concordant with the origins of diverse ancient clades previously recognized at ordinal or superordinal taxonomic ranks. These include Notopaleognathae, Tinamiformes, the sister clade to Tinamiformes (in the MRL3 tree; an unnamed clade uniting Rheiformes, Casuariiformes, and Apterygiformes), Neognathae, Columbea, Passerea, Otidae (i.e., Otidimorphae + Strisores, sensu Wagler 1830, Jarvis et al. 2014), the sister clade to Otidae (in the MRL3 tree; the remainder of Neoaves), Aequornithes, Coraciimorphae, Psittaciformes, and Passeri (Supplementary Table 1 For every case, our approach identifies molecular shifts with 100% of the model weight when considering a shift's existence and location, indicating that shifts have a robust statistical signal. Several molecular shifts fall on the GC-AT axis (Figure 1, Supplementary Figure 7a-d) of compositional variation, with the most pronounced shifts occurring in exon data, followed by mtDNA, UTRs, and introns. We detect the highest frequency and magnitude of compositional shifts in our large exon dataset ( Figure 1). Nonetheless, there is no trend relating dataset size to the magnitude of inferred substitution parameters (e.g., the relatively small mtDNA data set indicates shifts with significant deviations in estimated equilibrium base frequencies). We observe similar deviations from empirical frequencies among coding (exon, mtDNA) when compared to non-coding sequences (intron, UTR) ( Figure 1, see Discussion). Thirteen phylogenetic regimes encompassing seventeen model shifts are required to explain heterogeneity in equilibrium base frequencies across genetic data types. Fifteen shifts are detected at nodes with stem ages within ~5 Ma of the K-Pg boundary [51]. Top: the aggregate signal of molecular model shifts across nuclear and mitochondrial data types identified by Janus, mapped onto the MRL3 supertree, with the ancestral regime "0" in black (*Otidae = Otidimorphae + Strisores, sensu Wagler 1830, Jarvis et al. 2014). Patterns of molecular model shifts across data types are summarized as 2x2 grids (Supplementary Figure 3, 7). Numeric labels at each grid position correspond to a particular shift in a specific data type. Pie charts summarize the positive detection rate ("P") for shifts in trait optima θ(t) across eight lifehistory traits relative to the null false positive ("FP") rate (e.g., ℓ1ou detection rate / (ℓ1ou detection rate + false positive rate), under AICc; see Supplementary Figure 4). Below: estimated magnitude of shifts in equilibrium base frequency relative to the empirical base frequencies for a given taxon partition for each data type, ordered by dataset size (also see Supplementary Figure 7a-d). Edges with wellsupported shifts in metabolic allometry are labeled with an asterisk, with the most substantial support observed for Coraciimorphae (pp = 98%, see results and Figure 3).

Life-history evolution
In general, analyses of life-history data support the hypothesis that molecular model shifts coincide with shifts in the evolutionary optima of life-history traits. In this context, optima [θ(t)] are equilibrium values within trait space that a lineage evolves towards under stabilizing selection and genetic drift. The detection rate for shifts in θ(t) occurring on candidate edges, normalized relative to the false positive rate under mvBM, is consistently high (Figure 1); all candidate cases are supported under AICc or pBIC criteria, but not both in a few cases (e.g., median of 76.2-87.2%, see Supplementary Appendix for additional detail). Analyses with a Random Forest machine learning approach are consistent with those from the model-based analyses reported above (further detail reported in the Supplementary Appendix). Avian developmental mode ("ChickPC1") [14], followed by adult body mass, are consistently the most important predictors of taxon partitions identified by molecular model shifts (AUC = 0.94). By contrast, traits reflecting substrate or dietary preferences are ranked as relatively unimportant, except for granivory, which is ranked fourth after average clutch size ( Within Neoaves, molecular model shifts are broadly associated with θ(t) shifts toward increased altriciality at hatching or decreased adult body mass relative to θancestral (7/7 and 6/7, respectively, Figure 2). Although with an overall lower θ(t) than θanc, Aequornithes and Psittaciformes show derived increases in body mass θ(t), along with derived shifts toward increased altriciality θ(t) (Figure 2, right). Outside Neoaves, developmental mode θ(t) within Palaeognathae is not clearly associated with molecular model shifts. For body mass, however, θ(t) for Tinamiformes is similar to θanc., while its unnamed sister clade shows a marked increase in θ(t) relative to θanc. (Figure 2). Notably, an alternative set of analyses estimating θ(t) separately for Struthio+root indicates all molecular shifts, including those within Palaeognathae, are associated with derived decreases in body mass or ChickPC1 θ(t) (Supplementary Figure 6 and supplemental discussion). These patterns support the hypothesis that molecular model shifts indicate evolutionary shifts in life-history θ(t). Left: permutation-based variable importance for life-history traits reported in [52] and [14]. With a Random Forest classifier, we identified variation in avian developmental mode ("ChickPC1" from [14]) and adult body mass as closely associated with taxon partitions recognized by Janus in an analysis of exon data (Supplementary Appendix). Right: estimates of trait optima θ(t) when anchored with molecular model shifts from nuclear genetic data (100 parametric bootstraps; colors and labels match Figure 1). Background distributions (light gray) indicate expected values of θ(t) at T = 0. For reference, vertical lines mark θancestral and phylogenetic levels within Neoaves (e.g., demarking lineages with probable pre-and post-Cretaceous originations). Molecular model shifts are often associated with shifts toward increased altriciality at hatching or decreased adult body mass (also see Supplementary Figures 5, 6).

Metabolic allometry
Bayesian models of metabolic allometry indicate that deviations from a prior expectation of ¾ power law scaling are associated with molecular model shifts close to the K-Pg boundary. Modal estimates for slope (βmass) range from 0.68 (~⅔) to 0.84 (~⅘), and intercept (β0) from -4.32 to -3.17 ( Figure 3, Supplementary Table 1). These values are similar to those previously estimated for avian and mammalian subclades [53][54][55]. Compared to life-history traits like mass or developmental mode ( Figure 2, Supplementary Figure 5, 6), the evolution of metabolic scaling parameters appears more uncertain, with broad posterior estimates in some cases (e.g., Palaeognathae, Figure 3). Nevertheless, these analyses reveal that the origins of several K-Pgassociated subclades within Neoaves (e.g., Passerea) coincide with a shift toward overall lower body mass, as well as lower slope and higher intercept terms (e.g., under 10 kg as noted in [53,56]). However, shifts in the allometric optima of individual candidate subclades are not always in a consistent direction ( Figure 3).
Seven edges detected at a 10% posterior probability cutoff reflect K-Pg-associated clade originations ( Figure 3, Supplementary Table 1). Under a more conservative threshold, only three candidate edges are detected with moderate to strong support, including the unnamed sister clade of Otidae (pp~38%), Columbea (pp~39%), and Coraciimorphae (pp~98%) (Figure 1). The diverse clade Coraciimorphae is the only candidate group for which molecular model shifts are detected across all nuclear genetic data types. Overall, metabolic parameter estimates are consistent with the hypothesis that allometric shifts in avian metabolism may be associated with molecular model shifts.  Table 1). Median parameter estimates are indicated with white dots. Horizontal dashed lines indicate the prior mean (prior density intervals are shown on the right vertical axis for reference); dotted lines indicate the range of median parameter estimates. A shift toward lower slope and higher intercept characterizes lineages with probable post-Cretaceous origin (e.g., Passerea).

Discussion
Unraveling interactions among significant events in Earth's history and macroevolutionary patterns is a fundamental and persistent problem in evolutionary biology. Our study investigates patterns of sequence evolution and the adaptive radiation of birds near the K-Pg boundary. We demonstrate how proximity to the K-Pg boundary increases the probability of molecular model shifts (Figure 1, Supplementary Figure 1a-d), linking a major mass extinction to the macroevolution of the avian genome. By anchoring a series of phylogenetic comparative models with molecular model shifts, we find evidence that shifts in the mode of genomic evolution were likely concurrent with shifts in the evolutionary optima θ(t) of important avian life-history traits (Figure 2), as well as shifts in metabolic allometric slope βmass and intercept β. Broadly, model shifts in genomic sequences identify shifts toward increased altriciality or smaller adult body mass, consistent with the hypothesis of a K-Pg-associated "Lilliput effect" [4].
Our examination of metabolic allometry provides insight into the biological consequences of size evolution; transitions toward smaller sizes are correlated with weaker scaling relationships between metabolic rate and body mass (e.g., along the Neoaves-Passerea backbone; Figure 3, [4]). This pattern implies that the energetic costs of evolutionary increases in size are reduced in clades with a smaller average body mass. In a scenario like the K-Pg extinction, during which networks of ecological competitors were reset [57], the survivorship of clades with smaller body sizes-and weaker associations between metabolic rate and body mass-may have thus facilitated the evolution of a range of physiological strategies in the early Cenozoic (e.g., [58][59][60]). This deduction is consistent with theoretical and empirical advances which predict transitions toward harsher environments with increased extrinsic mortality drive the evolution of lower metabolic scaling exponents because of selection to maximize lifetime reproduction [16,61]. These associated phenomena lead to earlier maturation and faster growth [16], aligning with our inference of increased altriciality associated with the K-Pg extinction (e.g., [14,62,63]).

Recognizing early bursts
Lineages can enter novel adaptive zones during diversification due to ecological, geographic, or phenotypic opportunities [43,64]. The aftermath of mass extinctions, especially those of short duration, may present all three classes of opportunities, resulting in recovery faunas that experience "early bursts" of lineage and character diversification [45,65]. If opportunities for diversification become more constrained as ecospace fills, rates of morphological evolution and lineage accumulation should decline, with the fastest rates of change restricted to a short interval following the mass extinction event [11,12,66]. Accordingly, we expect initially high rates of evolution to generate outsized disparity early in post-extinction adaptive radiations.
An exclusive focus on rates of change, however, (e.g., [13]) precludes recognition of early burst patterns based on other sources of disparity [67]. Our approach diagnoses a molecular early burst of disparate processes of nucleotide sequence evolution. It is, therefore, conceptually more similar to paleontological approaches examining patterns of disparity [46,67,68] than it is to comparative techniques estimating rates of change in quantitative or molecular characters [4,13,69]. While many studies have quantified early burst phenomena through patterns of morphological evolution or rates of lineage diversification, we show that ancient diversification events may impart a signature of genomic disparity that remains detectable in surviving lineages for tens of millions of years. In this context, the molecular model shifts we identify reveal "genomic fossils" that describe canalized macroevolutionary regimes, within which sequence evolution is aligned along contrasting axes of substitutional variation.

A novel dimension of avian adaptive radiation
The detection of numerous model shifts within a ~5 Ma interval of the K-Pg boundary ( Figure  1) is likely not a coincidence -and instead reveals a canonical "early burst" pattern in which ecological disparity and lineage disparity accumulated rapidly in the early history of crown birds [45,70]. Given that the uncertainty in estimated molecular divergence times typically exceeds 5 Ma [1, 3, 4, 50, 51, 71], a conservative interpretation of available divergence time estimates does not reject the hypothesis that these events were closely linked (Figure 1, Supplementary Figure  1). We emphasize that Janus does not have access to information about the absolute timing of divergence events; it is, therefore, striking that shifts cluster in temporal proximity on a welljustified time-calibrated phylogeny [51] (Figure 1, Supplementary Figure 1).
Consistent with this pattern, models of quantitative trait evolution estimate short phylogenetic half-lives (t1⁄2 = ln(2)/α, e.g., t1⁄2 body mass ~ 0.18 [0.17-0.52] Ma, Figure 2). Such short intervals imply a scenario with rapid character displacement followed by a relatively stationary process (assuming a median generation length of three years [52] reflects a t1⁄2 of ~ 60,000 generations). This interpretation is prompted by a conservative view of the fossil record, which currently only supports limited crown bird diversification in the Late Cretaceous (e.g., [4,6,50]), and short circum-K-Pg internodes. Our results, therefore, support the hypothesis that developmental and life-history traits were canalized early in crown bird evolutionary history [2,14,72], partitioning avian bauplan and higher-taxa (e.g., [67,73,74]). Short t1⁄2 may also partly explain deviations between estimated equilibrium and empirical base frequencies in these data ( Figure 1). Considering the coincidence of molecular model shifts with shifts in the evolutionary optima of life-history traits (e.g., Supplementary Figure 4), both phenomena likely reflect integrated evolutionary responses to the post-Cretaceous adaptive landscape and a permanent shift to new adaptive zones (e.g., [75]).
While many mechanisms link substitution rates to the life-history spectrum [30], we have limited intuition about how the mode of molecular evolution may relate to life-history variation. One idea proposed to explain variation in DNA compositional heterogeneity links a recombination process known as GC-biased gene conversion (gBGC) to generation length, effective population size, and co-varying life-history traits [10,[76][77][78][79][80][81]. Increases in Ne [8] and decreased generation lengths [4] could plausibly explain some of the patterns we observe (e.g., higher Ne is predicted to increase the efficacy of gBGC [79]). The fact that we observe many substitution model shifts in exon data is also consistent with the hypothesis that transcriptionally active regions experiencing more recombination may be more subject to gBGC [82,83]. Lastly, our observation that coding regions have considerably greater deviations between estimated equilibrium and empirical base frequencies could be due to functional constraints, codon usage bias, recombination, or selection (see supplementary analysis of codon usage and Supplementary Figure 7a-d).
Although we do not investigate it directly, our results imply a mechanism to explain the "data type" effect in studies of avian phylogeny [1, 71,84,85], in which phylogenetic analyses of coding or non-coding nucleotide data recover conflicting signals of phylogeny. Timeheterogeneous processes of molecular evolution are expected to introduce systematic error in attempts to resolve phylogeny through model-misspecification (e.g., [6,84,86]). We note that the data type exhibiting the weakest deviations between empirical and estimated equilibrium base frequencies are non-coding introns (Figure 1). Thus, our results support the conclusions of [71,84], who favor non-coding data in the inference of avian phylogeny. At the same time, our results show that non-coding data are likely not immune from these issues. Given that numerous model shifts are detected for both coding and non-coding sequences in a region of the avian phylogeny with uncertain topological relationships [71,87], future progress on avian phylogeny will hinge on our ability to rule out or better accommodate these potential sources of error.

Conclusion
Although high-throughput sequencing has clarified the history of many vertebrate clades, early evolutionary events within ancient lineages of crown birds-a group comprising over 10,000 extant species-have long been shrouded in mystery. Our study illuminates how the diversification of modern birds in close association with the end-Cretaceous mass extinction was characterized by rapid shifts across several axes of life-history variation. Directional selection on these parameters across the end-Cretaceous mass extinction event, such as selection towards increased altriciality or decreased adult body mass, appears to have shifted the processes of genome evolution. Ultimately, our results reveal a novel dimension of how one of the most significant events in Earth's history -the Chicxulub bolide impact -structured the evolutionary potential of surviving lineages and gave rise to the spectacular diversity of birds.

Nuclear Sequence Data Collection and Processing
We re-assembled an existing short-read sequence targeting 394 gene regions across 198 bird species and two crocodilian outgroups from Prum et al. [3]. These data were initially collected using target-capture of anchored hybrid enrichment (AHE) loci [88], a set of single-copy regions semi-conserved across vertebrates. We analyzed the existing raw sequencing reads with a common pipeline designed to extract phased exons. First, we removed low-quality regions and adaptor sequences using Trimmomatic v0.36 [89] and merged overlapping reads using FLASH v1.2.11 [90]. We assembled reads for each sample using Trinity v2.11 [91]. We then annotated assemblies by comparing assembled contigs to target loci using blat v36x2 [92]. To ensure we annotated orthologs, we retained only contigs with a reciprocal best-hit match to a target locus. To identify intron-exon boundaries, we used exonerate 2.4.0 to compare the nucleotide sequences of annotated loci to the protein sequences for the exons of each locus, based on protein-coding data and annotations from the zebra finch (Taeniopygia guttata, genome assembly bTaeGut1_v1). This approach assumes that intron-exon boundaries are conserved across the avian radiation. Occasionally, mapping of the exon sequence to the nucleotide sequence was discontinuous, suggesting the presence of an intervening non-coding region. In such cases, we retained the highest-scoring contiguous stretch of sequence only.
To identify variable sites, we mapped cleaned reads back to annotated contigs using bwa v0.7.17-r1188 [93] and used GATK v4.1.8 to mark duplicates [94]. We called variants on this alignment using GATK HaplotypeCaller and filtered it to only retain variants with coverage >20x and quality >20. Using this high-quality variant set, we recalibrated the base quality scores in the alignment files using GATK. We then called variants and phased them using HaplotypeCaller. Finally, we exported diplotypes and phased haplotypes per intron and exon in a coding region, masking any sites with coverage <2✕. Ultimately, we captured 453 exons, 573 introns, and 213 untranslated regions.
Before alignment, we applied a series of sequential filtering steps to remove remaining short or low-quality fragments. We removed (1) leading and trailing N characters from each fragment, and resulting sequences that were zero length (see locus-filtering.R script), (2) fragments with > 40% N characters, (3) fragments that were < 50 bp long, and (4) whole loci that lacked coverage for at least 10% of the taxa in the dataset.
We aligned phased exon sequences with the Multiple Alignment of Coding Sequences (MACSE) ALigning, Filtering, and eXporting pipeline (ALFIX) [95]. MACSE-ALFIX chains together several programs that perform reading frame aware alignment with MACSE and subsequent alignment filtering with HmmCleaner [96] to remove non-homologous sequence fragments. We aligned phased non-coding sequences with Fast Statistical Alignment (FSA) [97]. We calculated alignment statistics using AMAS [98]. We used trimAl [99] to evaluate the effect of 5%-30% alignment column occupancy filtering on alignment length and the loss of parsimony informative sites. We ultimately filtered our non-coding alignments to require a minimum column occupancy of 5% (i.e., 95% of the sequences in an alignment are allowed to contain a gap for a given site). This procedure, which we believe is conservative, increased the signal-tonoise ratio in these data by removing stretches of unaligned nucleotides (characteristic of FSA alignments) while retaining most of the informative data [100] (also see [101]). Unfiltered alignments and the final filtered dataset are provided as Supplementary Data.

Mitochondrial Sequence Data Collection and Processing
We ran Mitofinder 1.4 [102] to identify the mitochondrial regions from the previously assembled contigs. For reference mitogenomes, we used complete mitogenomes available in GenBank (Supplementary Table 3). When available, we used a reference from the same order (though for Passeriformes, we used different references for oscines [Passeri] and suboscines [Tyranni]); in a few cases, it was necessary to use a reference from a closely related order (Supplementary Table  3). We then extracted the 13 protein-coding genes and 2 rRNAs from the mitofinder output (final_genes.fasta file). In some cases, limited mitochondrial data were recovered (Supplementary Table 3). In those cases, we searched GenBank for the same or a phylogenetically equivalent species that could be substituted. When no suitable alternative was available from GenBank, we also used mitogenomes assembled from the raw data collected as part of other studies ( [103][104][105][106][107], and Braun et al. in prep). To increase data coverage in five cases, we generated chimeric sequences using available GenBank data from multiple individuals of the same species (Supplementary Table 3).
Once a set of sequences had been assembled, we performed an initial analysis using these data, combined with a larger set of mitogenomes (Kimball et al., in prep) to ensure sequences were correctly identified (placed phylogenetically with expected relatives) and did not exhibit unusually long branch lengths, which might suggest assembly errors. To do this, we ran an initial alignment using MAFFT 7.294b [108] using default parameters. This alignment was then analyzed in IQ-TREE 2.1.2 ([109]) using the GTR+I+G4 substitution model with 1000 ultrafast bootstrap replicates ( [110]). Lastly, we regenerated alignments for the present study, using the same procedure described above for nuclear coding and non-coding data.

Phylogenetic framework
To avoid issues of circularity related to inferring molecular patterns and phylogenetic topology from the same molecular dataset and to control for stochastic resolutions of Neoaves (e.g., [87]), our focal analyses use the MRL3 supertree [51] as a topological constraint. This topology balances the signal of phylogeny among several recently published avian genomics datasets, and resolves the seven major higher-level clades identified by [84], as well as most intra-ordinal clades that are robustly supported by [3]. It is also in line with a growing number of studies that have suggested that early diversification events within the avian crown group were associated with the K-Pg boundary [1, 3,50,51]. As inference of avian phylogeny is an active area of research [71], we explored how patterns of gene-tree discordance (e.g., [1]) may confound inference of molecular model shifts or potential statistical associations with the K-Pg boundary (Supplementary Appendix, Supplementary Figure 1b).
Molecular model shift analysis with Janus takes a rooted input phylogram with branch lengths in substitution units (method details reported in Supplementary Appendix, [10]). To generate starting trees for Janus, we used our reprocessed datasets and estimated maximum likelihood branch lengths with IQ-TREE v 2.1.1 [109,111]. For each data type, we applied an optimal partition model selected with the MFP+MERGE approach in IQ-TREE [112,113], with each locus defined as the unit for partitioning. We estimated molecular branch lengths separately for exons, introns, untranslated regions, and mtDNAs, but kept the topology fixed across datasets.
The Janus algorithm fits models that describe patterns of substitutions irrespective of the absolute timing of divergence events (Supplementary Appendix). Therefore, temporal patterns must be evaluated on a reference timeline. To interpret our model-shift results on a time-calibrated phylogeny, we used congruification [114] with treePL [115] to apply the divergence date estimates from the reduced taxon set analysis presented in [51] to the phylogenetic branch length estimates derived from the present study. The well-constrained divergence estimates from [51] are broadly congruent with those reported across several phylogenomic analyses of independent datasets [1, 3,50]. These estimates reject the hypothesis that many modern avian clades originated in the Cretaceous and centers the diversification of most superordinal variation within ~ +/-5 Ma of the K-Pg boundary (Figure 1). Therefore, our interpretations are conditional on this general divergence time scenario, an area of active research [4,6,50] (see discussion).

Fitting time-heterogeneous models with Janus
For nuclear genomic data, we considered the signal across three concatenated datasets (exons, introns, and UTRs). For mitochondrial data, we considered three alternative datasets (all data combined, protein-coding genes combined, and rRNAs combined). Our focus on data type mirrors recent developments implicating this axis of genomic variation as a primary source of phylogenetic incongruence [1, 71,84]. We fit time-heterogenous substitution models to each dataset with Janus (Supplementary Appendix, commit 8952e31d, https://git.sr.ht/~hms/janus). By default, Janus optimizes a model which enables shifts in equilibrium base frequencies across a provided phylogram. We set each search to accommodate rate heterogeneity across sites according to a discretized gamma distribution (-g) and to assess model weights for the existence (-u) and location (-l) of model shifts. Simulations indicate that this combination of parameter options has high power (e.g., a negligible false-positive rate) to detect the phylogenetic position of molecular model shifts (Supplementary Appendix, [10]). Considering our genetic dataset's taxonomic sample, we set the minimum clade size to >|= 4 (-m 4). Thus, the set of possible shift configurations reflects 101 internal nodes spanning ~77 Ma (square markers in Figure 1; e.g., with postorder traversal, any node (excluding the root) with >|= 4 descendant edges).

Life-history and metabolism datasets
To assess how the configuration of molecular model shifts detected with Janus may be related to life-history variation, we considered how life history might vary across multiple dimensions (e.g., [4,116]). We assembled two life-history datasets to minimize the amount of missing data in each analysis. The first dataset focused on quantitative life-history traits and was compiled from the [52]. These data include body mass, modeled generation length, latitude centroid, mean clutch size, annual adult survival, age at first breeding, maximum longevity, and categorically coded variables for diet, habitat, and diurnality and migratory status. We also included a metric of avian developmental mode ("ChickPC1") that describes variation in hatchling state along an altricial to precocial spectrum [14]. These data reflect exact species matches relative to those in the re-assembled nuclear genetic dataset.
The second dataset reflects energetic constraints on life-history variation and includes basal metabolic rates (BMR) expressed in watts and associated body masses. Metabolic rates broadly scale as a ~¾ power law function of organism mass and reflect rates of energy flow in and through organisms [16,54,[117][118][119][120]. [15] previously considered the hypothesis that allometric scaling parameters relating BMR and body mass have evolved across the vertebrate tree of life. We apply the same general approach to our sample of avian metabolic diversity (below). We first collected available BMR records from the AnAge senescence database Build 14 [121]. For most of the exact species in the present dataset (and most avian species in general), conspecific BMR data have not been measured. Therefore, we conducted an extensive literature search for each avian family in the molecular dataset and filled in many missing entries by identifying phylogenetically equivalent matches (e.g., at the genus level) for which BMR and mass data were available.
Several downstream analyses required complete datasets, so we used two methods to generate unbiased estimates of missing values under a multivariate Brownian motion process (mvBM). In the case of the larger eight-dimensional breeding ecology dataset, we used Rphylopars [122] to fit a variance-covariance matrix (VCV) and to estimate values for missing entries. In the case of the two-dimensional metabolic scaling dataset, we used mvMORPH [123] to compare the fit of alternative multi-regime, mvBM models based on the model shift points identified by Janus. In the latter case, the values of the imputed data were virtually identical across alternative models (e.g., R 2 > 0.98), so we selected the model with the lowest AIC score to use for downstream analyses.

Analysis of life-history data
Using multiple approaches, we investigated the degree to which patterns of life-history variation reflect distinct evolutionary regimes that coincide with molecular model shifts. Several methods have been developed to automatically generate evolutionary hypotheses by identifying an optimized configuration of evolutionary models describing variation in the process of trait evolution (e.g., [124,125]), but few are expressly multivariate (e.g., [49,126]). We investigated model heterogeneity across our high-dimensional life-history dataset with the bootstrapping approach implemented in the software ℓ1ou [49]. ℓ1ou uses a phylogenetic lasso method to identify points on a phylogeny where a trait's optimum value θ(t) has shifted, assuming α (the "pull" toward the optimum or adaptation rate) and σ 2 (the Brownian diffusion rate parameter) are fixed across the tree. The ℓ1ou approach is extended to multiple traits by assuming that traits shift their optimum simultaneously and in the same location on the tree [49].
Conveniently, ℓ1ou allows the researcher to specify a set of candidate edges for the lasso approach to consider for shifts in θ(t). This attribute allows us to articulate the specific models we want to compare. We ran ℓ1ou with a constrained set of candidate edges reflecting the 12 candidate shift edges identified across analyses of different molecular data types (Figure 1,  Supplementary Table 1). Thus, for a given ℓ1ou analysis of this type, ℓ1ou can infer 0-12 shifts in θ(t). We repeated this procedure with the AICc and pBIC information criteria, as recommended by [49], and used 100 bootstrap replicates to assess the positive detection rate for each candidate edge. We emphasize that our intention is not to identify every case where a lifehistory shift may have occurred across avian phylogeny; our goal is to assess how much statistical support exists for shifts in life-history trait optima that coincide with shifts identified in our analysis of molecular data.
We validated our results through comparison to a null distribution of shift detections reflecting the false positive rate under multivariate Brownian motion (e.g., without shifts in θ(t).
Using the eight-dimensional variance-covariance matrix (VCV) estimated by RPhylopars [122], we simulated 500 null datasets using the function simRatematrix in the R package ratematrix [127]. We then analyzed each simulated dataset with ℓ1ou as previously specified. For each candidate edge, we used Fisher's exact test [128] to assess whether the frequency of positive shift detections in the empirical dataset was significantly greater (one-tailed p-value = 0.05) than the null false-positive rate observed across simulated datasets. Tables of p-values and odds ratios are reported as supplementary material (Supplementary Figure 4, Supplementary Table 2).
Lastly, we investigated which life-history traits are most closely associated with molecular model shifts using a machine learning approach implemented in the tidymodels framework [129] (Supplementary Appendix). To explore these results, we fit fixed shift OUM (shifting θ(t), with fixed α and σ 2 ; equivalent to that used by ℓ1ou) models using OUwie [130]. We used 100 parametric bootstrap replicates to estimate model parameter uncertainty. We also simulated 1000 datasets under each fitted OUM model to visualize the expected distribution of trait values at T (time) = 0 ( Figure 2). These models, therefore, reflect the expected shifts in trait optima which coincide with bipartitions identified in analyses of molecular model shifts.

Analysis of metabolic rate data
To assess whether molecular model shifts are associated with shifts in patterns of metabolic scaling, we assessed support for coincident shifts in metabolic allometry. We utilized the Bayesian phylogenetic framework implemented in the R package bayou 2.0 [15,124]. bayou applies a reversible jump Markov Chain Monte Carlo approach to detect the magnitude, number, and phylogenetic position of model shifts. Using bayou, we implemented an allometric regression model which relates BMR and body mass logarithmically, and for which slope βmass and intercept β0 evolve under a multi-regime Ornstein Uhlenbeck (OU) process. Here, model shifts reflect shifts in the optimum of the evolutionary allometry between BMR and body mass.
Using the rjMCMC approach in bayou, we estimated the posterior probability of an allometric shift occurring along the 12 candidate edges identified by Janus. Under a Poisson prior, we specified the mean number of shifts across the phylogeny reflecting 2% of the total edges in the tree (λ = 8) with equal probability. In this context, maximal posterior probability indicates an increase over the prior probability by ~50%. We ran each analysis across three replicate chains for 10 million iterations, sampling every 1000 iterations. Given that we did not have consistent estimates of measurement error, we followed the approach of [15] and explored alternative analyses assuming a Standard Error of 0.1 or 0.01 for BMR and recovered a negligible impact (not shown).
Our priors for α, 2 , βmass, β0 reflect half-Cauchy or Gaussian expectations: Replicate analyses with rjMCMC identified numerous shifts in the slope and intercept at a posterior probability cutoff of 0.1, after discarding the first 40% of samples as burn-in. Following rjMCMC runs, we re-estimated model parameters on a fixed configuration model reflecting all molecular shifts identified with Janus ( Figure 3). We assessed model convergence by examining Gelman and Rubin's R statistic [131] and effective sample sizes across chains and parameters.

Data Availability
All data or software code generated for the present manuscript is available at the corresponding author's GitHub repository https://github.com/jakeberv/avian_molecular_shifts and archived at Zenodo pending peer review. Janus has been implemented in Golang and C and is available at https://git.sr.ht/~hms/janus and https://git.sr.ht/~hms/hringhorni. the authors have applied a Creative Commons Attribution (CC BY) license to any Author