## Abstract

Lobe-fins transformed into limbs during the Devonian period, facilitating the water-to-land transition in tetrapods. We traced the evolution of well-articulated skeletons across the fins-to-limbs transition, using a network-based approach to quantify and compare topological features of fins and limbs. We show that the topological arrangement of bones in pectoral and pelvic appendages evolved in parallel during the fins-to-limbs transition, occupying overlapping regions of the morphospace, following a directional trend, and decreasing their disparity over time. We identify the presence of digits as the morphological novelty triggering topological changes that discriminated limbs from fins. The origin of digits caused an evolutionary shift toward appendages that were less densely and heterogeneously connected, but more assortative and modular. Disparity likewise decreased for both appendages, more markedly until a time concomitant with the earliest-known tetrapod tracks. Last, we rejected the presence of a pectoral-pelvic similarity bottleneck at the origin of tetrapods.

## INTRODUCTION

The evolution of tetrapod limbs from fish fins is heralded as one of the most important vertebrate morphological and functional transitions (*1*–*8*). Establishing what makes an appendage a fin or a limb is key to properly characterizing the fins-to-limbs transition (*3*). Functional criteria are of limited use because of the general consensus that limbs first evolved to move underwater (*5*, *9*). Developmental and paleontological studies place the distinction between fins and limbs in the most distal region, which bears the carpals/tarsals and digits in limbs and the radials and dermal lepidotrichia in fins (*3*, *10*). The distinction between fins and limbs blurs when we look at the lobe-fins of transitional tetrapodomorphs, such as *Panderichthys* or *Tiktaalik* (*11*–*14*). Both sarcopterygian fins and limbs share a division of the appendicular skeleton into three endoskeletal domains (*15*), of which the most distal domain differs most between sarcopterygian fishes (i.e., branching radial bones) and tetrapods (i.e., autopod with a mesopod and digits). Although in the past, researchers have disagreed about whether a zeugopod-mesopod boundary (wrist/ankle) (*16*, *17*) or the presence of digits alone (*3*) is sufficient to define limbs, the current convention is to define “true” limbs as appendages with digits (*6*). Although the anatomical organization of the distal radials and the autopod is superficially similar (i.e., a series of skeletal elements joined proximodistally) (*3*, *18*) and they share a common genetic control or “deep homology” (*19*), their anatomical similarity has never been assessed quantitatively. Moreover, pectoral and pelvic lobe-fins evolved into limbs in tandem during the fins-to-limbs transition due to the recruitment of a common developmental genetic toolkit (*2*, *20*). Because pectoral and pelvic appendages were originally different in their anatomy (*21*)—and still differ in the genetics of their girdles’ development (*22*)—we would expect to see a mix of evolutionary parallelisms/convergences (homoplasy) and divergences, as shared and specific developmental programs and biomechanical functions intertwined with each other during the fins-to-limb transition. Such a mix might result from compromises between these “evo-devo” and “evo-biomechanical” constraints.

As appendages evolved, the degree of similarity between the pelvic and pectoral appendages varied through time. Various authors have proposed alternative bottlenecks during evolution for the pectoral-pelvic similarity [reviewed in (*7*, *23*)]; these evolutionary bottlenecks represent times when pectoral and pelvic appendages showed the greatest anatomical similarity to each other (i.e., their morphologies showed less disparity). On the basis of skeletal and muscular anatomical and developmental features (*21*, *24*–*26*), pectoral-pelvic similarity bottlenecks have been proposed for the origins of ray-finned fishes, coelacanths, tetrapodomorphs, and tetrapods. In a recent study comparing the musculoskeletal network-anatomy of whole appendages (*23*), we found evidence for a pectoral-pelvic similarity bottleneck at the origin of sarcopterygians [as proposed in (*3*, *7*, *21*, *26*)] but not at the origin of tetrapods (as these same studies proposed). However, our previous work focused only on the network-anatomy of extant taxa (*23*). To further test the presence of a pectoral-pelvic similarity bottleneck at the origin of tetrapods for the network-anatomy of the skeleton, here, we have analyzed a broader sample of extinct sarcopterygian fishes and early tetrapods across the fins-to-limbs transition, including *Sauripterus*, *Sterropterygion*, *Gogonasus*, *Eusthenopteron*, *Cabonnichthys*, *Mandageria*, *Panderichthys*, *Tiktaalik*, *Acanthostega*, *Ichthyostega*, *Tulerpeton*, *Eryops*, *Balanerpeton*, *Celtedens*, *Seymouria*, *Westlothiana*, *Pantylus*, and *Hyloplesion*, for which the fully articulated pectoral and/or pelvic anatomy is reasonably well known. A prediction of the hypothesis of a pectoral-pelvic similarity bottleneck at the origin of tetrapods is that taxa closer to the split of Tetrapoda within sarcopterygians—where the bottleneck is—will have a greater pectoral-pelvic similarity (or lower disparity) than taxa that are farther away from the bottleneck.

