Rapid phenotypic evolution following shifts in life cycle complexity

1. Introduction

Patterns of selection can vary across ontogeny and impart strong conflicts on the evolution of traits [1–4]. This is epitomized in organisms with complex life cycles that may require an individual to perform in starkly different environments at each life stage [5,6]. Metamorphosis between life stages can be an effective retooling mechanism prior to an ecological transition and can ‘adaptively decouple’ the compromises of different selective regimes [6–8]. However, in some organisms, many tissues and structures may remain relatively unchanged even when the rest of the body undergoes an extreme transformation [9–12]. Unaltered traits that are expressed across ontogeny may be limited to a generalized phenotype that is effective under multiple selective regimes. For such traits the loss of ontogenetic complexity is predicted to alleviate constraints and permit trait evolution [6,12], but phenotypic responses at shallow time scales have not been evaluated in a phylogenetic context.

Clades that exhibit variation in life cycle stages are ideal systems to test how ontogeny can influence trait evolution [12]. Most amphibians exhibit a biphasic (metamorphic) complex life cycle with an aquatic larval stage followed by metamorphosis into a more terrestrial adult form [13]. Living in aquatic versus terrestrial environments should impose strongly contrasting demands on the evolution of anatomy and physiology. However, many lineages of salamanders forego metamorphosis and maintain their aquatic larval morphology and ecology throughout ontogeny, a life cycle simplification referred to as larval form paedomorphosis [14,15]. The repeated independent loss of the terrestrial stage through paedomorphosis provides a compelling system to test the impacts of life cycle complexity on trait evolution.

Macroevolutionary analyses show that life cycle simplifications from a biphasic to paedomorphic life cycle have probably had a dramatic effect on salamander body form and vertebral column diversification [12], but have presumably done so under a wide range of ecological scenarios. Phylogenetic comparative analyses of closely related lineages can facilitate evolutionary interpretations because unaccounted variation should be minimized. The Oklahoma salamander (Eurycea tynerensis) exhibits life cycle variation among populations. Most populations have the ancestral biphasic life cycle, while others exhibit paedomorphosis [16,17]. These alternative life history modes are highly correlated with local streambed substrate [17,18]. Biphasic populations occur on compact substrates and bedrock, whereas paedomorphic populations dwell within the relatively large interstitial spaces between coarse gravel, which permits permanent access to subsurface water [17–19]. Populations of E. tynerensis show high degrees of genetic isolation [17]. As a consequence, populations are more closely tied to the same selective regime (stream habitat), and local adaptations are not diluted by gene flow [20–22].

Here we develop a population-level phylogeny of E. tynerensis based on 190 nuclear genomic regions derived from anchored hybrid enrichment, and use it as a framework to test patterns of trait evolution. Specifically, we compare the fit of alternative evolutionary models to test whether salamander life cycle shifts influence the evolution of an important feature of vertebrate locomotion, the trunk skeleton. We test whether metamorphosis and a terrestrial life stage constrains the rate of evolution of the number of vertebrae in the trunk skeleton of E. tynerensis. This hypothesis predicts that shifts to a permanently aquatic lifestyle should remove terrestrial constraints and increase the rate of vertebral column evolution compared to biphasic lineages. Alternatively the specialized, permanently aquatic lifestyle of paedomorphic populations may have reduced rates of vertebral column evolution due to strong selection for a body form that can navigate within the interstitial spaces of a streambed. This study serves as an example of how changes in life cycle complexity can result in rapid phenotypic evolution.

2. Material and methods

(a) Specimens and vertebral column data

We examined the vertebral columns of 525 cleared and stained E. tynerensis from 28 populations (14 biphasic and 14 paedomorphic; figure 1; electronic supplementary material, table S1). Biphasic populations are more widespread across the Ozarks, and our sampling included representatives from throughout the range (figure 1), major genetic lineages [17], and most of the independent transitions to paedomorphosis. We only included 14 biphasic populations in our analyses in order to maintain balanced sampling of regimes (life cycle modes), which can be helpful when comparing the fit of trait evolutionary models [23]. At most locations E. tynerensis only exhibit a single life cycle mode [17,18]; however, our sample includes one locality (9) where paedomorphic (p9) and biphasic (b9) lineages co-occur in sympatry with little to no gene exchange.

