Data from: Lungless tadpoles breathe fresh air into hypotheses for tetrapod lung loss and trait regain

Supplementary Materials - Phillips et al., 2025 Table of contents Supplementary methods - page number 1. Data acquisition methods - 2 Dissections - 2 Contrast enhanced microCT - 2 Air breathing as evidence of lungs - 2 Literature search - 3 Assignment of ecomorphological categories - 3 2. Phylogeny methods - 5 Figure S1 - 6 3. Phylogenetic comparative methods and implementation a. BayesTraits extended methods - 5 b. Full list of BayesTraits models tested - 8 c. Bayes Factor estimation and interpretation - 11 d. CorHMM methods - 11 e. Stochastic character mapping and ancestral state estimation - 12 f. Phylogenetic generalized linear models - 12 Supplementary Results - 13 1. Supplementary Table 1 (Table S1) Legend - 13 2. Full BayesTraits model comparisons - 14 Table S2: summary of all BayesTraits models - 14 3. Visualizations of different evolutionary histories for different BayesTraits models - 15 to 19 Figs. S2, S3, S4, S5, and S6 4. Full CorHMM outputs - 20 to 24 Table S3: summary of CorHMM models 5. Phylogenetic generalized linear model results - 25 to 26 Table S4: summary of phyloglm models References - 27 on 2 Supplementary Methods 1 - Data Acquisition Methods Dissections Specimens of 315 species were dissected under a dissecting scope while pinned to a silicone mat at the bottom of a plastic container that was filled with 70% ethanol, 10% formalin, or water, depending on curatorial preferences. To access the lungs, we cut into the ventral skin and pulled back the ventral musculature to reveal the spiral intestine. If the lungs are large, they can sometimes be visualized immediately with no further damage to the specimen (e.g., Fig. 1A). In other cases, small lungs and lung buds could only be visualized by removing the intestine (e.g., Fig. 1F) and the smallest lung buds required careful examination of the most anterior section of the gut tube (e.g., Fig. 1G). In rare cases (N = 3), lung buds could not be reliably identified at all. Some specimens examined by PHD were previously stained with alcian blue solution for musculo-skeletal studies, with no damage to the lungs. We scored tadpole lungs as “absent” if lungs could not be found, if the lung buds were solid and obviously non-functional, or if they were less than 5% the length of the abdominal cavity (Figure 1). Contrast enhanced micro-computed tomography (CEμCT) scanning We also obtained lung presence/absence data from CEμCT scanning for 23 species. Specimens of these species were stained with either phosphotungstic acid in 70% Ethanol or an aqueous Lugol’s iodine stain and then reconstructed into voxels for a variety of other projects. If lung tissue could be clearly identified in the scan, then a judgment of present (i.e., functional) vs. absent (i.e., nonfunctional) was made. However, if lung tissue could not be identified or the lungs appeared ambiguous, we performed a dissection to confirm lung status. Those species are counted as “dissection” in Supplemental Table 1, which shows the evidence used for each taxon included in the analyses. Air breathing as evidence of lungs We also kept track of air breathing observations as a way to test for lung presence, which we observed in 27 species. However, in every case but one (Cruziohyla silviae tadpoles observed 3 breathing near Gamboa, Panama by JRP), lung status was confirmed either through dissections, CEμCT scan, or literature reports. Literature search etc. In addition to our primary dataset of dissections and CEμCT Scanning (338 species), lung presence/absence data was obtained from previously published studies, books, and photographs as well as correspondences with experts in specific taxa. The presence/absence of lungs was outright mentioned in some cases, while in others the presence of lungs was inferred by the reported occurrence of air-breathing behavior. Other studies made no remark on lungs, but presented photographs that showed the lungs either in dissection or in life, as many tadpoles are highly transparent, and the lungs can be visualized from the outside. We also used 25 photographs of tadpoles that clearly showed inflated lungs through the body wall as evidence of lung presence. We were able to confirm reported lung status in the literature with dissections of 23 taxa and an additional species via CEμCT.. Supplemental table 1 includes the sources, evidence type used, and sample sizes for all taxa examined. We sampled as widely as possible, prioritizing the inclusion of as many families and genera as possible. In total, we were able to sample 252 of 376 genera with free-living larvae and 44 of 47 possible families. The frog families Brachycephalidae, Brevicipitidae, Caligophrynidae, Ceratobatrachidae, Ceuthomantidae, Craugastoridae, Eleutherodactylidae, Neblinaphrynidae, and Strabomantidae are so-called “direct developers” and lack free-living larvae, meaning they could not be included in this project. Supplemental Fig. 1 shows the distribution of sampled taxa within a sample of 5,242 anurans with non-viable taxa (genera with no known free-living stage) shown in orange. Assignment of ecomorphological categories We assigned tadpoles to a set of overlapping categories based on the general ecomorphological guilds of (Altig and Johnston 1986). We first split tadpoles into aquatic and terrestrial. Most free-living 4 tadpoles are aquatic, but some are either fully terrestrial or “semi-terrestrial”. Fully terrestrial tadpoles are often endotrophic (not free-feeding), being parent-associated and/or living in nests (e.g., Sooglossus and Allobates nidicola). Semi-terrestrial tadpoles are considered a lotic guild and consist of tadpoles that live primarily on wet rock faces alongside streams (56). We characterized all of these tadpoles as “terrestrial”, and all other tadpole forms as “aquatic”. We next split tadpoles into lentic and lotic categories. We scored breeding location as lentic when the description referred to a pond, pool, puddle or similar water body, and as lotic when the description referred to a stream, torrent, rapid, or rocky surface associated with such environments. We scored species that were reported to breed in both habitat types (N=5) as lentic, because their tadpoles are generally more similar to lentic tadpoles (Altig and Johnston 1989). “Lotic, semi-terrestrial” tadpoles were counted as lotic for four-state analyses and as terrestrial in eight-state analyses. Within lotic tadpoles, we further categorized tadpoles as either “specialized” or “other lotic”. To clarify, specialized is a shorthand for two specific types of specialization, there are many highly specialized, lotic tadpoles that would be included as “other” in our analysis. We included tadpoles as specialized in cases where they either had some form of suction apparatus (suctorial) or were specialized for burrowing (fossorial). Suctorial tadpoles use suction to attach to rocks in fast-flowing torrents, including the previously defined guilds: lotic, suctorial (mouth-suckers) and the bizarre lotic, gastromyzophorous guild, which use a belly-sucker (78, 79). Fossorial tadpoles are long and worm-like, living under and within dense substrates like sand, gravel, and leaf mats (80, 81). We included the ecomorphological guilds “lotic, fossorial” and “lotic, psammonic” (sand-feeding) under the single umbrella of “fossorial” as these groups all live buried under solid substrata. We considered tadpoles specialists if they had been placed in any of the above guilds by (23) or others. In rare cases where a tadpole had not been placed in any guild, we assigned it ourselves (suppl. Table 1), based on ecological and morphological evidence. All other lotic guilds were scored as “other”. 5 2 - Phylogeny Methods We prioritized including as many sampled taxa in our phylogenetic analyses as possible (Suppl. Fig. 1). In cases where we sampled taxa for lung data and the taxon was not present in (Portik et al. 2023), we attempted to place those taxa in the tree by substituting for a present taxon we did not sample. We accepted a substitution of a sampled taxon for one present in the tree under three scenarios. First, if the sampled taxon that is missing in the tree is the only member of its genus sampled, we assumed generic monophyly and substituted the sampled taxon for a random member of the genus present in the tree. Second, if the sampled taxon has been analyzed phylogenetically by another study such that it could be aligned with a taxon present in the tree lacking lung data, we substituted the sampled species for the species in the tree. Third, if we sampled two members of a Genus, and only one was present in the tree, we randomly chose a third, present species in the genus to substitute for the sampled species. This third justification assumes generic monophyly, and produces trees with reasonable bifurcating relationships, but potentially alters “true” branch lengths. We believe the value of adding additional species to the analysis outweighs the potential negative effects of incorrect branch lengths, especially because both ecology and lung status are often uniform within genera, and there is already uncertainty in branch length estimates. Before making any substitutions, we first checked whether each genus was monophyletic in the published tree, and if not, we did not assume monophyly. See supplementary Table 1 for a full list of which taxa were substituted for each tree, with citations. 3 - Phylogenetic Comparative Methods implementation A - BayesTraits Extended Methods We used the program BayesTraits (BT) (V.5.0 windows version) (Pagel 1994; Pagel and Meade 2006, n.d.) to estimate and test Mk models. We used a custom-built R framework to analyze and visualize BayesTraits outputs. Code used to prepare BayesTraits datasets and analyze the outputs are available at https://github.com/phillipsjar/Phillips_etal_2025_tadpole_lung_loss/ and BayesTraits command files and actual outputs are available as supplemental files here: https://doi.org/10.5061/dryad.3tx95x6sk. Q matrices are included below. 6 For all analyses, we used the Bayesian mode with a rjMCMC sampler with a uniform prior from 0 to 100 for all potentially free parameters. Following the BT manual, we scaled trees down because the branch lengths in (Portik et al. 2023) correspond to millions of years rather than base pair substitutions, which is the assumed scale for BT. We scaled trees by running sets of full models without a scaling factor, and then scaling the tree by three factors of ten (.1, .01, .001). We chose the scale that resulted in the highest parameters being estimated above 10 and not leading to a truncation at the upper end of our priors (100). Scaling proportionally increases mean parameter estimates, preventing BT from failing to distinguish low rates with zero. The BT manual (p. 10) recommends scaling branch lengths to a mean of .1 to avoid this problem, but we preferred to scale intentionally, such that rates can easily be 7 back-transformed later on. We scaled the 4-state and 8-state models independently, but any models compared with 2logBayesFactors were scaled to the same amount. After testing different scales, we chose to scale the four-state models by a factor of 0.001 and the eight-state models by a factor of 0.01. For each run of the 4-state models, we used 2,500,000 iterations per run with a burnin of 500,000 and three runs per model. For each run of the eight-state models, we used 110,000,000 iterations and a burnin of 20,000,000 iterations. Eight-state models were much more complex, and so slower to run and converge, so we greatly increased burnin and iterations. We assessed within-run convergence by examining trace plots, and across-run convergence by comparing estimated marginal likelihoods and visually comparing violin plots of posterior data. We selected models for visualization that converged both within and across runs, presenting the best individual run (as determined by the highest estimated log marginal likelihood) in Fig. 4. We calculated highest density intervals as credible intervals for all estimated parameters using the R package HDInterval, v.0.2.4 (Meredith and Kruschke 2020). The more complicated, eight-state models of character evolution that we used required some fossilizations to reliably produce reasonable estimates of character evolution. We found that when unrestricted, the eight-state models often produced solutions that reconstructed the ancestors of most clades as being terrestrial. Because both specialized lotic and terrestrial tadpoles are associated with specific larval apomorphies not found across large parts of the anuran tree, we ruled this solution to not be biologically feasible, and so restricted several key nodes of the tree to be either lentic or generalized lotic, with equal probabilities of being lunged or lungless in either state. We did this for the most recent common ancestor (MRCA) of all frogs, the MRCA of Neobatrachia, the MRCA of the Dendrobatidae family, the MRCA of the Leptodactylidae family, and the MRCA of the Cyclorhamphidae and Alsodidae. The only effect of these fossilizations was to polarize transitions from aquatic tadpoles to terrestrial tadpoles, an assumption strongly supported by biological inference. 8 B - Full list of models tested and Q matrices We compared a set of potentially reasonable evolutionary models of lung evolution in BT for both the four-state and eight-state analyses. Among the four-state analyses, we began with a full dependent model, which included eight potentially unique parameters (lung loss and lung gain in both lentic and lotic habitats, and transitions between lentic and lotic habitats in both the lunged and lungless conditions). We also tested a series of nested, more restricted versions this full model, including a model in which lentic lung loss was restricted to zero, lung regains in lentic and lotic habitats were set equal, lung regains were set to zero, and both lung regains and lentic lung loss were set to zero. We also tested an “independent” model in which lung evolution in any direction was set equal for both habitat types and habitat evolution in either direction was set equal regardless of lung status, leading to only four potentially unique parameters. Among the eight-state models, we again tested a “full” model, which included up to 28 unique parameters. This model included four rates of lung loss (for each habitat classification: lentic, lotic non-specialized, lotic specialized, and terrestrial), an additional four rates of lung regain in each habitat, and then a series of potential transitions between habitat types for both the lunged and lungless statuses. We restricted habitat transitions between lentic and specialized lotic, believing a more realistic scenario requires an evolutionary transition into lotic non-specialized. However, we did allow transitions from all lentic and lotic habitats directly to terrestrial. We also tested a set of more restricted eight-state models, including a model in which lotic, non-specialized loss was set to zero, a model for which lentic lung loss was set to zero, one model with all lung regains were set to zero, another with lung regains all set to a single parameter value, another with lentic lung regains set to zero, and a final model with all regains set to zero and non-specialized lotic lung loss set to zero. 9 Q matrices Four state: Independent Model 1 - lungless lotic 2 - lunged lotic 3 - lungless lentic 4 - lunged lentic 1 - 𝛼2 (gain lungs) 𝛼1 (lotic to lentic) - 2 𝛽2 (lose lungs) - - 𝛼1 (lotic to lentic) 3 𝛽1 (lentic to lotic) - - 𝛼2 (gain lungs) 4 - 𝛽1 (lentic to lotic) 𝛽2 (lose lungs) - Dependent Model 1 - lungless lotic 2 - lunged lotic 3 - lungless lentic 4 - lunged lentic 1 - q12: gain lungs (lotic) q13: lotic to lentic (lungs) - 2 q21: lose lungs in lotic - - q24: lotic to lentic (no lungs) 3 q31: lentic to lotic (no lungs) - - q34: gain lungs (lentic) 4 - q42: lentic to lotic (lungs) q43: lose lungs (lentic) - Four-state model restrictions Dependent: no restrictions Regains set equal: q12 = q34 Regains prevented: q12 = q34 = 0 No lentic loss: q43 = 0 No regains and no lentic loss q12 = q34 = q43 = 0 Independent q12 = q34 ⇒ 𝛼1; q21 = q43 ⇒ 𝛽2; q13 = q24 ⇒ 𝛼1; q31 = q42 = 𝛽1 16 Supplemental Figure 3: 100 blended stochastic character maps of the full four-state dependent evolutionary model with ancestral state estimates, where rates of lung and habitat evolution vary depending on the state of the other trait. Note circled nodes, which differ between independent and dependent models. Q matrix from median values of parameter estimates in the best individual run of the full dependent model (run 2 in Supp. Table 2) and tree is a prune maximum likelihood phylogeny from (Portik et al. 2023). 17 Supplemental Figure 4: 100 blended stochastic character maps of the best four-state dependent evolutionary model with ancestral state estimates, where rates of lung and habitat evolution vary depending on the state of the other trait and rates of regain and lung loss in lentic habitats are set to zero. Note circled nodes, which differ between independent and dependent models. Q matrix from median values of parameter estimates in the best individual run of the best dependent model (run 1 in Supp. Table 2) and tree is a prune maximum likelihood phylogeny from (Portik et al. 2023). 18 Supplemental Figure 5: 100 blended stochastic character maps of the full eight-state evolutionary model (with ancestral state estimates), which includes terrestrial tadpoles as well as distinguishes between lotic tadpoles specialized for suctoriality and fossoriality from “other” lotic tadpoles. Q matrix from median values of parameter estimates in the best run of the full eight state model (run 2 in Supp. Table 2) and tree is a prune maximum likelihood phylogeny from (Portik et al. 2023). 19 Supplemental Figure 6: 100 blended stochastic character maps of the best eight-state evolutionary model with ancestral state estimates, where lung regains are impossible and lung loss cannot occur in lotic tadpoles not specialized for suctoriality or fossoriality. Q matrix from median values of parameter estimates in the best run of the full eight state model (run 1 in Supp. Table 2) and tree is a prune maximum likelihood phylogeny from (Portik et al. 2023). 20 4 - Full CorHMM Outputs from R Rates are unscaled, so based on the scale of the tree’s branch lengths, which are in millions of years. Therefore, rates refer to expected number of changes per million years. Example Q matrices with explanatory rates: Non-Hidden Markov Model (only 1 set of rates across the entire tree): lungless lentic lunged lentic lungless lotic lunged lotic (1,R1) (2,R1) (3,R1) (4,R1) 1,R1 - gain lungs (lentic) lentic to lotic (no lungs) - 2,R1 lose lungs (lentic) - - lentic to lotic (lungs) 3,R1 lotic to lentic (no lungs) - - gain lungs (lotic) 4,R1 - lotic to lentic (lungs) lose lungs (lotic) - Hidden Markov Model (2 possible rate matrices across the tree, each set of transition rates colored): lungless lentic lunged lentic lungless lotic lunged lotic lungless lentic lunged lentic lungless lotic lunged lotic (1,R1) (2,R1) (3,R1) (4,R1) (1,R2) (2,R2) (3,R2) (4,R2) 1,R1 - gain lungs (lentic) lentic to lotic (no lungs) - Rate 1 to Rate 2 - - - 2,R1 lose lungs (lentic) - - lentic to lotic (lungs) - Rate 1 to Rate 2 - - 3,R1 lotic to lentic (no lungs) - - gain lungs (lotic) - - Rate 1 to Rate 2 - 4,R1 - lotic to lentic (lungs) lose lungs (lotic) - - - - Rate 1 to Rate 2 1,R2 Rate 2 to Rate 1 - - - - gain lungs (lentic) lentic to lotic (no lungs) - 2,R2 - Rate 2 to Rate 1 - - lose lungs (lentic) - - lentic to lotic (lungs) 3,R2 - - Rate 2 to Rate 1 - lotic to lentic (no lungs) - - gain lungs (lotic) 4,R2 - - - Rate 2 to Rate 1 - lotic to lentic (lungs) lose lungs (lotic) - 21 Individual outputs for each tested model, with estimates rounded to four decimal places. simplified_independent_model_fit Fit -lnL AIC AICc Rate.cat ntax -367.8601 739.7201 739.7444 1 497 Rates lungless lentic lunged lentic lungless lotic lunged lotic (1,R1) (2,R1) (3,R1) (4,R1) 1,R1 - 0.0066 0.0066 - 2,R1 0.0024 - - 0.0066 3,R1 0.0024 - - 0.0066 4,R1 - 0.0024 0.0024 - simplified_hidden_Markov_independent_model_fit Fit -lnL AIC AICc Rate.cat ntax -352.5233 717.0465 717.218 2 497 Rate s lungless lentic lunged lentic lungless lotic lunged lotic lungless lentic lunged lentic lungless lotic lunged lotic (1,R1) (2,R1) (3,R1) (4,R1) (1,R2) (2,R2) (3,R2) (4,R2) 1,R1 - 0.0431 0.0431 - 0.0357 - - - 2,R1 0.0109 - - 0.0431 - 0.0357 - - 3,R1 0.0109 - - 0.0431 - - 0.0357 - 4,R1 - 0.0109 0.0109 - - - - 0.0357 1,R2 0.0093 - - - - 0.0014 0.0014 - 2,R2 - 0.0093 - - 0.0002 - - 0.0014 3,R2 - - 0.0093 - 0.0002 - - 0.0014 4,R2 - - - 0.0093 - 0.0002 0.0002 - simplified_correlated_model_fit Fit -lnL AIC AICc Rate.cat ntax -338.0255 684.0511 684.1324 1 497 Rates lungless lentic lunged lentic lungless lotic lunged lotic (1,R1) (2,R1) (3,R1) (4,R1) 1,R1 - 0.0001 0.0112 - 2,R1 0.0001 - - 0.0077 3,R1 0.0112 - - 0.0069 4,R1 - 0.0077 0.0069 - 22 simplified_hidden_Markov_correlated_model_fit - best model (lotic lung loss bolded, lentic lung loss in italics) Fit -lnL AIC AICc Rate.cat ntax -323.3033 666.6067 667.0594 2 497 Rate s lungless lentic lunged lentic lungless lotic lunged lotic lungless lentic lunged lentic lungless lotic lunged lotic (1,R1) (2,R1) (3,R1) (4,R1) (1,R2) (2,R2) (3,R2) (4,R2) 1,R1 - 1.00E-09 0.0278 - 0.0376 - - - 2,R1 1.00E-09 - - 13.7856 - 0.0376 - - 3,R1 0.0278 - - 0.0041 - - 0.0376 - 4,R1 - 13.7856 0.0041 - - - - 0.0376 1,R2 0.0058 - - - - 1.00E-09 0.0032 - 2,R2 - 0.0058 - - 1.00E-09 - - 1.00E-09 3,R2 - - 0.0058 - 0.0032 - - 0.0087 4,R2 - - - 0.0058 - 1.00E-09 0.0087 - independent_model_fit Fit -lnL AIC AICc Rate.cat ntax -353.0246 714.0492 714.1305 1 497 Rates lungless lentic lunged lentic lungless lotic lunged lotic (1,R1) (2,R1) (3,R1) (4,R1) 1,R1 - 0.0007 0.0079 - 2,R1 0.0017 - - 0.0079 3,R1 0.0118 - - 0.0007 4,R1 - 0.0118 0.0017 - 23 hidden_Markov_independent_model_fit Fit -lnL AIC AICc Rate.cat ntax -326.5059 673.0119 673.4645 2 497 Rate s lungless lentic lunged lentic lungless lotic lunged lotic lungless lentic lunged lentic lungless lotic lunged lotic (1,R1) (2,R1) (3,R1) (4,R1) (1,R2) (2,R2) (3,R2) (4,R2) 1,R1 - 0.0012 0.1349 - 0.0234 - - - 2,R1 0.008 - - 0.1349 - 0.0234 - - 3,R1 1.00E-09 - - 0.0012 - - 0.0234 - 4,R1 - 1.00E-09 0.008 - - - - 0.0234 1,R2 0.006 - - - - 0.0005 1.00E-09 - 2,R2 - 0.006 - - 1.00E-09 - - 1.00E-09 3,R2 - - 0.006 - 0.0526 - - 0.0005 4,R2 - - - 0.006 - 0.0526 1.00E-09 - correlated_model_fit Fit -lnL AIC AICc Rate.cat ntax -331.4721 678.9441 679.2392 1 497 Rates lungless lentic lunged lentic lungless lotic lunged lotic (1,R1) (2,R1) (3,R1) (4,R1) 1,R1 - 1.00E-09 0.0127 - 2,R1 0.0004 - - 0.0068 3,R1 1.00E-09 - - 0.0022 4,R1 - 0.0158 0.0056 - 24 hidden_Markov_correlated_model_fit Fit -lnL AIC AICc Rate.cat ntax -322.9828 681.9656 683.3966 2 497 Rate s lungless lentic lunged lentic lungless lotic lunged lotic lungless lentic lunged lentic lungless lotic lunged lotic (1,R1) (2,R1) (3,R1) (4,R1) (1,R2) (2,R2) (3,R2) (4,R2) 1,R1 - 1.00E-09 0.0141 - 0.0147 - - - 2,R1 1.00E-09 - - 0.0154 - 0.0147 - - 3,R1 1.00E-09 - - 0.0126 - - 0.0147 - 4,R1 - 0.0455 0.0034 - - - - 0.0147 1,R2 0.0046 - - - - 1.00E-09 0.0151 - 2,R2 - 0.0046 - - 1.00E-09 - - 0.0009 3,R2 - - 0.0046 - 0.0155 - - 1.00E-09 4,R2 - - - 0.0046 - 1.00E-09 0.0096 - Supplementary Table 3: summary of all CorHMM models CorHMM Model AIC deltaAIC number of rates Independent 714.05 47.44 4 indep_HMM 673.01 6.40 10 Dependent 678.94 12.33 8 dep_HMM 681.97 15.36 18 indep_simplified 739.72 73.11 2 indep_HMM_simplified 717.05 50.44 6 dep_simplified 684.05 17.44 4 dep_HMM_simplified 666.61 0 10 25 5 - Phylogenetic generalized linear model results Full outputs from R, with coefficients combined into a single table below Model 1: correlation of breeding habitat (lentic/lotic) with lung presence and absence phyloglm(formula = lung ~ ecology, data = data_aqua, phy = full_tree, method = c("logistic_MPLE"), start.beta = NULL, start.alpha = NULL, boot = 1000, full.matrix = TRUE) AIC logLik Pen.logLik 204.75 -99.37 -97.86 Method: logistic_MPLE Mean tip height: 179.3411 Parameter estimate(s): alpha: 0.001150445 bootstrap mean: 0.0008360564 (on log scale, then back transformed) so possible downward bias. bootstrap 95% CI: (0.0001061738,0.006334944) Model 2: correlation of ecomorphological guild (specialized lotic guilds/other types) with lung presence and absence phyloglm(formula = lung ~ Spec_lotic, data = data_aqua, phy = full_tree, method = c("logistic_MPLE"), start.beta = NULL, start.alpha = NULL, boot = 1000, full.matrix = TRUE) AIC logLik Pen.logLik 145.54 -69.77 -69.08 Method: logistic_MPLE Mean tip height: 179.3411 Parameter estimate(s): alpha: 0.0002628877 bootstrap mean: 0.0002568441 (on log scale, then back transformed) so possible downward bias. bootstrap 95% CI: (0.000102581,0.002927915) Model 3: full model testing the correlation of breeding habitat, ecomorphological guild, and terrestriality with lung presence and absence phyloglm(formula = lung ~ ecology + Spec_lotic + terrestrial, data = data_full, phy = full_tree, method = c("logistic_MPLE"), start.beta = NULL, start.alpha = NULL, boot = 1000, full.matrix = TRUE) AIC logLik Pen.logLik 160.46 -75.23 -71.99 Method: logistic_MPLE