To compare the skeletal anatomy of appendages in extinct and extant forms across the fins-to-limbs transition, and to better characterize anatomical parallelism/convergence and divergence between pectoral and pelvic appendages, we focused our analysis on their anatomical organization or network-anatomy, that is, the topological arrangement/pattern of skeletal elements of fully articulated appendages (Fig. 1A). This level of abstraction allowed us to compare evolutionary changes that are not amenable to quantification using other morphometric methods due to the large disparity of forms and presence/absence of parts between pectoral and pelvic fins and limbs (*23*, *27*, *28*). Furthermore, the abstraction retains biological meaning in that the contacts between skeletal elements reflect potential direct developmental and biomechanical interactions, for example, ontogenetic sequences of ossification or embryonic interaction, and joint reaction forces or ranges of motion. Using a network-based approach (*23*, *29*, *30*), we modeled the skeleton of fully articulated appendages as networks, in which nodes code for bones and links code for physical contact in a standardized resting pose. We compared the evolution of eight network-based topological variables (see Methods for details) in a phylogenetic context (Fig. 1B) to test whether (i) there are topological differences between fins and limbs and between pectoral and pelvic appendages, (ii) pectoral and pelvic anatomy followed convergent/parallel or divergent modes of evolution during the fins-to-limbs transition, and (iii) there was an evolutionary bottleneck in pectoral-pelvic similarity in tetrapods. We tested these hypotheses by comparing the occupation of appendicular morphospace, estimating shift(s) of evolutionary regimes, describing the evolution of disparity through time (DTT), and testing bottlenecks with phylogenetic regressions.

## RESULTS

### Topological discrimination of appendages

The network-anatomy of the appendicular skeleton varied for each taxon and between pectoral and pelvic regions (Table 1). We used a principal components analysis (PCA) to visualize global patterns of topological variance across anatomical networks and a permutational multivariate analysis of variance (PERMANOVA) to test for differences between fins and limbs, pectoral and pelvic appendages, and extinct and extant species.

The first two PCA components explained 79.3% of the total topological variation among appendages (Fig. 2A). The first axis of variation (60.1%) broadly discriminated between (i) more modular (higher P) sparsely connected (lower D) appendages, such as limbs, and (ii) less modular (lower P) densely packed (higher D) appendages, such as lobe-fins. The second axis of variation (19.2%) broadly discriminated between (i) more regular appendages, in which bones tend to have the same number of articulations (lower H) and which contact bones with a similar number of articulations (high A), such as limbs and anatomically plesiomorphic lobe-fins, and (ii) more heterogeneous appendages, in which bones have a varying number of articulations (higher H) and which preferentially contact bones with a different number of articulations (lower A), such as in the anatomically derived lobe-fins of *Neoceratodus*.

The PERMANOVA showed a statistically significant difference in topological variability between fins and limbs (*F*_{1,44} = 15.77, *P* = 9.9 × 10^{−5}; Fig. 2B). The assortativity of appendages (i.e., tendency of bones to contact bones with a similar number of connections) was the main discriminator between fins and limbs, which occupied opposite positions along the assortativity-axis of variation. Limbs had larger and positive values, whereas fins had lower and negative values. We can explain this difference by the presence of the autopod in limbs and, more specifically, by the presence of digits. Phalanges and, to a lesser extent, carpal/tarsal bones tend to articulate with other autopodial bones with a similar number of articulations. For example, phalanges connect in a proximodistal series to other phalanges or metacarpals/metatarsals so that most phalanges have two articulations (one proximal and one distal) to other phalanges that also have two articulations (hence, A increases). A similar pattern may occur among carpal/tarsal bones because of their nearly polygonal shapes. However, the PERMANOVA showed no significant difference in topological variability between pectoral and pelvic appendages (*F*_{1,44} = 0.28, *P* = 0.88; Fig. 2C) and between extinct and extant species (*F*_{1,44} = 2.24, *P* = 0.09; Fig. 2D). Pectoral and pelvic appendages did not occupy different areas of the morphospace, which indicates that they share a similar topological organization. Last, although extinct and extant taxa are statistically indistinguishable, they appear to occupy slightly different areas of morphospace: Extant species varied equally along PC1 and PC2 axes, while the variance in extinct species is primarily concentrated along the PC1 axis and barely varied along the PC2 axis.

### Evolution of pectoral and pelvic appendages