• Download figure
• Open in new tab
• Download PowerPoint

Our samples consisted of 120 previously cleared and stained museum specimens, and 405 of our own preparations of museum specimens or material that we collected from 2001 to 2015. Salamanders we collected were anaesthetized by submersion in a 0.5% solution of tricaine methanesulfonate (MS-222), and tissues were collected and stored at −80°C as genetic vouchers. The remainder of the specimens were fixed in 10% buffered formalin. Soft tissues of preserved specimens were cleared and skeletons were stained generally following a procedure adapted from [24]. Formalin-fixed specimens were hydrated in distilled water, and cartilage was stained with an ethanol/glacial acetic acid solution containing alcian blue 8XG. Tissues were then macerated (cleared) in a saturated solution of sodium borate and trypsin. Bones were stained in a solution of 0.5% KOH solution containing alizarin red. The specimens were de-stained through a graded series of KOH/glycerol solutions, and ultimately stored in 100% glycerol.

The number of trunk vertebrae between the atlas and sacrum was counted at least twice using compound microscopes with bottom illumination (figure 2). The atlas of a salamander is a single cervical vertebra and is immediately posterior to the skull. The sacrum is also typically a single vertebra in salamanders and is located in the pelvic region, clearly identifiable by the presence of enlarged sacral ribs. Many salamander species have a low frequency of ‘double sacra’, where two adjacent vertebrae in the pelvic region each carry a single sacral rib on opposite (lateral) sides [25]. We observed this condition in less than 1% of specimens (five out of 525), which were removed from the analyses. At least 10 specimens were examined per population, and all analyses were performed on population means (with standard errors) and population modes.

• Download figure
• Open in new tab
• Download PowerPoint

Temperature is known to influence the number of trunk vertebrae in many ectothermic vertebrates, including salamanders [25,26]. Paedomorphic and biphasic populations of E. tynerensis occur in streams with different thermal regimes. However, based on the seasonal cycle of secondary sex characteristics and patterns of vitellogenesis (in the field and laboratory), E. tynerensis breed during late summer with oviposition occurring in the fall. Ozark streams generally have similar thermal regimes in the fall (ranging from approximately 16 to 20°C; Treglia et al. 2014–2015, unpublished data). The number of trunk vertebrae is a product of embryonic somitogenesis and specification of the pelvic region [25,27]. In E. tynerensis this process occurs less than a month after oviposition, and therefore the embryos of both paedomorphic and biphasic populations are likely to be experiencing similar temperature regimes during vertebral column development in the wild. Nevertheless, we examined 134 specimens that were raised from oviposition under common garden conditions to evaluate whether our observed patterns of wild collected E. tynerensis are probably inherited or simply an effect of differences in temperature or other factors experienced in the field. Three to seven mothers from each of three biphasic and three paedomorphic populations (electronic supplementary material, table S2) were mated in the laboratory with multiple males from their same wild population. Offspring were raised from the time of oviposition under a common temperature in the laboratory (18 ± 2°C).

(b) Phylogenomics

We used anchored hybrid enrichment [28–30] to isolate 190 nuclear genomic regions from representative individuals of each of the 28 populations (lineages) of E. tynerensis. Our anchored sequence capture protocol was described in detail in [30]. Briefly, genomic DNA was extracted using Qiagen DNeasy extraction kits, and sheared using a Covaris M220 focused-ultrasonicator. Illumina adapter indices were added to the fragmented DNA following Illumina TruSeq DNA LT Protocol. The concentrations of indexed DNAs were determined using a Qubit and equal proportions were pooled for sequence capture. Of note, 1365 unique 100 bp RNA MyBaits (MYcroarray, Ann Arbor, MI) were used to enrich several hundred nuclear exonic regions from across the genome (see [30] for details). Libraries were sequenced with 300 or 500 cycle version 2 cartridges on an Illumina MiSeq at the University of Tulsa.

