Genome size influences adaptive plasticity of water loss, but not metabolic rates, in lungless salamanders

Natural selection operates on individual organisms by acting on all their functional phenotypes simultaneously. Because many ecologically important phenotypes, e.g. physiological performance, are themselves complex integrations of multiple underlying traits, characterizing the functional role of component traits driving performance phenotypes can offer a simplified method to identify patterns of adaptation in nature (Adler et al., 2014). To do this, we must link variation in traits to functional phenotypes of known adaptive importance. One of the most fundamental performance phenotypes is metabolic rate, the sum total of all cellular respiration within an organism. Metabolism is a critical indicator of energy turnover that limits many other biological functions (Brown et al., 2004), making the functional contributions of traits to variation in metabolism of broad ecological importance (Glazier, 2018; Uyeda et al., 2017). The most well documented such trait is body size, which co-varies with metabolism across taxa (Brown et al., 2004; Gillooly et al., 2001), yet we know that variation in metabolic rate is influenced by factors beyond simply body mass (Glazier, 2005; White and Kearney, 2013). Thus, to fully understand mechanisms underlying physiological diversity will require integrating morphological and environmental variation influencing performance phenotypes to test the importance of functional diversity underlying adaptive performance in nature.

One potentially important trait influencing metabolism is cell size. Cell size correlates strongly with genome size (C-value) (Gregory, 2001a). Although the mechanism linking these traits remains unresolved, the prevailing hypothesis holds that C-value causatively determines cell volume, by extension exerting a nucleotypic effect on phenotypes related to cellular transport and differentiation rates (Gregory, 2001b). This effect on cellular physiology can scale to the level of the organism, for example C-value appears to be related to mass-specific metabolism in some homeotherms (Gregory, 2002; Vinogradov, 1995) and to cellular differentiation rate (Sessions and Larson, 1987) and duration of embryonic development in amphibians (Jockusch, 1997; Litvinchuk et al., 2007). Theoretical work further suggests that C-value moderates the allometric scaling of metabolic rate across taxa, with metabolism in lineages with variable C-values demonstrating shallower allometric slopes that deviate further from isometry than in lineages without C-value diversity (Kozlowski et al., 2003). Despite these associations, it remains an open question whether variation in C-value is itself an adaptive target of selection, or represents neutral variation subject to random drift (Mueller, 2015). Given the importance of performance phenotypes in limiting many biological processes (Kearney, 2012; Pettersen et al., 2019; Riddell et al., 2018), the possibility that C-value variation influences physiology highlights a need to identify potential mechanisms by which genome size contributes to observed patterns of organismal function and diversity in nature.