We assessed the potential for parallel/convergent and divergent changes in topology for pectoral and pelvic appendages by estimating shifts in evolutionary regimes (SURFACE) and analyzing DTT. The SURFACE analysis on PC1 and PC2 estimated a shift in mean values at the root branch of sarcopterygians for the complete sample of pectoral appendages (Fig. 3A), thus grouping together lobe-fins and limbs and singling out ray-fins. The signal-to-noise ratio (SNR) of this estimated pattern was higher than one (PC1 = 1.26; PC2 = 4.68), which indicated a high effect size of both variables in discriminating groups and adequate power to detect shifts. Comparisons of alternative evolutionary models using Akaike information criteria (AIC) weights showed that a one-regime Ornstein-Uhlenbeck (OU) model best explains the evolution of pectoral appendages (table S1); however, we could not discard that a multi-rate Brownian motion (BM) model could also explain it because AIC values were very close. When comparing fin-bearing taxa with limb-bearing taxa, a phylogenetic MANOVA did not detect a statistically significant difference between them (*F*_{2,21} = 4.79, *P* = 0.238). A similar pattern was found when we analyzed only those taxa for which we had a pelvic correspondence in the sample (fig. S1). The SURFACE analysis for the complete sample of pelvic appendages also estimated a single shift of mean values at the origin of sarcopterygians (Fig. 3B). However, the SNR of the estimated pattern was uneven for each variable (PC1 = 0.69; PC2 = 2.12), which indicated a lack of effect size and power to discriminate between groups for PC1. Comparisons of alternative models using AIC weights showed that an OU model with multiple optimal means (θ) better fitted the evolution of pelvic appendages (table S1). However, a multi-rate BM model yielded similar goodness-of-fit results and could not be discarded. When comparing fin-bearing taxa with limb-bearing taxa, a phylogenetic MANOVA detected a statistically significant difference between them (*F*_{2,16} = 16.85, *P* = 0.0024). A similar pattern was found when we analyzed only those taxa for which we had data for both pelvic and pectoral appendages in the sample (fig. S1B). The congruence of results added support to an estimated shift in topological organization prior to the fins-to-limbs transition at the origin of sarcopterygian appendages. However, it seemed to be possible to discriminate between fins and limbs in pelvic appendages. Overall, limitations in our sample size and methods require cautious interpretations of the estimates of fitted of evolutionary shifts, although the SNR suggested that there was enough power to detect them.

Mean relative subclade topological DTT decreased for pectoral and pelvic appendages alike (Fig. 4). The DTT profile of the complete samples showed a statistically significant deviation from the 95% confidence envelope for pectoral [MDI (morphological disparity index) − 0.41, *P* = 3 × 10^{−4}] and pelvic appendages (MDI = −0.41, *P* = 1 × 10^{−4}), especially after a point in time roughly concurrent with the age of the oldest purported tetrapod tracks (Fig. 4, light green dashed line). The negative MDI values indicated that lineages had nonoverlapping topological variables and that their mean relative subclade disparity was lower than expected. Because pectoral and pelvic samples differed in some taxa, we also analyzed a subsample of both appendages to compare them. The subsampled DTT showed similar patterns between the pectoral and pelvic appendages, with a change in disparity roughly coincident with the age of the oldest tetrapod tracks (fig. S2); the only difference was that the subset showed an almost invariant disparity after the early Middle Devonian time frame.

### Pectoral-pelvic similarity bottleneck

To test the similarity bottleneck hypothesis, we performed phylogenetic generalized least square (PGLS) regressions of pectoral-pelvic disparity against time within the clade Tetrapoda for each of the topological variables. As a proxy for similarity, we used disparity, calculated as the absolute residual of pectoral and pelvic network variables relative to the identity line (lower disparity = higher similarity). The bottleneck hypothesis predicted regression slopes significantly greater than zero for disparity against time. Because the bottleneck marked the point of minimum disparity (maximum similarity), taxa far from the bottleneck would have had higher disparity (lower similarity) than taxa close to the bottleneck. None of the PGLS regressions were statistically significant (Fig. 5), which meant that the topological arrangement of the pectoral-pelvic appendages did not support a similarity bottleneck at the origin of Tetrapoda. We confirmed this result by performing phylogenetic Spearman’s correlation nonparametric tests, which showed no statistically significant correlation (rho’s *P* > 0.05) between each parameter’s pectoral-pelvic disparity and time from the bottleneck.

### Alternative anatomical interpretations