CLC Genomic Workbench was used to map non-duplicate Illumina reads of each sample to 190 genomic regions (templates). These regions included one or more exons and intervening/flanking intron sequences, had high capture efficiency, high data coverage, and clearly captured orthologous genomic regions across all specimens. Positions with greater than 5× coverage were included in the alignment and positions with an alternative nucleotide of greater than 25% were scored with ambiguity codes (likely heterozygous positions). These totalled 133 241 aligned base positions (bp; average approximately 700 bp per genomic region), with an average coverage of 67% (approx. 89 300 bp) coverage per sample.

The phylogeny and divergence times (branch lengths) were estimated using a concatenated approach in BEAST v. 2.4 [31]. While patterns of population divergence are best estimated by summarizing individual phylogenies of unlinked genes [32–34], concatenated approaches can be effective when phylogenetic signal is low [34]. In our dataset many of our 190 genomic regions had relatively low variation (few variable sites), which challenges the reconstruction of individual gene trees. Nevertheless, phylogenetic analyses of more than 10 subsets of 100 randomly sampled genomic regions (approx. half of the dataset) produced largely congruent phylogenies, especially with respect to relationships among geographically proximal populations. Analyses were based on uncorrelated lognormal molecular clocks and constant population size priors. There are no known fossils of E. tynerensis, therefore we used a previously estimated date for the deepest divergence among populations of this taxon (approx. 5 million years) [12]. We applied a normally distributed calibration prior for the coalescence of all populations with a mean of 5 million years ago and standard deviation of 1 million years, which reflects the distribution of divergence time estimates for this node [12]. Analyses are based on relative branch lengths of the chronograms, and would be the same regardless of the overall time scale. The analysis was run for 10 million generations with trees saved every 10 000 generations (total 1000 trees). Stationarity of likelihood values was evaluated in Tracer v. 1.5 [35] and the first 50% of generations (500 trees) were conservatively discarded as burn-in, which was well beyond stationarity. We used a posterior distribution of 500 post-burn-in trees for comparative analyses.

(c) Modelling vertebral column evolution

Phylogenetic comparative methods are recognized as critical procedures when testing for relationships among traits or traits and the environment [36]. Gene exchange among populations can violate the assumptions of independence when applying comparative methods within ‘species’ [37]. However, while our 28 populations are considered part of a single taxon (E. tynerensis) they show strong geographical genetic structure with limited gene flow [17]. Even though this taxon is comprised closely related lineages, it is still important to account for relatedness when testing for evolutionary correlations [36–38].

Phylogenetic analyses of variance (phyloANOVA) [39] in the R package phytools [40] were used to test whether alternative life cycle modes (biphasic versus paedomorphic) exhibit significant differences in the number of trunk vertebrae among populations while considering phylogenetic relatedness. phyloANOVAs were performed using both means and modes of trunk vertebral counts for each population (electronic supplementary material, table S1).

We used OUwie v. 2.1 [23] to test whether rate of evolution (σ²) and optimum (θ) number of trunk vertebrae differ between biphasic and paedomorphic lineages of E. tynerensis under Brownian motion (BM) or Ornstein Uhlenbeck (OU) processes [23,41]. Differences in these parameters among life cycle modes were tested by comparing the following six models: (1) BM1: the simplest model based on BM with a single θ and a single σ² across the phylogeny; (2) BM_σ²: BM with a single θ but allowing for different σ²s for biphasic and paedomorphic regimes; (3) BM_θσ²: BM with different θs and σ²s for biphasic and paedomorphic regimes; (4) OU1: OU process with a single θ, σ², and selection strength parameter (α); (5) OU_θ: OU process allowing for different θs between paedomorphic and biphasic regimes, but with single α and σ²; and (6) OU_θσ²: OU process allows for different θs and σ²s, but a single α. Models with varying αs among regimes are also available in OUwie, but this parameter can be difficult to estimate accurately [23,42]. In order to reduce the number of parameters in our model comparisons we did not include these models (see also comments in the Results section).

Model comparisons were performed on 1000 stochastic character maps (make.simmap function) of life cycle mode (biphasic and paedomorphic) across the 500 post-burnin chronograms. We compared AICc scores (that account for small samples size) and AICc weights (w_i) [43] to assess the relative fit of models. Lower AICc indicates that the model is a better fit, and generally ΔAICc of less than two is considered an equal fit [43].