The lungless salamanders (Urodela: Plethodontidae) present an excellent system in which to test for effects of genome size on organismal performance. This lineage exhibits exceptionally broad variation in C-values (∼10–76 pg) (https://www.genomesize.com), offering substantial natural variation to test for signatures of selection on genome size. Critically, all respiratory gas and water transport in this group occurs by diffusion across the skin surface, making the structure and function of the skin strong targets of selection (Feder and Burggren, 1985). The obligately permeable and moist skin in these salamanders (Feder and Burggren, 1985; Shoemaker and Nagy, 1977) leaves individuals at high risk of desiccation, with selection favoring skin with a higher resistance to water loss (Riddell et al., 2018). Skin thickness (the diffusion path length) and composition (which determines diffusivity) influence the resistance of the integument to gas and water transport (Feder and Burggren, 1985; Lillywhite, 2006), and in this capacity C-value may influence resistance. All else being equal, larger cells in high C-value taxa will produce thicker skins, resulting in lower transport rates (Feder and Burggren, 1985; Piiper, 1988), or will otherwise necessitate compensatory changes in skin morphology to maintain adequate transport (Lillywhite, 2006). Although resistance must remain generally low to enable cutaneous gas exchange in this group (Feder, 1983), relatively higher resistances hold clear adaptive value in reducing the risk of desiccation (Riddell and Sears, 2015; Riddell et al., 2018). The potential link between C-value and resistance to water loss in the Plethodontidae remains untested, and presents a clear test case for a potentially adaptive role of genome size variation in nature.

Here, we explore the consequences of C-value diversity on physiological performance in the lungless salamanders, testing for physical effects of C-value at the individual level, which may demonstrate an adaptive role of C-value in cutaneous gas and water transport across species. We quantified body and genome size variation in eight species to test for their combined functional contributions to performance under different environmental conditions. Using a combination of linear models, comparative methods and path analyses, our work explored the question: is C-value a functionally relevant trait in plethodontid physiology, and if so, what are the mechanisms and environmental context through which it moderates physiological performance? We hypothesized two potential mechanisms by which C-value can influence plethodontid physiological rates. First, a general nucleotypic effect of cell-specific metabolism may result in lower metabolic rates in larger C-value individuals for a given body size (Gregory, 2002). Second, C-value may moderate diffusion lengths for water evaporating out of the body, resulting in adaptive differences in resistance to water loss. By integrating morphological and environmental variation underlying critical performance phenotypes, this work takes an explicit trait-based approach to test the importance of C-value diversity underlying adaptive performance in nature.

All field, laboratory and experimental techniques described here were approved by the Cornell University Institutional Animal Care and Use Committee (IACUC protocol 2016-0081). We collected eight species of salamanders [Eurycea bislineata (Green 1818), Eurycea longicauda (Green 1818), Plethodon cinereus (Green 1818), Plethodon dunni (Bishop 1934), Plethodon glutinosus (Green 1818), Plethodon larselli (Burns 1954), Plethodon vandykei (Van Denburgh 1906) and Plethodon vehiculum (Cooper 1860)] from multiple natural populations (Table S1). Plethodon are direct-developing terrestrial species, and Eurycea are semi-aquatic with a biphasic life history (aquatic larvae that metamorphose into semi-aquatic adults). These species range broadly in body mass from P. cinereus (mean 0.71 g in our dataset) to P. glutinosus (mean 5.03 g), as well as substantial variation in C-value (range 22.64–69.3 pg). Mass and C-value are not co-linear in this sample (P=0.936), allowing us to test for their interaction. We maintained animals in the laboratory at two temperatures (10 and 15°C), housed individually in deli cups on a sponge saturated with dechlorinated water and fed crickets weekly ad libitum. Animals were acclimated for a minimum of 1 week before inclusion in physiological experiments, and were measured a minimum of 4 days after eating or a prior trial to ensure all measurements were taken in a minimally stressed, post-absorptive state.

Experimental trials consisted of repeated measurements over 30 min. We removed animals from their housing containers and placed each inside the respiration chamber constrained within a steel mesh tube oriented facing the incurrent airflow. This ensured all animals were immobile with the entire body surface area exposed to the airstream during steady-state measurements, effectively controlling for airflow and boundary layer effects contributing to water loss rates. Immediately after placing an animal in the chamber, we set the chamber humidity to the appropriate VPD and recorded atmospheric conditions inside the chamber every ∼40 s, continuously monitoring the respiration rate. In all cases, respiration dropped continuously (range 5–64 min, median 13.4 min) before stabilizing at a resting rate. After 5 min of stability (determined visually), we continued the trial for an additional 25 min. We saved only stable measurements over this 30 min period for later analyses, additionally discarding all measurements that fell outside our ±0.015 kPa boundaries. To test for environmental effects on physiological performance, we conducted four experiments per individual, factorially designed across two temperatures (10 and 15°C) and two VPD treatments (0.3 and 0.6 kPa). We randomized treatments with respect to temperature first and then VPD to prevent order effects and maximize randomness while avoiding repeated re-acclimation between every trial.

We tested for the influence of environment and morphology on physiology using linear mixed effects models in R (http://www.R-project.org/). We tested body mass and C-value (collectively referred to as morphological variables within our analyses) plus temperature and VPD against six physiological responses: gross and specific metabolic rate, gross and specific EWL, CWL and r. For total resistance models we included mass as a covariate in all cases. Total resistance is the sum of two component resistances: the resistance of the integument itself and the boundary layer, a layer of relatively still air surrounding the animal. Because boundary layer resistances co-vary with body size (Gates, 1980), this approach effectively controls for boundary layer effects contributing to total resistance, and focuses on effects attributable to physiological variation in cutaneous resistance.

We followed a combination of model building and model simplification, first testing environmental and morphological variables separately to survey a priori any potential effects of C-value. For our environmental treatments, we tested four models for each physiological response: the interaction between temperature and VPD (Temp×VPD), both effects in tandem (Temp+VPD), and each predictor independently. All environmental models were tested across our complete dataset, which includes multiple measurements per individual across treatments. To account for multiple testing, we included individual as a random effect in all models. Because we are interested in physiological differences between morphological types, and not between species per se, we included species as a second random effect in all models. This approach thus focuses not on the effect of phylogeny on performance or on the effect of morphological predictors through evolutionary time, but rather on the physical effects of body mass, C-value and environmental variation on organismal physiology at the individual level.

To identify environmental conditions under which C-value may be functionally relevant, we tested morphological models at four treatment levels: across the complete dataset, at each temperature, at each VPD, and within each treatment (Temp×VPD) combination. Across the full dataset, we tested each response against morphological predictors while controlling for treatment combination, plus random effects for species and individual. Models tested at each temperature or VPD included random effects for species and individual while controlling for VPD or temperature, respectively. When tested within treatment combinations, models used only a random effect for species. For each response–subset combination, we compared four models: the interaction between body size and genome size (Mass×C-value), both in tandem (Mass+C-value), and each independently.

For each response variable, we constructed a more complex full dataset model incorporating all significant predictors from our morphological and environmental models, plus random effects. We selected a final model for each response by removing non-significant effects from these larger models through stepwise model simplification. These final models are presented in Table 1, with the complete table of all tested models and model results in Table S2. We tested a total of 216 morphological and 24 environmental models before constructing final models. We assessed each model's support by calculating the corrected Akaike information criterion (AICc) of each and comparing values among nested models (Table S2).

While our linear models account for environmental and individual-level variation contributing to trends in physiological performance, this approach cannot account for phylogenetic history in contributing to performance. To control for phylogenetic non-independence among species and analyse our performance data in a comparative framework, we conducted phylogenetic generalized least squares (PGLS) regressions. We used DNA sequences from previously published studies (Fisher-Reid and Wiens, 2011; Newman et al., 2016) to construct a multi-locus phylogeny of 60 plethodontid species, plus two outgroup species of Amphiuma. A similar phylogeny for 55 of these species was previously published by Newman et al. (2016). We used all sequences used by Newman et al., plus homologous sequences for P. dunni, P. larselli and E. longicauda. Because E. bislineata was used as an outgroup by Newman et al. but is a focal taxon in this study, we included sequences from Amphiuma means and Amphiuma tridactylum as new outgroups. All sequences were downloaded from GenBank (Clark et al., 2016), aligned using MAFFT version 7 (Katoh and Standley, 2013), and concatenated for a final alignment of 4078 nucleotides over six nuclear loci: BDNF (707 bp), GAPD (664 bp), ILF3 (281 bp), MLC2A (253 bp), POMC (481 bp), RAG1 (1467 bp) and RHO (225 bp). A complete list of sequences is available in Table S10 in Dryad (dryad.t76hdr7zh). We assembled our phylogeny in a Bayesian framework using MrBayes version 3.2.7 (Ronquist et al., 2012), employing the same sequence partitions, models of nucleotide evolution and software settings as used previously by Newman et al. (2016) and Fisher-Reid and Wiens (2011). We conducted two runs of 1 million generations using four chains sampled every 1000 generations. After discarding the first 100,000 samples as burn-in, we assessed convergence in Tracer version 1.7.1 (Rambaut et al., 2018) by ensuring effective sample sizes were above 200 before accepting the consensus tree.

For PGLS analyses, we trimmed taxa from the consensus tree to produce a phylogeny of just our eight focal species. We calculated mean predictor and response values for each species across our entire dataset and at each experimental subset, using subset-specific trait values for each model. In all cases, we included a covariance structure based on Brownian evolution across the trimmed phylogeny and fitted models by maximizing the log-likelihood. For morphological models, we used C-value estimates as in our linear models, and mean body mass for each species from our experimental measurements. Because PGLS analyses allow only one measurement per species (Garamszegi, 2014), this approach cannot make within-species comparisons between our experimental treatments. Thus, for environmental models, we substituted experimental treatments with climate profiles generated for each species by randomly sampling geographic coordinates within each species distribution and extracting environmental data from the sample. We extracted mean annual temperature (BIO1) from WorldClim (http://www.worldclim.org; Fick and Hijmans, 2017) and aridity index (AI) from CGIAR-CSI (http://www.cgiar-csi.org; Zomer et al., 2008), extracting one sample per 400 km² [range 39 (P. larselli) to 4722 (P. cinereus)], and used the mean of extracted values for each species in our model dataset. All geospatial processing was conducted in the R packages ‘raster’ (https://CRAN.R-project.org/package=raster) and ‘sp’ (Bivand et al., 2013) using a 30 arc-second (∼1 km²) spatial resolution.

Using our phylogeny, we conducted PGLS regressions in R using ‘ape’ version 5.0 (Paradis and Schliep, 2019). We constructed PGLS models equivalent to each of our morphological linear models, substituting experimental temperature and VPD treatments with species-level climate profiles for BIO1 and AI. While our linear models test direct environmental effects on individual-level variation in performance, these PGLS models test for the contribution of the natural home environment to each species’ performance, as well as evolutionary associations between morphological variables and physiological responses. As above, we constructed final models across the complete dataset with all significant morphological and environmental effects, reaching a final model through a last round of model simplification. We tested 216 morphological and 24 environmental models before constructing final models, evaluating model support by AICc comparison.

We constructed further linear mixed models to test for effects of C-value diversity on the relationship between body size and gross metabolic rates. Species in the genus Plethodon vary substantially in genome size, in contrast with minimal C-value diversity in Eurycea (http://www.genomesize.com). To test if these differences result in different allometric scaling of metabolism among genera (Kozlowski et al., 2003), we modelled metabolic rate as a function of the interaction between mass and genus using the R package ‘lme4’ (Bates et al., 2015). We tested models across our entire dataset, and at each temperature and VPD, using species and individual as random variables in all models. Because our linear models found no interaction between temperature and VPD (see Results), we did not test individual treatment combinations.

We tested for differences between genera in both the slope and variance around the slope [‘weights=varIdent(form=∼log(mass)|genus)’]. To test for differences in allometric scaling among genera, we compared model slopes between genera using the ‘emtrends’ function in the R package ‘emmeans’ (https://CRAN.R-project.org/package=emmeans). To test whether C-value diversity contributes significant variance to the relationship between mass and metabolism, we constructed a second model per data subset without genus-specific variances, and compared each model pair using likelihood ratio tests. We assessed parameter significance of each by correcting for the false discovery rate within each response variable (five model pairs each). Although the theoretical model tested here (Kozlowski et al., 2003) applies specifically to metabolic rate, shared physical constraints on gas and water transport in plethodontid salamanders and the collinearity of these two responses in our dataset (P<2×10⁻¹⁶) make it reasonable to expect a similar effect of C-value on water loss rates. Thus, we investigated the allometry of water loss by repeating the same analyses on allometric scaling and variance using EWL as a response.

We conducted path analyses to quantify mediation effects among our morphological and physiological variables. These models tested three responses: gross, specific and cutaneous water loss rate. We tested each response under each VPD and temperature treatment, and across our complete dataset. In each case, models had the same general structure: mass had a direct path to the response variable and a path mediated through total resistance. C-value was also mediated through the resistance parameter but lacked a direct path to the response. While similar in concept to our linear models testing physical links between morphology and water loss rates, these path models explicitly address the mechanisms underlying these relationships. This path structure tests the hypothesis that C-value may influence resistance to water loss through its impact on skin morphology, and quantifies the relative contributions of both mass and C-value as predictors of water loss rates, both directly and as mediated by resistance.

We log-transformed all variables, testing a total of 15 path models across all responses and treatments. All models were constructed and analysed in the R package ‘lavaan’ (Rosseel, 2012) using a maximum likelihood estimator with robust standard errors. We assessed the significance of path parameters using a Bonferroni correction for each response (five models each). Because lavaan only provides significance estimates to three decimal points, we manually assigned all estimates of 0.000 a value of 0.00049 (the largest value that rounds down to 0.000) before P-value correction to provide a maximally conservative estimate of significance. Full descriptions of all models can be found in Table S6 in Dryad (dryad.t76hdr7zh). For each response, we compared the total effects of mass and C-value between temperature and VPD treatments (Table S7 in Dryad: dryad.t76hdr7zh), again with modified P-values and Bonferroni corrections between treatment combinations (two comparisons). Because all path predictors and responses are on a log scale, we interpreted model slopes in terms of percentage change (Yang, 2012).

To investigate temperature effects evident in our linear models (see Results), we compared r values between 10 and 15°C using a t-test across all individuals and for each species. Because the different life histories of our focal taxa (terrestrial Plethodon versus semiaquatic Eurycea) are likely to influence their physiological resistance to water loss, we also compared r thermal dependency for each genus. We quantified the temperature dependence of r using the thermal ratio (R₁₀) for r following Bennett (1984). This metric quantifies the degree of change in a response with changing temperature, where R₁₀>1 is a positive thermal dependence, R₁₀<1 is a negative dependence, and R₁₀=1 indicates no effect of temperature. We calculated one R₁₀ per group, using group mean r values at 10 and 15°C. We then used PGLS to test for an effect of C-value on temperature dependence by regressing C-value and R₁₀ while controlling for mass (Table S9 in Dryad: dryad.t76hdr7zh). We ran models across the full dataset, and at each VPD. To test for a signature of adaptive functionality of C-value, we ran two PGLS models per subset: one assuming Brownian motion (BM, random drift), and another assuming an Ornstein–Uhlenbeck evolutionary model (OU, adaptive evolution), and compared each model pair using a likelihood ratio test.

In our environmental models, temperature and VPD significantly predicted all physiological responses in tandem, but without an interaction effect between them (Table S2). Among our morphological models, when tested across the entire dataset, body mass was the best supported predictor for all responses. Mass best predicted gross and specific respiration across all treatments (Fig. 1), and best predicted the majority of water loss models (gross and specific EWL, CWL), with occasional exceptions. C-value showed no significant effect on metabolism, or on gross or specific water loss across the full dataset; however, it predicted CWL and total resistance to water loss at 15°C (Fig. 2). For both these responses, C-value alone was the best supported model when tested at either VPD combination at 15°C, but tested across both VPD treatments we found equivocal evidence for an interaction between C-value and mass, as C-value and Mass×C-value were statistically indistinguishable by AICc. When compared with our separate environmental and morphological models (Table S3), final models were the best supported model for all responses (Table 1). All top models except total resistance included mass, temperature and VPD as predictors. For r, the best supported model included the interaction between mass and C-value plus temperature and VPD.

Body mass was the best supported morphological predictor for all PGLS models when tested across the complete dataset (Table S4). Similar to our linear models, we also found evidence for an effect of C-value on CWL and r at 15°C (Fig. 2). For gross and specific EWL at 15°C, the best supported model was either Mass+C-value in tandem or Mass alone; however, this was indistinguishable from Mass+C-value and occasionally Mass×C-value. For r at 15°C, C-value alone was best supported. For environmental models, mean annual temperature was the best supported predictor for all responses except specific respiration, which was best predicted by aridity index. Final models were variable among responses (Table S5). When compared with morphological and environmental models, final models were either the best supported or indistinguishable from the best supported model for all responses except gross metabolic rate, which was best predicted by mass alone (Table 1).

We found no difference in scaling coefficients for metabolism between Plethodon and Eurycea, regardless of data subset (P=0.429–0.870) (Table S6 in Dryad: dryad.t76hdr7zh; Fig. S1). Gross water loss also scaled equivalently between genera (P=0.265–0.731). Variance around the slope was consistently higher in the diverse C-value genus Plethodon for all metabolism models (mean ratio: 1.241/1), but only one of five water loss models (mean ratio: 0.989/1). Following a P-value correction, none of these differences was significant for either metabolic rate (P=0.100–0.731) or water loss (P=0.675–0.813).

C-value had a positive effect on total resistance in all path models (Table S7 in Dryad: dryad.t76hdr7zh). Across our full dataset, this amounted to a 0.11% increase in r per 1% increase in C-value (P=0.025). The significance of this effect varied across treatments; C-value consistently predicted total resistance across the full dataset (P=0.025), but never at 10°C (P=1) or 0.3 kPa (P=0.175). Mass had a positive effect on r in all cases (0.14% across the full dataset), consistent with covariation between mass and boundary layer resistance. For all responses, this influence was significant across our complete dataset (P=0.002), but when divided among treatments was only significant at 10°C and never at 15°C.

C-value had a negative effect on all water loss responses (mean=−0.15%) (Fig. 3), but significance of this effect was highly variable across treatments. C-value predicted each response at 15°C and 0.6 kPa, but never at 10°C or 0.3 kPa. This effect remained significant across the complete dataset for all responses (P=0.025). The total (direct+indirect) effect of mass was positive on gross water loss rates (0.55%) and negative on specific (−0.45%) and cutaneous rates (−0.15%) (Fig. 3B). This effect was consistently significant across the complete dataset (P<0.003), but was not significant for CWL at 15°C. For all water loss responses, neither mass nor C-value showed significant differences in total effect when compared between temperatures (P=0.064–0.074) or VPD treatments (P=1) (Table S8 in Dryad: dryad.t76hdr7zh).

Total resistance values increased with temperature (R₁₀>1) for all species except P. glutinosus (Table 2; Fig. 4A). Across the complete dataset, r was significantly different between temperature treatments (t=3.95, P=0.0001), but this trend was driven by patterns in Plethodon. When split by genus, r was significantly different between temperatures in all subsets for Plethodon, but never for Eurycea. All Plethodon except P. glutinosus showed a significant difference in at least one VPD subset, but this trend held across the full dataset for only three species (P. cinereus, P. larselli and P. vandykei). Our PGLS models found an effect of C-value on R₁₀ (Fig. 4B; Table S9 in Dryad: dryad.t76hdr7zh) across the full dataset using an Ornstein–Uhlenbeck covariance structure, but this model was not distinguishable from Brownian motion.

Organisms are highly integrated biological units, and the functional interactions that underlie their performance are extremely complex. In this study, we tested for a role of genome size contributing to physiological performance in lungless salamanders, for which we find support under specific conditions. While we found no relationship between C-value and metabolic rates or any effect of C-value variation on physiological allometry, our data do demonstrate that larger genomes are associated with reduced water loss in lungless salamanders exposed to warmer temperatures – a potentially adaptive function of genome size that has not previously been reported. Although identifying a consistent model for this relationship was challenging, our results repeatedly demonstrated that C-value – in some capacity, whether interacting with mass, in tandem with mass, or acting alone – exerts an effect on water loss physiology at 15°C. Discerning a mechanism for this effect requires a careful accounting of our analyses and synthesis with established theory and empirical evidence.

Any effects of genome size on physiology are clearly secondary to the influence of body size. Consistent with a robust literature documenting the importance of body mass to metabolic rates (Brown et al., 2004; Gillooly et al., 2001) and water loss (Peters, 1986; Spight, 1968) across taxa, mass best predicted all physiological responses in all linear and PGLS models across our full dataset, and consistently exerted a greater total effect than C-value in our path models. Unsurprisingly, the influence of body mass in plethodontid physiological performance is clear. How genome size moderates this effect is more complicated to assess. C-value appears to play no role in moderating physiological allometry in these taxa (Fig. S1; Table S6 in Dryad: dryad.t76hdr7zh), as both variance and allometric slopes for gross respiration and water loss were statistically equivalent between Plethodon and Eurycea despite their differences in C-value diversity. We found no effect of genome size on metabolic rates either, as models containing C-value had universally low support across linear and PGLS models. This result is consistent with previous studies demonstrating no relationship between genome size and metabolism in salamanders specifically (Gregory, 2003; Licht and Lowcock, 1991) and contributes to an emerging consensus that these traits are unrelated across vertebrates generally (Gardner et al., 2020; Uyeda et al., 2017), and suggests that C-value does not exert a nucleotypic effect on metabolism in this lineage as we hypothesized.

In contrast to metabolism, water loss physiology appears to be influenced by C-value in this group. Gross and specific EWL, CWL and r all showed support for an effect of C-value at 15°C, although in all cases our linear models lacked a clear best model. Multiple models containing C-value were nonetheless better supported than models for mass alone, suggesting a functional role of C-value in driving water loss rates under warm conditions. Despite this, only our top linear model for r had a C-value parameter, and our top PGLS models all lacked C-value. Together, these results imply a potentially subtle functional role of C-value that is not universal across physiological responses, environmental conditions or phylogeny. Instead, larger C-values may play a role specific to moderating water loss that: (1) acts specifically via resistance to water transport and (2) is environment dependent, with greatest effect under warmer conditions. This inference is supported by both our linear and PGLS models at 15°C. Although our path models also found a significant r-mediated effect of C-value on all water loss phenotypes across the full dataset, as all three model sets agreed in identifying an effect at 15°C, we tentatively infer that larger C-values help reduce water loss specifically at higher temperatures.

Our analyses raise the possibility that C-value variation plays a functional component role in plastic increases in r at high temperatures. Resistance increased with temperature in our experimental trials (Fig. 4A), coincident with decreases in all our water loss responses, a pattern indicative of plasticity in r across temperatures. Given the selective strength of desiccation risk in plethodontids (Baken et al., 2020; Feder, 1983) and the ecological importance of minimizing water loss (Feder, 1983; Riddell et al., 2018), this result suggests that plasticity in r is adaptive at higher temperatures (Ghalambor et al., 2007). The significant mediation of C-value in our path models at 15°C thus implies that genome size mechanistically influences increased resistance to water loss under warmer conditions. Our analyses cannot address the specific mechanism underlying this relationship, but the biophysics of resistance indicates one possibility. Total resistance is the sum of two component resistances in series: the resistance of the integument to water transport (r_i), and that of the boundary layer surrounding the animal (r_b). r_i is a biophysical and physiological property dependent upon skin structure and composition, and r_b is a physical property dependent on environmental conditions (largely air flow), body mass, shape and orientation to air flow (Gates, 1980). Thus, any strictly physiological (non-behavioral) control of resistance will by necessity be mediated through the cutaneous (r_i) component of r. Because we constrained our study organisms to prevent behavioral adjustment of r_b, after controlling for differences in body mass, measured variation in r can be attributed to differences in cutaneous resistance. Our path models demonstrate a significant relationship between r and C-value, therefore suggesting a mechanistic link between them within plethodontid skin.

Recent work has revealed the importance of temperature-induced adaptive plasticity in r_i in plethodontids (Riddell et al., 2018; Riddell et al., 2019), but how C-value influences this phenotype has until now been unexplored. The apparent mechanism for plastic resistance is by regulating cutaneous perfusion, in which vasoconstriction and the eventual breakdown of skin capillaries shunts blood away from the body surface to minimize water loss to the external environment as both temperature (Riddell et al., 2019) and dehydration severity increase (Brown, 1972). The effects of C-value we observed at our higher experimental temperature thus raise the possibility that C-value is only functionally relevant in less-perfused skin. Thus, we hypothesize that under stressfully warm conditions, as capillaries recede from or vasoconstrict in the integument, skin composed of larger cells may provide a defense of last resort against further rapid desiccation.

Theoretical and empirical work in amphibians support our hypothesis. Theory predicts that thicker skins will reduce mass transfer rates across the integument (Piiper, 1988), while denser capillary networks are empirically associated with higher transport rates (Malvin, 1988). Capillary recruitment influences ‘effective’ skin thickness, as vasoconstriction in the dermis increases the distance between the body surface and the nearest perfused capillaries exchanging material with the environment (Lillywhite, 2006). Thus, under normal skin perfusion, high water content in the capillaries of the skin should promote relatively high rates of water loss through the integument (Burggren and Vitalis, 2005). Under these high fluxes, C-value variation between species may produce negligible differences in their capacity to resist water loss. However, as cutaneous capillaries constrict to increase the skin's effective thickness, flux rates will depend more upon the diffusion of water from skin cells themselves. Under these conditions, larger cells with higher C-values may more effectively limit water transport through the integument before evaporation from the body surface, either due to their lower diffusion rates (Goniakowska, 1973) or their effect on skin composition (Lillywhite, 2006). Skin thickness (the depth of the dermal layer) per se therefore may not underlie transport rates as we originally hypothesized, but rather water loss may be constrained by an interaction between effective skin thickness and cell size in tissues between peripheral capillaries and the body surface. While C-value is clearly not the primary driver of variation in water loss rates, we hypothesize that this mechanism may play a role in reducing marginal evaporative loss for terrestrial plethodontid salamanders under stressful conditions.

Other factors besides effective skin thickness may contribute to observed water loss patterns. Water evaporating from the body surface may be lost after diffusing through skin cells, or from cutaneous glandular mucous secretions produced to avoid loss from the epidermis itself (Lillywhite and Maderson, 1988). Studies of variation in mucous secretion in salamanders are lacking, but in frogs mucous is predominantly secreted by species that bask in direct sunlight (Lillywhite and Licht, 1975). Given that lungless salamanders are restricted to moist microhabitats away from direct sun (Feder, 1983), mucous secretion is unlikely to be a meaningful contributor to water loss rates in our sample. Alternatively, lower water loss rates may result from higher lipid concentrations within, or secreted onto, the integument. Lipid secretions can enormously improve resistance to water loss, but also reduce cutaneous gas exchange (Lillywhite, 2006). Strict reliance on cutaneous respiration in the Plethodontidae thus makes lipid secretion another unlikely strategy for managing water loss. Variation in skin lipid content may contribute to variation in r, but the strength of this relationship is, to our knowledge, unexplored in this group. While we cannot rule out this possible mechanism, it is unclear how lipid content would plastically increase r at higher temperatures, or in species with larger C-values. Covariation between genome size and lipid content in the skin could still produce our observed results, as could differences in skin diffusing capacity if C-value influences cell packing and other structural changes under desiccation. However, without having quantified these traits and without strong literature in this family, the effects of these or other unmeasured potential covariates with C-value remain speculative. As it stands and given our data, our best hypothesis holds that, whatever the mechanism, observed plasticity in r is underlain by variation in genome size itself.

Selection on C-value for resistance to water loss is likely to be strongest in terrestrial environments. C-value correlates with cellular differentiation rate (Sessions and Larson, 1987) and development time (Jockusch, 1997; Litvinchuk et al., 2007) in salamanders, with strong selection in biphasic, semi-aquatic species, for whom ephemeral habitats favor fast development and associated small C-values to metamorphose before habitats become dry (Lertzman-Lepofsky et al., 2019). This timing constraint is absent in direct-developing terrestrial species, in which selection instead acts to minimize the risk of desiccation under exposure to air, making an association between C-value and r more likely in this group. We found this pattern in our R₁₀ data, where the terrestrial genus Plethodon consistently showed adaptive increases in r with temperature (R₁₀>1 for five of six species), but biphasic Eurycea showed none (Table 2). In terrestrial taxa, selection for this plasticity should be strongest where desiccation risk is greatest, placing tighter constraints on C-value evolution in drier climates. Intriguingly, the contribution of random drift to variation in plethodontid genome size is most evident in the Neotropics (Itgen et al., 2019), an environment defined by selective release from desiccation risk (Baken et al., 2020). Thus, selection on genome size is probably stronger in temperate species more prone to evaporative water loss.

Our results collectively support a new and specific potential context for C-value function. That is, larger genomes appear to increase resistance to water loss in temperate terrestrial lungless salamanders exposed to increased temperature. While this physical effect of C-value is clear in our data, its adaptive functionality is less so. Future tests of this putative function will need to quantify both r and the extent and depth of capillarization and lipid deposits in the skin, compared between species with adequate variation in genome size across a greater thermal breadth to better characterize reaction norms for effective skin thickness and r. The thermal dependence of r in our dataset is modest (per species mean R₁₀=1.41), and our path analyses found no distinction between effects of C-value across temperatures, although only marginally so (P=0.064–0.070). Similarly, we did find an effect of C-value on R₁₀ (Fig. 4B) under our full dataset OU model (P=0.042; Table S9 in Dryad: dryad.t76hdr7zh), but were unable to distinguish this model from Brownian motion. Alongside our occasional difficulty in selecting top models using AIC, these consistently subtle results suggest that some of our analyses may be missing a clearer pattern of temperature-dependent plasticity due to lack of power, either through small sample size, species sampling, measurement error, or experimental design. By only sampling two temperatures, we cannot conclude whether r increases linearly or by some other function with increasing temperature, such as a sudden increase in resistance above a certain temperature threshold. Tests across a wider range of temperatures will increase power to identify temperature effects, and will ultimately reveal the overall thermal reaction norm of resistance to water loss, helping us identify conditions under which physiological performance is most adaptive. Measuring potential covariates of C-value may increase power further and would clarify the mechanism through which genome size influences performance. These issues underscore the need to further test for mechanisms of C-value influencing water loss and desiccation risk in the Plethodontidae. Indeed, C-value may affect analogous physiological performance phenotypes by constraining transport rates across exchange surfaces in other taxa, potentially revealing a previously unappreciated functional consequence of genome size variation. By identifying this research avenue, this study presents a template for integrating data on morphological, physiological and environmental variation underlying an important performance phenotype that will help quantify the role of functional diversity underlying adaptive plasticity in nature.

We thank Reed Ojala-Barbour, Marc Hayes and Jarrett Johnson for assistance in the field, and Kim Sparks for technical assistance with experiments. We also thank Lynn Johnson and Stephen Parry of the Cornell University Statistical Consulting Unit for invaluable assistance with linear and path models, respectively. Thanks also to Eric Riddell and two anonymous reviewers, whose comments greatly improved earlier drafts of the manuscript.