In building anatomical networks of appendicular skeletons, we included only bones and articulations whose presence we were confident about (“minimal networks”). However, 18 of 46 appendages modeled had one or more articulations for which we were uncertain, mostly involving mediolateral contacts between bones in extinct taxa. We accounted for variations due to the presence/absence of these articulations by analyzing network models that also included these uncertain articulations (“extended networks”). Regardless of minor changes in the values of some topological variables, the main results and interpretation for extended networks were the same as for minimal networks (see details in figs. S12 to S15 and table S6). This indicated that our analysis could accommodate informed variations in network models, due to alternative interpretations of the anatomy of appendages, without nullifying the main results and conclusions of this study.

## DISCUSSION

The evolution of tetrapod limbs from fish fins occurred through a series of anatomical changes, including the loss/gain of girdle elements, acquisition of wrist/ankle joints, and the development of digits (*5*). Here, we demonstrated that digits had the greatest impact on the evolution of the topological anatomy of the appendicular skeleton. Topological variables discriminated between digit-bearing taxa and radial-bearing taxa, which occupied distinct regions of the morphospace (Fig. 2). Studies based on different appendicular skeletal traits also discriminated between digit-bearing taxa and radial-bearing taxa (*31*). From a topological viewpoint, digits increased the number of bones (N) and articulations (K) of appendages but in such a way (in contrast to carpal/tarsal topological arrangements of bones) that both the relative density of articulations (D) and the heterogeneity of number of articulations of bones (H) decreased. At the same time, phalanges articulating one-on-one proximodistally increased the overall path length (L) of appendages and the assortativity of bones (A). Last, the formation of new digit modules increased the parcellation (P) of limbs compared to lobe-fins, making them more modular (*23*). The origin of digits was concomitant with the parallel evolution of pectoral and pelvic limbs toward a region of the morphospace (Fig. 2A, top right), where digit-related features predominated (high N, K, A, and P; low D and H). These findings are not a surprise if one acknowledges the presence of digits as a morphological and evolutionary novelty (*27*, *32*–*34*), regardless of shared developmental histories or “deep homologies” between digits and radials (*35*–*39*). Novelties that increase the potential topological morphospace provide an opportunity for greater diversification (*40*). The separation of digit-bearing taxa in the morphospace is in line with the idea that the evolution of limbs represented a truly novel anatomical reorganization (*41*). This novelty has well-recognized functional implications: By partitioning fins into more modular limbs, there should have been more potential for localized functional specializations of bones, joints, muscles, and more, including differential functions of pectoral and pelvic appendages (*9*). Such potential would also arise from the modified degrees of freedom of the limb, transformed from somewhat homogeneous, more flexible fins into fewer, more stiffened distal limb joints (i.e., in the mesopodium and autopodium) (*42*). This is a logical speculation, although our analysis of bone topological networks is unable to account for details of the lepidotrichia, which overall should be less stiff and more numerous than distal limb joints. To our surprise, we could not detect these topological differences as an evolutionary shift through the phylogeny of early tetrapods (Fig. 3). Here, our estimations find only an evolutionary shift at the origin of sarcopterygians that set lobe-fins and limbs apart from ray-fins. This result suggests a conserved, constrained pattern among sarcopterygian appendages, which share some organizational features in the proximal regions, such as the topological arrangement of the girdles, stylopod, and parts of the zeugopod (*15*).

The topological disparity of pectoral and pelvic appendages also decreased during the time span between the origin of sarcopterygians and the origin of amniotes (Fig. 4). Disparity decreased more markedly until approximately the time of the earliest described tetrapod trackways [e.g., the Zachełmie tracks from the Middle Devonian (*43*), similar in timing to other records predating substantial body fossils (*44*, *45*)] and then stabilized. This pattern, although based on a modest sample size as mandated by the existing fossil record, is congruent with an evolutionary stabilization of the disparity of appendicular anatomical organization near the origin of tetrapods. If we consider these earliest tetrapod trackways part of the fins-to-limbs transition, our results would suggest that the transition indirectly decreased the morphological variation, which may have constrained the evolution of different topologies in limbs (see, for example, the morphospace occupation of limbs in Fig. 2C). This conclusion would agree with a dynamic compromise between possibly different functional demands in pectoral and pelvic appendages during the water-to-land transition (*9*) and a shared developmental program constraining the evolvability of limbs (*46*) [constraints that are still strong even in more deeply nested tetrapod lineages like primates (*47*)]. Differing constraints and perhaps compromises have also been proposed to explain a decrease of disparity in the lower jaw of tetrapodomorphs across the water-to-land transition (*48*). However, while functional trade-offs are likely for the feeding versus locomotor systems of stem tetrapods, this decoupling remains to be studied for the craniocervical region, which was likely more constrained in tetrapodomorphs by mechanical interactions between the pectoral appendages, axial column, and skull.