Support for differences in a given parameter between paedomorphic and biphasic lineages is based on substantial improvement in fit for a model that allows the parameter to differ between the selective regimes. For example, if the rate of evolution differs between paedomorphs and biphasics then at least one of the models that allows for differences in rate (BM_σ², BM_θσ², OU_θσ²) among regimes should be a substantially better fit (lowest AIC) than all models that estimate a single rate across the tree (BM1, OU_θ). 95% confidence intervals (CIs) were calculated for parameters of the best-fitting models.

We also evaluated the validity of parameters through simulation using the ‘multiOU’ function in phytools [40]. Using the parameters (θ, σ² and α) from the best-fitting model (OU_θσ²) we simulated 100 trunk vertebral datasets and fit the simulated data to BM and OU models using OUwie [23]. This was performed to assess the consistency of parameter estimates and relative fit of the best model compared with unsimulated trunk vertebral data.

The intentionally balanced sampling of paedomorphic and biphasic populations results in an ambiguous ancestral life history of E. tynerensis; however, when we consider that the most closely related Eurycea are biphasic (or ancestrally biphasic) then the ancestral population of E. tynerensis was also probably biphasic or polymorphic [44,45]. Therefore, we performed analyses considering three different root states: biphasic, paedomorphic or both life histories as equally probable.

3. Results

(a) Population-genomic-based phylogeny of E. tynerensis

The 190-nuclear phylogeny was largely consistent with previously published phylogenies based on portions of mtDNA and one nuclear gene [16,17] (figure 3; electronic supplementary material, figure S1). Prior analyses show three distinct genetic lineages from the eastern Ozarks (b6 and b8), southwestern Ozarks (b13 and b14) and western Ozarks (all other populations sampled here). A comprehensive phylogenomic analysis of E. tynerensis will be addressed elsewhere, but the analyses presented here based on 190 nuclear genomic regions (and assuming approximately 5 Ma coalescence) shows at least two major clades with several geographically structured subclades that are at least a few million years divergent. It also estimates that our sampling includes up to 20 lineages that are more than 1 million years divergent. Most importantly, the phylogenomic-based analyses agree with previous analyses (based on mtDNA) [16] that paedomorphic and biphasic populations are not reciprocally monophyletic, and include several transitions between life cycle modes during the Pleistocene (less than 2.5 Myr ago).

• Download figure
• Open in new tab
• Download PowerPoint

(b) Variation and models of trunk vertebral evolution

Paedomorphic lineages have on average 20.88 ± 0.469 trunk vertebrae (mode 21) while biphasic lineages have an average of 20.02 ± 0.135 (mode 20). The range of average trunk vertebral numbers is also almost 3.5 times greater in paedomorphic than in biphasic lineages (range: biphasic = 19.75–20.29 = 0.54; paedomorphic = 19.92–21.80 = 1.88; figure 3; electronic supplementary material, table S1). Phylogenetic ANOVAs of population mean numbers of trunk vertebrae show that these differences are significant (phylogenetic ANOVA: F_(1,26) = 43.11; p < 0.0001).

Offspring raised (from oviposition) under common laboratory conditions had statistically the same average and variance as their parents' wild population (lineage) indicating a heritable component to trunk vertebral number over at least one generation (electronic supplementary material, table S2, and figure S2). This includes paedomorphic lineages with both low (p3 = 20.15) and high (p13 = 21.05) average numbers of trunk vertebrae, as well as offspring of sympatric paedomorphic (p9) and biphasic (b9) lineages. We also found no differences in numbers of trunk vertebrae between males and females for either life cycle mode (electronic supplementary material, table S2 and figure S3).