Previous studies have suggested the presence of bottlenecks in pectoral-pelvic anatomical similarity during the evolution of vertebrates (*7*, *21*, *24*–*26*). One similarity bottleneck was proposed at the origin of tetrapods coincident with the fins-to-limbs transition (*7*, *26*). This bottleneck was supported by data from general morphological features (e.g., shape and size similarities and presence/absence of homologous bones) (*3*) and overall configuration and number of muscles (*26*). Our disparity versus time regression tests (Fig. 5) rejected the presence of this bottleneck for the eight topological variables measured for taxa across the fins-to-limbs transition. These new tests override our previous tentative support for this bottleneck from the analysis of the network-anatomy of extant taxa alone (*23*). Our result, based on the network-anatomy or topological organization of the skeleton, differs from previous observations based on different skeletal and muscular features. One possible resolution to this apparent contradiction is that bones and muscles may have had a decoupled evolution during the fins-to-limbs transition, mirroring the idea that bones and muscles can respond differentially in time and magnitude to evolutionary pressures (*49*–*51*). Lingering challenges for this question include the differing nature of data in bottleneck analyses [e.g., phylogenetic characters (*3*) and muscular attachments (*26*)] relative to this study, in addition to other complex traits not yet considered in such analyses, such as lepidotrichia.

Our study of the network-anatomy of appendages during the fins-to-limbs transition revealed an overall parallelism in the evolution of pectoral and pelvic appendages during this time, shaped greatly by the origin of digits. Digits were a morphological novelty that significantly changed the topological features of appendages, discriminating limbs from fins and even from transitional forms. The presence of digits set apart appendages that were less densely and heterogeneously connected but more assortative and modular. Digits evolved in the context of a general decrease in topological disparity among pectoral as well as pelvic appendages, which may have had an impact on the subsequent evolution of tetrapods in terms of function, behavior, and ecology.

## METHODS

### Anatomy of extinct taxa

We examined the skeletal anatomy of the pectoral and pelvic appendages in 18 extinct taxa: *Sauripterus taylory* Hall 1843; *Sterropterygion brandei* Thomson 1972; *Eusthenopteron foordi* Whiteaves 1881; *Cabonnichthys burnsi* Johanson & Ahlberg 1997; *Mandageria fairfaxi* Johanson & Ahlberg 1997; *Gogonasus andrewsae* Long 1985; *Panderichthys rhombolepis* Gross 1941; *Tiktaalik roseae* Daeschler, Shubin & Jenkins, 2006; *Acanthostega gunnari* Jarvik 1952; *Ichthyostega* sp. Säve-Söderbergh 1932; *Tulerpeton curtum* Lebedev 1984; *Balanerpeton woodi* Milner & Sequeira 1994; *Eryops megacephalus* Cope 1877; *Celtedens ibericus* McGowan and Evans 1995; *Seymouria baylorensis* Broili 1904; *Westlothiana lizziae* Smithson and Rolfe 1990; *Pantylus cordatus* Cope 1881; and *Hyloplesion longicostatum* Fritsch 1875. Our resources included museum collections, photographs, and literature descriptions (see details in the Supplementary Materials). These taxa were selected because they have articulated specimens with complete pectoral and/or pelvic appendages for examination and/or described in the literature. We only considered incomplete or disarticulated materials when a full, rigorous reconstruction of the appendage was available in the literature. Last, we decided to exclude from the analysis those appendages for which the complete skeletal anatomy could not be confidently reconstructed due to a large number of missing elements, namely, pectoral appendages of *Acanthostega*, *Ichthyostega*, and *Westlothiana* and pelvic appendages of *Sauripterus*, *Gogonasus*, *Cabonnichthys*, *Mandageria*, and *Tiktaalik*.

### Anatomy of extant taxa

We examined the skeletal anatomy of pectoral and pelvic appendages in nine extant taxa. Six of them were recently described elsewhere (*23*): *Polypterus senegalus* Cuvier 1829, *Latimeria chalumnae* Smith 1939, *Neoceratodus forsteri* Krefft 1870, *Ambystoma mexicanum* Shaw 1789, *Salamandra salamandra* Linnaeus 1758, and *Sphenodon punctatus* Gray 1842. In addition, we built network models for *Iguana iguana* Linnaeus 1758, *Didelphis virginiana* Kerr 1792, and *Mus musculus* Linnaeus 1758 (see details in the Supplementary Materials). We selected these extant taxa because there were available dissection data and because they bracket the fins-to-limbs transition (i.e., rootward and crownward relative to Tetrapoda/Amniote).

### Network modeling

We built unweighted, undirected network models for the appendicular skeleton, where nodes coded for bones and links connecting nodes coded for physical articulation or contacts between two bones. Network models included the girdle and fin/limb skeleton. For the girdles, we considered all skeletal elements present or presumed as present: In pectoral girdles, these may include interclavicle, clavicle, supracleithrum, anocleithrum, cleithrum, and scapulocoracoid; in pelvic girdles, these may include the hip bones fused (pelvis) or divided into two or three parts (ilium, pubis, and ischium). For fin and limb skeletons, we considered all endochondral elements with a sufficient degree of ossification to be directly observed, as well as those elements for which there was enough indirect evidence (for example, an articular surface in another bone). We decided to exclude peripheral dermal elements, such as lepidotrichia and scales, from the fin network models for two main reasons. First, it is often impossible to precisely identify their physical contacts to other elements in fossil taxa; second, their absence in digit-bearing taxa adds noise to the comparison of the skeletal topology between fins and limbs using network analysis.

We coded the articulations among bones following detailed descriptions of each taxon (see the Supplementary Materials). When in doubt, we considered physical contiguity and adjacency as presence of articulation, which allowed us to code for contacts between bones in fossils that did not preserve details of the articular surface due to lack of preservation. Nevertheless, it was sometimes difficult to discern the presence/absence of a given contact between two bones in fossil taxa. We tackled this uncertainty at the modeling and analysis levels. At the modeling level, we tested uncertainty by building two types of networks for each appendage: a minimal network that includes the contacts with high certainty and an extended network that includes also potential but more uncertain contacts (see the Supplementary Materials). This assesses whether different criteria may affect the evolutionary patterns reported. At the analysis level, we tested uncertainty by performing a robustness test of network parameter values under the assumption of random noise or sampling error (see below).

### Anatomical network analysis

We characterized the architecture of fins and limbs using eight topological variables (network parameters): number of nodes (N), number of links (K), density of connections (D), mean clustering coefficient (C), mean path length (L), heterogeneity of connections (H), assortativity of connections (A), and parcellation (P). In short, parameters N and K are counts of the number of bones and physical contacts among bones, respectively. D measures the actual number of connections divided by the maximum number possible (it ranges from 0 to 1); D increases as new bones evolve if they form many new articulations, and decreases otherwise. C measures the average of the ratio of a node’s neighbors that connect among them so that two nodes connected to the same node are also connected (it ranges from 0 to 1); in the appendicular skeleton, triangular motifs can form by adding mediolateral articulations to the most commonly present proximodistal ones. L measures the average number of links required to travel between two nodes (minimum 1); L will increase, for example, by the presence of serial bones articulating one on one proximodistally. H measures the variability in the number of connections of nodes as the ratio between the standard deviation and the mean of the number of connections of all nodes in the network (minimum 0); appendages where each bone has a different number of articulations will have a high H, while if bones have the same number of articulations (i.e., forming a regular pattern) the appendage will have a low H. A quantifies the extent to which nodes with the same number of connections connect to each other (positive if nodes with the same number of connections connect to each other, and negative otherwise); when this happens, A is positive, whereas the inverse tendency means that A is negative. Last, P measures the degree of modularity of the network (it ranges from 0 to 1); appendages with more network modules and with bones evenly distributed among modules will have a high P. See the Supplementary Materials for further mathematical details. We measured topological variables in R (*52*) using functions from the package igraph (*53*).

### Parameter robustness

We tested the robustness of topological variables to potential errors in assessing the presence of bones and articulations by comparing the observed values to a randomly generated sample of 10,000 noisy networks for each anatomical network. We created noisy networks by randomly rewiring the links of the original network with a 0.05 probability, which results in introducing a 5% artificially generated error. Then, we compared the observed values of empirical networks to the sample of noisy networks. In each case, we tested the null hypothesis that observed values are equal to the sample mean. We rejected the null hypothesis with α = 0.05 if the observed value is in the 5% end of the distribution of simulated values (table S2; “TRUE,” cannot reject H_{0}; “FALSE,” reject H_{0} with α = 0.05). We tested a total of 272 values (34 networks × 8 parameters): 268 fell within the confidence intervals and scored TRUE in the test. The exceptions were for *Neoceratodus* pectoral path length, *Neoceratodus* pelvic clustering coefficient and path length, and *Didelphis* pelvic parcellation. Because the anatomy of *Neoceratodus* and *Didelphis* derived from our own dissections, these few cases of rejection of the null hypothesis for these parameters can be attributed to the difficulty for a random noise process to produce realistic dissection errors (for example, by coding the femur as not articulating with the pelvis).

### Phylogenetic relationships