Ornstein Uhlenbeck models that allow for different rates of trunk vertebral evolution and different optima between biphasic and paedomorphic life cycles were the best fit (OU_θσ²; table 1). These models substantially outperformed all other models (ΔAICc ≥ 15.6, w_i = 0.9995). Paedomorphic lineages were estimated to have a higher optimum number of trunk vertebrae (θ: biphasic = 20.02 ± 0.008; paedomorphic = 20.92 ± 0.021) and a dramatically higher rate of trunk vertebral evolution (σ²: biphasic = 0.007 ± 0.002; paedomorphic = 2.942 ± 0.197). Analyses based on data simulated under the parameters of the best-fit model (OU_θσ²_α) show patterns similar to analyses of the observed data (electronic supplementary material, figure S4). In other words, paedomorphs have a significantly higher optimum (θ: biphasic = 19.98 ± 0.011; paedomorphic = 20.88 ± 0.023) and rate of trunk vertebral evolution than biphasics (σ²: biphasic = 0.271 ± 0.066; paedomorphic = 3.051 ± 0.623). The estimated differences in rate based on simulated data are substantially reduced, but there is still an estimated 11-fold increase in rate of trunk vertebral evolution in the simplified paedomorphic regime.

We recovered the same best-fit model (OU_θσ²; electronic supplementary material, table S3) when we analysed population modes (instead of means). The only differences were that the selection parameter (α) was higher and the difference in rate (σ²) between paedomorphs and biphasics was even higher. The results were also the same regardless of whether we allow the ancestral life cycle mode of E. tynerensis to stochastically vary, set it to biphasic, or set it to paedomorphic (table 1; electronic supplementary material, table S4 and S5).

The strength of selection parameter was relatively high (α = 8.03 ± 0.532) [23], which may be expected given the strong shifts in optima between life cycle modes. Preliminary analyses showed that models that allowed for multiple αs (OU_θα and OU_θσ²_α) were a poor fit or did not substantially improve upon simpler models. There was also tremendous variation in α estimates across the stochastic character maps, and although α was consistently higher on average for biphasic lineages (indicating higher selection strength) the overlap in standard errors were large (e.g. α: biphasic = 5.17 ± 0.466; paedomorphic = 5.09 ± 0.468) and this difference was not significant (p = 0.4082).

4. Discussion

(a) Rapid evolution following changes in life cycle complexity

Complex life cycles are ubiquitous across metazoans [5,6], but patterns of phenotypic evolution immediately following changes in life cycle complexity have rarely been investigated. By examining vertebral column evolution among closely related salamander lineages that vary in life cycle complexity, we show that the phenotypic optimum and rate of evolution can respond rapidly to changes in number of life cycle stages. Paedomorphic (aquatic-only) lineages of E. tynerensis have a higher optimum number of vertebrae, and a dramatic increase in rate (variance) compared to lineages with biphasic (aquatic-to-terrestrial) life cycles (table 1 and figure 1). We also show that the number of trunk vertebrae in E. tynerensis is heritable (electronic supplementary material, table S2 and figure S2), reinforcing that the observed variance is a product of selection (as opposed to phenotypic plasticity).

Life cycle transitions have occurred repeatedly in E. tynerensis over the last five million years, with some shifts within the last one million years. The recency of life cycle transitions, the close association with vertebral evolution and also the high strength of selection estimate for vertebral evolution suggest that the phenotypic response has been relatively swift and often. Recent macroevolutionary analyses also showed that life cycle simplification was associated with dramatic increases in the rate of salamander body form and vertebral column evolution, with some of the most disparate forms observed in lineages that have been paedomorphic since the Mesozoic [12,46]. This echoes the patterns that we observed during the Pleistocene (less than 2.5 Myr ago) among lineages of E. tynerensis (figure 3). Taken together, our results provide evidence that phenotypic evolution can rapidly respond to changes in selective regime after transitions in life cycle complexity, and provides an early snapshot of macroevolutionary trends in salamander vertebral column evolution.

(b) Functional response to changes in ontogenetic complexity

The evolution of traits that are expressed across multiple life stages without significant transformation (metamorphosis) may be encumbered by functional requirements of each stage [1–4,12]. Therefore, the relaxation of constraints should offer opportunities to restore evolutionary flexibility. There is evidence of decreased performance when aquatic larvae and terrestrial metamorphosed salamanders locomote through opposing environments [7,47]. Metamorphosis may decouple salamander life stages, but such ontogenetic modifications have limitations. Across salamanders there have been clear shifts in both optimum and rate of phenotypic evolution between lineages with one-part versus two-part life cycles (this study) [12]. The higher rate of vertebral column evolution in paedomorphic lineages of E. tynerensis (and other obligately paedomorphic salamanders) is consistent with the hypothesis that simplification can remove constraints of a complex ontogeny. In other words, if biphasic lineages are constrained to trunk skeletons that are a compromise between the aquatic and terrestrial stages, then the removal of this constraint in aquatic-only paedomorphic lineages should result in shifts in optima and increases in variance of previously compromised traits.

Tetrapods exhibit extensive variation in trunk form [48–51]. Elongate trunk skeletons have repeatedly evolved in association with fossorial and aquatic ecologies [48]. Therefore, differences in patterns of vertebral column evolution between life cycle modes may reflect a functional response to changes in aquatic-only versus aquatic-to-terrestrial selective regimes. Ancestral salamanders exhibit a basic tetrapod body plan with a relatively short trunk containing approximately 13–17 vertebrae [12,52–54]. ‘Eel-like’ body forms occur in several paedomorphic families (amphiumids, sirenids, some proteids, batrachosauroidids and scapherpetonids), which have up to four times as many trunk vertebrae as the average biphasic species [12,46,55]. This illustrates the frequent association between salamander body elongation and permanently aquatic life cycles. Similar to higher-level patterns, we found that paedomorphic E. tynerensis show a significant increase in their optimal number of trunk vertebrae (table 1 and figure 3).

Understanding the primary selective force that drives the evolution of salamander elongation will require additional experimentation, ideally in a comparative context. The most cogent hypotheses involve locomotion and fecundity. Salamanders swim via anguilliform locomotion, which entails undulating the trunk and tail [56], similar to the locomotion of actinoperygean eels [57,58], so the evolution of a permanently aquatic salamander may select for a more ‘eel-like’ optimum number of trunk vertebrae [12]. Expanded numbers of trunk vertebrae occur in some (but not all) aquatic and terrestrial salamanders that burrow. In the case of E. tynerensis, paedomorphic populations dwell within the interstitial spaces of chert gravel [17–19]. Therefore, an elongated trunk may be just as advantageous for living in narrow spaces as it is for swimming [48,58].

The volume of the body cavity can limit the carrying capacity of yolked eggs or developing offspring [59], and in some species females have higher numbers of trunk vertebrae than males [25,60]. Based on analyses of scincid lizards, Griffith [61] hypothesized that trunk (and abdominal) elongation may evolve to maximize capacity for eggs or offspring, especially in diminutive species. Terrestrial slender salamanders (genus Batrachoseps) show distinct increases in presacral vertebrae of females, but no differences in body size between the sexes [25]. We found no differences in trunk vertebral numbers between males and females (electronic supplementary material, figure S3), but this still remains a viable hypothesis. The repeated transitions in life cycle among the closely related lineages of E. tynerensis are an excellent system for testing the selective pressures on tetrapod body forms following simplification to a permanently aquatic life cycle.

(c) Evolution of vertebral column elongation

Somitogenesis and elongation of the axial skeleton occur early in vertebrate development [27]. Relative lengthening of the trunk skeleton can arise through the addition or elongation of vertebrae [12,62,63]. These variables are somewhat correlated in salamanders, but some lineages have changed patterns of elongation without significant alterations to numbers of trunk vertebrae [12]. Perhaps one of the most unusual examples are neotropical salamanders (tribe Bolitoglossini), which mostly have a static number of trunk vertebrae (14), but display substantial variation in trunk form [12,62,63]. In this clade, worm-like body forms evolved through the addition (genus Oedipina) as well as elongation (some Pseudoeurycea) of trunk vertebrae [62,63]. The trunk vertebrae of E. tynerensis are proportionally similar among lineages, and the primary mode of elongation appears to be from the addition of trunk vertebrae.

It is important to consider the possibility that shifts in patterns of variation may result from pleiotropic interactions between somitogenesis and metamorphosis. For example, the loss of metamorphosis may break otherwise integrated molecular constraints [64–66], and as a consequence indirectly allow for an increase in variance of other structures such as the vertebral column. Another co-regulatory explanation is that variation in trunk vertebral numbers are governed along a molecular concentration gradient and metamorphosis is only achieved beyond some molecular threshold (e.g. [67]). However, at the individual level, paedomorphic E. tynerensis express a range of trunk vertebral numbers (19–22), which strongly overlap with those of biphasic individuals (19–21; figure 3). In fact, diverse trunk vertebral numbers are observed among paedomorphic lineages of Eurycea and paedomorphic salamanders in general [12].