We assembled a phylogenetic tree for our study taxa using the supertree recently published by Ruta and co-workers (*54*). Following Ruta *et al*.’s work, we calibrated the tree branches using the “equal” method defined by Lloyd and colleagues (*55*, *56*) to adjust zero-length and internal branches, as implemented in the package paleotree (*57*) for R. This method is the most conservative to time-calibrate a phylogeny: It makes few assumptions about divergence times and relies exclusively upon first appearance dates (*54*). Temporal ranges of taxa [first and last appearance date in million years (Ma)] follow those of Ruta *et al*. (*54*), the Paleobiological Database (https://paleobiodb.org/), and TimeTree (www.timetree.org) (tables S4 and S5). We constrained tree calibration by assigning minimum dates for known internal nodes (clades or splits) based on molecular inferences and fossil dates from the literature (*54*). The exact dates for first and last appearance of taxa and for internal nodes are available within source code included. Note that when required by the analysis, our phylogenetic tree was pruned to only include those taxa of interest.

### Analysis of topological variation

We ran a PCA of topological variables by a singular value decomposition of the centered and scaled measures using the function prcomp in the R built-in package stats (*52*). We used PCA components to test whether the anatomical organization of the appendicular skeleton differed (i) between fins (without digits) and limbs (with digits), (ii) between pectoral and pelvic appendages, and (iii) between extinct and extant taxa. We performed a PERMANOVA over 10,000 permutations using the function adonis in the R package vegan (*58*). PERMANOVA used a permutation test with pseudo-*F* ratios on the Euclidean distances of the matrix of PCA components to test the null hypothesis that the centroids and dispersion were equivalent for each group comparison. Rejection of the null hypothesis meant that the network topology differed between the groups compared.

### Evolutionary modeling

We estimated the occurrence of evolutionary shifts in the topological organization of appendages in our phylogenetic tree using a SURFACE (*59*) analysis of the first two PC components for pectoral and pelvic appendages independently. SURFACE estimates change of evolutionary regimes—i.e., changes in the strength (α) and rate (σ) of evolution and in the optimal mean (θ) in the branches and/or clades of a tree—from multivariate data and a non-ultrametric tree. SURFACE uses an OU stabilizing selection model of evolution, which allows changes in the rate of evolution and optimal means of variables. If present, this method identifies homoplasy: two clades with the same regime. Given the small sample of appendages in our comparisons, we deemed it necessary to calculate the power of the SURFACE analysis as an indicator of reliability in the accuracy of estimated patterns. Because in OU models power is dependent by strength, rate, and optimal mean combined, effect-size measures offer a better prediction of power than sample size (*60*). We calculated the power of the estimated regimes using the SNR, which is defined as *T* is the total depth of the phylogeny. High power can be inferred when SNR > 1. To further validate the estimated regimes, the output of SURFACE was then fitted to alternative evolutionary models: a one-rate BM; a multi-rate BM; an OU with fixed strength, rate, and mean; an OU with fixed strength and rate, and multi-mean; and an OU with fixed strength, and multi-rate and multi-mean. Fitted models were compared using AIC. Last, we statistically tested the resulting evolutionary models of evolution with a phylogenetic MANOVA to confirm that the clades identified had a different regime corresponding to different groups with different means. The combination of estimation, fitting, and testing allowed us to build confidence that the evolutionary patterns found were reliable if they converged on the same result. Evolutionary modeling was carried out in R using functions from the packages surface (*59*), mvMORPH (*61*), and geiger (*62*) for the estimation, fitting, and testing, respectively.

### Disparity through time

To examine how topological disparity changed over time, we calculated the mean relative subclade DTT on pectoral and pelvic appendages separately. Following a previous study on mammalian neck anatomy (*63*), we first obtained the covariation of topological variables by performing independent PCAs for pectoral and pelvic networks. Next, we calculated the mean subclade disparity on the PC scores using the function dtt in the R package geiger (*62*). The higher the disparity, the higher the variance within subclades (i.e., lower conservation) and the lower the variance between subclades (*64*, *65*). We tested the statistical significance of the observed disparity with a randomization inference test with 10,000 simulations under a BM evolution on our phylogeny. Function dtt also calculated the morphological disparity index, which quantified the overall difference in relative disparity of a clade compared to that expected under the null BM model. To offer a better pectoral-to-pelvic comparison, DTT analyses were also performed on the subsample of taxa for which we have both pectoral and pelvic appendages.

### Pectoral-pelvic similarity bottleneck

We tested the hypothesis of the existence of a pectoral-pelvic similarity bottleneck at the origin of Tetrapoda for each topological variable independently. For practical purposes, we used pectoral-pelvic disparity (lower disparity means greater similarity). Pectoral-pelvic disparity was calculated as the absolute residuals of pectoral and pelvic values on the identity line (or 1:1 line, a line with intercept = 0 and slope = 1) so that identical pelvic and pectoral appendages—maximal similarity—had a value of zero disparity. According to previous formulations of the bottleneck hypothesis for tetrapods, taxa before and after the split of tetrapods would have a greater pectoral-pelvic disparity (lower similarity) than taxa closer to the origin of tetrapods (*7*, *23*, *26*). Thus, the farther we go in time from this event, the greater the expected pectoral-pelvic disparity should be. To test this prediction, we performed a PGLS of the absolute pectoral-pelvic residuals on the 1:1 line against the time from the Tetrapoda branch of taxa having both appendages in the sample (*t*_{0} = 370.4 Ma). We used PGLS to test against the null hypothesis of a slope = 0, meaning no difference in pectoral-pelvic DTT. For our prediction to hold, PGLS needed to show a statistically significant positive regression slope. PGLS was computed in R using a standard generalized least square with an a priori correlation structure derived from the phylogenetic tree using the function corPagel of the package ape (*66*). Because some parameters may not follow some of the assumptions of PGLS, we also performed a phylogenetic Spearman’s correlation test (a nonparametric test with fewer assumptions) to test whether there was a relationship between pectoral-pelvic disparity and time from the bottleneck; a lack of correlation means no difference in pectoral-pelvic DTT.

## SUPPLEMENTARY MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/5/5/eaau7459/DC1

Table S1. Model fitting of SURFACE estimations (wAIC).

Table S2. Result of the robustness test.

Table S3. Summary of the community structure search in anatomical networks.

Table S4. First and last appearance dates used to calibrate the phylogeny.

Table S5. Internal nodes used to constrain the calibrated phylogeny.

Table S6. Values of network variables measured on each extended network model.

Fig. S1. Estimated evolutionary shifts using SURFACE.

Fig. S2. Topological DTT for pectoral and pelvic appendages.

Fig. S3. Visual representation of network parameters [from (*23*)].

Fig. S4. Diagnostic plots of pGLS model assumptions for the number of nodes.

Fig. S5. Diagnostic plots of pGLS model assumptions for the number of links.

Fig. S6. Diagnostic plots of pGLS model assumptions for the density.

Fig. S7. Diagnostic plots of pGLS model assumptions for the clustering coefficient.

Fig. S8. Diagnostic plots of pGLS model assumptions for the shortest path length.

Fig. S9. Diagnostic plots of pGLS model assumptions for the heterogeneity.

Fig. S10. Diagnostic plots of pGLS model assumptions for the assortativity.

Fig. S11. Diagnostic plots of pGLS model assumptions for the parcellation.

Fig. S12. Biplots of the first two PCA components of topological variables for pectoral (circles) and pelvic (triangles) appendages combined.

Fig. S13. Estimated evolutionary shifts using SURFACE.

Fig. S14. Topological DTT for pectoral and pelvic appendages.

Fig. S15. Testing the pectoral-pelvic similarity bottleneck hypothesis for Tetrapoda.

Supplementary Methods

Supplementary Results

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is **not** for commercial advantage and provided the original work is properly cited.

## REFERENCES AND NOTES

**Acknowledgments:**We thank J. Clack, T. Smithson, M. Coates, and R. Benson for their help in interpreting skeletal anatomy in fossil specimens; D. Brinkman and A. Heimer (Yale Peabody Museum of Natural History) for providing correct number and high-resolution images of

*Sterropterygion brandei*; J.R.H.’s team at the Royal Veterinary College for their input on early drafts; J. Smaers for teaching and insight on phylogenetic analyses; and the Transmitting Science team for making it possible. We also thank J. Cundiff (Museum of Comparative Zoology at Harvard) and E. Bernard (Natural History Museum of London) for their assistance in accessing specimens. S.E.P. also thanks the support by Harvard University funds.

**Funding:**This project was funded by the European Union’s Horizon 2020 research and innovation program under Marie Skłodowska-Curie grant agreement no. 654155.

**Author contributions:**B.E.-A., S.E.P., and J.R.H. built and revised the adjacency matrices for extinct taxa. B.E.-A., J.L.M., P.J., J.R.H., and R.D. built and revised the adjacency matrices for extant taxa. B.E.-A. designed the study, analyzed the networks, and wrote the manuscript. All authors discussed the results and revised the manuscript.

**Competing interests:**The authors declare that they have no competing interests.

**Data and materials availability:**All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from B.E.-A. at borja.esteve{at}upf.edu. The data and code that support the findings of this study are available from Figshare at https://figshare.com/s/4c3a1dd62d1f55728a0e.

- Copyright © 2019 The Authors, some rights reserved; exclusive licensee American Association for the Advancement of Science. No claim to original U.S. Government Works. Distributed under a Creative Commons Attribution NonCommercial License 4.0 (CC BY-NC).