(d) Permanence of life cycle transitions on trait diversification

Traits can exhibit a surprising degree of stasis over deep timescales despite a wealth of evidence that the same traits can rapidly respond to recent changes in selection patterns [68,69]. This is likely because patterns of selection may fluctuate around a central mean over time, which can ultimately limit major trait deviations. Pronounced changes in trait disparity often result from permanent transitions to new adaptive zones [70–72].

Environmental factors are associated with vertebral variation among some closely related salamanders [25,73,74]. However, biphasic salamander lineages have maintained relative stasis in numbers of trunk vertebrae since the Late Jurassic [12]. Lineages with environmentally regulated life cycles that can facultatively express paedomorphosis [75] have also exhibited relative stasis in numbers of trunk vertebrae [12]. By comparison, the most disparate body forms and trunk vertebral numbers occur in salamander lineages with obligately paedomorphic life cycles [12]. Obligately paedomorphic salamanders have lost the ability to transform into a fully terrestrial salamander, even when exposed to thyroid hormone [76,77], which otherwise induces amphibian metamorphosis. The loss of metamorphosis appears to be reversible over several million years [44,45], but some lineages such as the families Amphiumidae, Cryptobranchidae, Proteidae and Sirenidae have probably been paedomorphic since the Mesozoic [46,52,53]. Therefore, transitions to obligate paedomorphosis, as opposed to facultative paedomorphosis, likely represent long-term shifts to a completely aquatic adaptive zone.

Larval-form paedomorphosis in plethodontid salamanders is largely obligate [77]. Facultative life cycles do occur in E. tynerensis, but most populations appear to be obligately paedomorphic or obligately biphasic [77]. There is also variability in sensitivity to thyroid hormone among paedomorphic lineages [78]. Interestingly, paedomorphic lineages with reduced sensitivity to thyroid hormone (e.g. p13) have more trunk vertebrae than paedomorphic lineages with higher sensitivity (e.g. p3). Reduced sensitivity to thyroid hormone and more divergent numbers of trunk vertebrae may both represent a longer commitment to a completely aquatic lifestyle.

5. Conclusion

Our study shows that recent transitions in salamander life cycles result in substantial changes in patterns of phenotypic evolution. Lineages with a simplified aquatic-only life cycle have more trunk vertebrae and a higher rate of vertebral evolution, a pattern that reflects macroevolutionary trends across salamanders. If life cycle complexity constrains vertebral column evolution then our study serves as example of how relaxation of such constrains can result in a rapid shift in phenotype. Clades with variation in life cycle complexity offer powerful systems to study how changes in selective regimes during ontogeny shape patterns of phenotypic and molecular evolution.

Ethics

Salamanders were handled in accordance with IACUC protocols at the University of Tulsa (TU-0028 and TU-0029). Scientific collecting permits were issued by the states of Arkansas (108203252003162322, 031020055, 110520071, 041320091, 022720124, 021720143, 030220153, 022320164), Missouri (11516, 14242, 15997), and Oklahoma (1471, 3807, 4843, 5107, 5375, 5545–5548, 5960).

Data accessibility

DNA sequence alignments are available on Dryad digital repository, doi:10.5061/dryad.s02p0. Data summaries are listed in electronic supplementary material, tables S1 and S2.

Authors' contributions

R.M.B designed the study. R.M.B., J.G.P., S.D.M. and L.L. performed laboratory research and collected data. R.M.B., J.G.P., and N.M.L. analysed the data. R.M.B. wrote the manuscript with contributions and edits from all authors.

Competing interests

The authors declare no competing interests.

Funding

Funding for this research was provided in part by the University of Tulsa, Oklahoma Department of Wildlife and U.S. Fish and Wildlife Service Grant to R.M.B., and National Science Foundation Grants (DEB 1050322 and Oklahoma EPSCoR IIA-1301789) to R.M.B.

Footnotes

Electronic supplementary material is available online at https://dx.doi.org/10.6084/m9.figshare.c.3965352.