|Home | About | Journals | Submit | Contact Us | Français|
Anuran (frog) metamorphosis has long-served as a model of how thyroid hormones regulate post-embryonic development in vertebrates. However, comparatively little is known about urodele (salamander) metamorphosis. We conducted a detailed time-course study of induced metamorphosis in the Mexican axolotl (Ambystoma mexicanum) that probed metamorphic changes in morphology and gene expression in the skin. Using morphometrics, quantitative PCR, histology, and in situ hybridization we demonstrate that the development of transcriptional markers is fundamental to the resolution of early metamorphic events in axolotls. We then use linear and piecewise linear models to identify a sequence of morphological and transcriptional changes that define larval to adult remodeling events throughout metamorphosis. In addition, we show that transcriptional biomarkers are expressed in specific larval and adult cell populations of the skin and that temporal changes in these biomarkers correlate with tissue remodeling. We compare our results with other studies of natural and induced metamorphosis in urodeles and highlight what appear to be conserved features between urodele and anuran metamorphosis.
Amphibian metamorphosis is a well-studied example of a complex developmental process that is regulated by endocrine factors (Wilder, 1925; Nieuwkoop and Faber, 1967; Rosenkilde and Ussing, 1996; Shi, 2000; Brown and Cai, 2007). Thyroid hormone (TH) control of metamorphosis is broadly conserved across amphibians, and radioimmunoassay (RIA) data from anurans (frogs; e.g., Leloup and Buscaglia, 1977) and urodeles (salamanders; e.g., Larras-Regard et al., 1981; Alberch et al., 1986) support the idea that L-thyroxine (T4; relatively inactive form of TH) and 3,5,3′-triiodothyronine (T3; relatively active form of TH) markedly increase at metamorphic climax. TH (T3 or T4) is necessary and sufficient to induce metamorphosis in anurans (reviewed by Shi, 2000; Brown and Cai, 2007) and urodeles (Prahlad and DeLanney, 1965; Prahlad 1968; Rose, 1995b,c; Rosenkilde and Ussing, 1996) and its biological effects are mediated by nuclear receptors (thyroid hormone receptors α and β; TR-α and TR-β) that repress or activate transcription in a TH-dependent manner (Safi et al., 2004; Buchholz et al., 2006). However, while these and other general features of amphibian metamorphosis are broadly conserved (see Denver et al., 2002), there is considerable variation in the timing, duration, and remodeling patterns that occur across different taxa (Duellman and Trueb, 1994). For example, anurans completely resorb their tails during metamorphosis (Shi, 2000) while urodeles remodel and retain their tails as adults (Duellman and Trueb, 1994). As another example, hind limb formation and growth are intimately linked to metamorphosis in anurans (Shi, 2000), but occur months before metamorphic climax in urodeles (Rosenkilde and Ussing, 1996). Finally, some urodeles (but no anurans) are paedomorphic and altogether fail to undergo a conspicuous metamorphosis thus retaining larval morphological traits and completing their life cycles in water as aquatic adults (see Petranka, 1998).
While the endocrinology (e.g., Leloup and Buscaglia, 1977; Buchholz and Hayes, 2005; Larras-Regard et al., 1981; Alberch et al., 1986), histology (e.g., Heady and Kollros, 1964; Nieuwkoop and Faber, 1967; Fahrmann, 1971a,b,c; Alberch et al., 1985; Ohmura and Wakahara, 1998), and morphology (e.g., Taylor and Kollros, 1946; Nieuwkoop and Faber, 1967; Norman, 1985; Cano-Martinez et al., 1994; Rose, 1995a,b,c) of metamorphosis have been examined in anurans and urodeles, the vast majority of our knowledge about metamorphic gene expression comes from studies of anurans and Xenopus in particular (reviewed by Shi, 2000; Buchholz et al., 2006). This has enabled the conceptualization of models of Xenopus, and, to a lesser extent, Rana metamorphosis that integrate changes in morphology, histology, and gene expression (for examples see Ishizuya-Oka and Shi, 2007; Yoshizato, 2007). However, there have been few attempts to incorporate gene expression data into integrative models of salamander metamorphosis.
Ambystomatid salamanders are a logical urodele family to develop as metamorphic research models because they exhibit interesting life-history variation (Wilbur and Collins, 1973; Shaffer and Voss, 1996; Denver et al., 2002), have a maturing genomic resource infrastructure (Putta et al., 2004; Smith et al., 2005a,b; Monaghan et al., 2009), and are amenable to laboratory culture and conditions (Armstrong et al., 1989; Armstrong and Duhon, 1989; Frost et al., 1989). In particular, the Mexican axolotl (Ambystoma mexicanum) has a long history as a developmental model (Smith, 1989). The axolotl is a paedomorphic salamander that is part of a large species complex (i.e., the tiger salamander species complex) that consists of metamorphic, facultatively paedomorphic, and paedomorphic taxa that have recently diversified from a metamorphic ancestor (Shaffer and Mcknight, 1996; Shaffer and Voss, 1996; Weisrock et al., 2006). Despite rarely spontaneously metamorphosing in nature or the lab, the axolotl can be induced to undergo metamorphosis by adding T3 or T4 to rearing water. We typically induce metamorphosis in the axolotl with T4 because it is generally thought to be the primary iodine containing hormone released by the thyroid and delivered to other tissues where it is locally converted to T3 via deiodinase activity (Denver et al., 2002; Brown, 2005). Fifty nM T4 is an appropriate standard for studying induced metamorphosis in the axolotl because it is well below the threshold of toxicity (80 nM; Rosenkilde and Ussing, 1996), induces metamorphosis at a rate that is not accelerated by higher doses (Rosenkilde and Ussing, 1996), and maximizes rate without altering the sequence of morphological and transcriptional events observed at physiologically relevant concentrations (Page et al., 2008).
During TH-induced metamorphosis, axolotls undergo a series of morphological changes that include reductions in body mass and growth rate, and complete resorption of the tail fins, dorsal ridge, and gills. Comparisons of staging series based on these gross morphological events suggest that the pattern exhibited by TH-induced axolotls (Cano-Martinez, 1994) qualitatively resembles the pattern exhibited by naturally metamorphosing tiger salamanders (Ambystoma tigrinum; Norman, 1985; also see Rosenkilde and Ussing, 1996). Thus the axolotl represents a convenient and informative model system for examining metamorphic gene expression in the tiger salamander complex. In a recent study, we used microarray technology to correlate gene expression in the skin of TH-induced axolotls with metamorphic stage and T4 concentration (Page et al., 2008). These gene expression screens revealed that hundreds of genes are differentially regulated in response to T4 in a way that strongly correlates with morphological stage, but shows little or no correlation with the dose of T4 (5 or 50 nM) used to induce metamorphosis. While potential biomarkers of skin metamorphosis were identified, sampling across developmental time and individuals was too sparse to detail temporal changes in gene expression and morphology. Here, we report on a detailed time-course study of TH-induced metamorphosis of axolotl skin using a group of candidate biomarker genes, quantitative reverse transcriptase real-time PCR (Q-RT-PCR), histology, and in situ hybridization (ISH). We use linear and piecewise linear modeling to identify the onset, offset, and rate of change for morphological and molecular traits during metamorphosis and articulate an integrative model of TH-induced metamorphosis in A. mexicanum. We compare and contrast our results with results obtained from other studies of natural and induced metamorphosis in urodeles and highlight what appear to be conserved features of urodele and anuran metamorphosis.
Sixty Mexican axolotls were obtained from a single inbred cross and reared individually at 20–22 C in 40% Holtfreter’s solution. After hatching, larvae were fed freshly hatched brine shrimp (Artemia) napuli until they were three weeks of age. Beginning at three weeks of age, salamanders were fed California blackworms (Lumbriculus). At 120 days post-fertilization (Day 0), 50 nM T4 (Sigma, St. Louis, MO; T2376) was administered to the animals’ rearing water as described in Page et al. (2007). A ≈ 1 cm2 section of skin was removed from a dorsal region of the head (posterior to the eyes and anterior to the neck) of three individuals beginning on Day 0 (prior to initiating T4 treatment). Skin was then sampled from three individuals every other day for 30 days (i.e., Day 0, Day 2, Day 4, … Day 30). Sampling was also conducted on Day 32; however, only two individuals were sampled at this time point. On Days 0, 12, and 28, additional animals (4, 3, and 3 respectively) were sampled to obtain skin for histological staining and ISH. The remaining 50 tissue samples were used for Q-RT-PCR. All 60 animals were sacrificed at the time of tissue collection. On Day -1, all 60 animals were measured for mass and the seven other morphometric traits described by Fig. 1 (here out referred to as the pre dataset). Just prior to tissue collection, all 60 animals were measured for these same morphometric traits a second time (here out referred to as the post dataset). Morphological staging was conducted according to Cano-Martinez et al. (1994).
After T4 is administered to an axolotl’s rearing water, there is a latency period before metamorphic changes are observed (Rosenkilde and Ussing, 1996; Page et al., 2008). During this latency period, growth continues at a normal rate and metamorphic changes, such as weight loss and tissue resorption, are not evident. At the end of this latency period, there is an onset and changes in one or more traits become obvious as morphological metamorphosis proceeds. Similarly, because metamorphic change ultimately results in new adult morphology, many metamorphic traits also undergo an offset after which they no longer change because they have taken on their adult form. These onset and offset times can be conceptualized as transitions or breakpoints that connect linear equations with different slopes (rates of change). These slopes describe the direction and rate of change for a trait before the onset/offset and the direction and rate of change after the onset/offset (Fig. 2). Piecewise linear regression is a useful tool because it can enable one to estimate: (1) the time at which a transition (onset/offset) occurs, (2) the direction and rate of change before the transition, and (3) the direction and rate of change after the transition. In its simplest form, a piecewise linear model can be defined as follows: for xi ≤ Ψ, yi = β0 + β1xi + εi and for xi > Ψ, yi = β0 + β1xi + β2(xi − Ψ) + εi where yi is the value observed for a given trait from individual i, β0 is the baseline value for the trait, β1 is the slope or rate of change to the left of the transition, xi is the time at which individual i was sampled, β2 is the slope to the right of the breakpoint minus the slope to the left of the breakpoint (Fig. 2), Ψ is the transition time or breakpoint, and εi is the error associated with individual i (for a detailed discussion of piecewise regression see Toms and Lesperance, 2003). Thus, it follows that the rate of change (slope; m) to the right of the breakpoint can be obtained as follows: m = β1 + β2.
We performed piecewise regression as implemented by the “segmented” software package (Muggeo, 2008) that is freely available for the R statistical computing environment (www.r-project.org). All models in which morphological traits (t) were treated as response variables were conducted using the differences between the pre and post measurements. The following formula: Δt = tpre − tpost makes resorption (for most morphometric traits the phenomenon of interest) positive. Thus, positive values of Δt indicate that a trait is vanishing and negative values of Δt indicate that a trait is growing. For morphological traits where piecewise relationships were not evident, we fit linear regression models of the form yi = β0 + β1xi + εi where β0 is the baseline trait value, β1 is the slope, and ε is the error term (Zar, 1999). Such linear models provide estimates of the direction and rate of change for traits (slopes) that do not exhibit evidence of metamorphic onset or offset within the experimental interval. In all cases, intercept free models (β0 = 0) were fit so as to force lines through the origin as, by definition, all traits have initial values of zero days T4 administration (x-axis) and zero change between the pre and post datasets (y-axis).
We previously used microarray technology to identify hundreds of genes that are differentially expressed in response to T4 treatment (Page et al., 2007; Page et al., 2008). This gene list included keratin 6A (KRT6A), keratin 14 (KRT14), uromodulin (UMOD), calmodulin 2 (CALM2), epithelial membrane protein 1 (EMP1), glutamate ammonia-ligase (GLUL), and matrix metalloproteinase 1 (MMP1). To the best of our knowledge, only UMOD and MMP1 are known TH response genes from other amphibian taxa. UMOD-like genes have been identified from Xenopus laevis and are thought to contribute to a larval specific structure on the surface of the skin that has “protective” and “secretory” functions (Furlow et al., 1997). MMP1 is a gene that is known to be directly targeted by TH in Rana catesbeiana (Oofusa and Yoshizato, 1996; Sawada et al., 2001). In R. catesbeiana, MMP1 is involved in collagen turnover that occurs during tail resorption and skin remodeling (Oofusa and Yoshizato, 1991). The functions of the rest of these genes (KRT6A, KRT14, GLUL, EMP1, CALM2) during amphibian metamorphosis are largely unknown. We assign gene names to our assembled and curated expressed sequence tag (EST) based contigs (Putta et al., 2004; Monaghan et al., 2009) in accordance with their presumptive orthologies to human based on BLASTx searches of the human RefSeq protein database (E < 10−7; see Monaghan et al., 2007; Page et al., 2007). Thus, sequence homology leads us to assume that KRT6A and KRT14 encode cytoskeletal proteins that are constituents of intermediate filaments. Similarly, we assume that CALM2 encodes a protein involved in several aspects of calcium regulation and signaling (Friedberg and Rhoads, 2001). We also assume GLUL encodes an enzyme that plays important roles in nitrogen metabolism (Kumada et al., 1993). Finally, we assume EMP1 is associated with cell proliferation (Ben-Porath and Benvenisty, 1996) and epithelial differentiation (Marvin et al., 1995). In addition to genes that were identified by previous expression screens, we have chosen to investigate TR-α and TR-β due to the fundamental roles that these transcription factors play in TH signaling during metamorphosis (see Safi et al., 2004; Buchholz et al., 2006). Table 1 summarizes the genes that we have investigated as well as their presumptive functions.
Total RNA was extracted from 50 tissue samples using TRIzol (Invitrogen, Carlsbad, CA) according to the manufacturer’s protocol. RNA samples were further purified using Qiagen RNeasy mini columns. All samples were treated with RNase-Free DNase Sets (Qiagen, Valencia, CA) according to the manufacturer’s protocol. UV spectrophotometry and gel electrophoresis were used to quantify and qualify RNA preparations.
Relative transcript abundances for the genes listed in Table 1 were estimated via Q-RT-PCR according to the protocol described in Page et al. (2008). Briefly, primers (Table 1) were designed using Primer3 (Rozen and Skaletsky, 2000). All PCRs were 25 μL reactions that contained cDNA corresponding to 10 ng total RNA, 41 ng forward and reverse primers, and iQ SYBR-Green real-time PCR mix (Bio-Rad, Hercules, CA). Reaction parameters were 10 minutes at 50 C, five minutes at 95 C, 45 cycles of 10 seconds at 95 C followed by 30 seconds at 55 C, one minute at 95 C, and one minute at 55 C. For all reactions, melting curve analysis was performed to ensure that only a single product was amplified. PCRs were performed on a Bio-Rad iCycler iQ Multi-Color Real Time PCR Detection System. All plates contained template free controls (Bustin and Nolan, 2004). Primer efficiencies were estimated via linear regression and relative expression ratios (R) were calculated according to Pfaffl (2001). All expression ratios are relative to the average of Day 0 animals, and normalized to the expression of transcriptional intermediary factor 1 (Table 1; Page et al., 2008).
Just as many of the morphological changes that occur during induced metamorphosis can be conceptualized as offset times, so can many of the gene expression changes that occur during metamorphic skin remodeling. For example, larval specific and adult specific genes are often down and up-regulated respectively during metamorphosis to new adult expression levels. Thus, as is the case with morphological traits, there are transitions or breakpoints that connect the rate of change observed during the metamorphic period with the rate of change observed after metamorphic gene regulation has ended (i.e., metamorphic offset). Hence, the simple piecewise model described in Fig. 2 allows one to estimate several useful parameters including: β1 the rate of change in transcript abundance during the metamorphic period, Ψ the offset time of metamorphic gene regulation, and m the rate of change in transcript abundance after the offset of metamorphic gene regulation. This same piecewise model (Fig. 2) also enables the estimation of meaningful parameters for genes associated with metamorphic remodeling processes. Such remodeling genes are frequently transiently up or down-regulated before returning to baseline expression levels (Page et al., 2008). Thus, β1 provides an estimate of the rate of change in transcript abundance prior to reaching maximum/minimum expression levels, Ψ provides an estimate of the time at which maximum/minimum expression levels are reached, and m provides an estimate of the rate of return to baseline expression levels. In keeping with these ideas, we applied piecewise and linear regression, as described above for morphological traits, to our gene expression data. In order to decide between piecewise and linear models for each gene, we used Akaike’s information criterion (AIC; Venables and Ripley, 2002) which attaches a numerical value (smaller AIC is better) to each model that accounts for the tradeoff between model fit and model complexity (i.e., the number of parameters in a model). All models of gene expression treat days of T4 treatment as the predictor variable and log2 transformed R-values as the response variable. As with the morphology data, we used models that lack intercepts (β0 = 0) to estimate β1, Ψ, and m (m is estimated indirectly via knowledge of β1 and β2; see Fig. 2) as larval expression levels on a log2 scale (y-axis) are by definition zero and sampling began following zero days of T4 exposure (x-axis). By assuming that baseline expression for each gene is equal to the mean of the three individuals observed at Day 0 (i.e. that β0 = 0), we are able to obtain consistent estimates of β1, Ψ, and m that put the expression of all genes on a common scale. However, because our estimation of R for each gene at Day 0 rests on a small sample size, we use models that estimate β0 (baseline expression) from all 50 data points via least squares to conduct hypothesis testing (i.e., to make statements about differential regulation in a statistical sense). When we were unable to estimate piecewise models due to convergence problems in the model fitting routine, we used quadratic regression to test for the presence of curvilinear relationships that are similar in shape to the piecewise model described in Fig. 2. Quadratic models take on the following form: log2(R)i = β0 + γ1xi + γ2x2i + εi where β0 and ε are the baseline expression value and error term respectively and γ1 and γ2 are the linear and quadratic regression coefficients respectively. Quadratic models provide a well-established framework for testing for the presence of curvilinear relationships (Zar, 1999), but do not provide parameter estimates with interpretations as straight forward as the model described in Fig. 2.
In order to show that TH induces tissue remodeling events in the skin, we conducted hematoxylin and eosin staining at three time points: Days 0, 12, and 28. Skin samples were fixed at 4 C in 1x phosphate buffer saline (PBS), 4% paraformaldehyde (PFA) overnight. Following cryoprotection in 30% sucrose, tissues were sectioned at 20 μm using a Microm 500 HM cryostat (Richard Allen Scientific, Kalamazoo, MI) and stained with Gill’s hematoxylin #2 (Sigma, St. Louis, MO) and eosin Y (Amresco, Solon, OH). Microscopy was performed using an Olympus AX80 microscope and images were acquired with an Olympus DP70 camera (Olympus, Center Valley, CA).
We used ISH to verify Q-RT-PCR results and spatially localize gene expression at the cellular level. Tissues were collected and prepared as described above for histological staining. Day 0 and Day 28 sections were collected on the same SuperFrost Plus microscope slides (VWR, West Chester, PA) and analyzed simultaneously. ISHs were performed as described in Monaghan et al. (2007), with minor modifications. Tissues were not decalcified and sectioning was conducted at 20 μm rather than 16 μm. Images were captured using an AX-80 Olympus microscope (Olympus, Center Valley, CA). Hybridizations and washes were performed overnight at 63.5 C, except for CALM2 for which hybridizations and washes were performed at 58 C. Probes for each gene are shown in Table 1.
We observed morphological and histological changes that are consistent with TH-induced metamorphosis. The first Stage 1 (Cano-Martinez et al., 1994) animal was observed at Day 8 and individuals completed metamorphosis (Stage 4) between Day 28 and Day 32 (Fig. 3). We note that the time required for a given animal to reach a particular Cano-Martinez stage was quite variable (Fig. 3), suggesting that even among full siblings raised individually under identical conditions, there is variation in response to T4. Histological changes were consistent with Fahrmann’s (1971a,b,c) descriptions of metamorphic skin remodeling (Fig. 4). Briefly, larval cell types (Fig. 4A–C) such as apical cells and Leydig cells were removed while adult cell types (Fig. 4G–I) increased in abundance. By Day 12 (Fig. 4D–F), the skin was chimeric and consisted of larval and adult cell types. By Day 28, the epidermis consisted of a heavily keratinized stratified squamous epithelium similar to that observed in metamorphosed anurans and adult mammals (Fig. 4G–I).
Of the eight morphological traits investigated (Fig. 1), five (LTH, GL, mass, TL, and SVL) exhibited piecewise relationships and two (UTH and TH) exhibited linear relationships (Fig. 5A–G). As measured, DRL is not a continuous variable (Fig. 1) and therefore could not be modeled in a regression framework. However, graphical inspection indicates that the dorsal ridge is completely resorbed between Day 18 and Day 32 (Fig. 5H). Inspection of the slope of the linear regression on ΔUTH (Fig. 5A; β1 = 0.022 cm/day T4 treatment) and the left line segment of the piecewise model fit to ΔLTH (Fig. 5B; β1 = 0.012 cm/day of T4 treatment) suggest similar resorption rates for the upper and lower tail fins. Thus, it is not surprising that the rate of reduction in total tail height was also of a similar magnitude (Fig. 5C; β1 = 0.030 cm/day of T4 treatment). The offset time corresponding to complete lower tail fin resorption was estimated to occur following 19.73 days of T4 treatment; however, there is considerable uncertainty associated with this estimate (lower 95% confidence interval (LCI) = 14.69 days, upper 95% confidence interval (UCI) = 24.76 days). As expected, following complete resorption, the slope of the right line segment for ΔLTH was not statistically different from zero (m = −0.002, LCI = −0.009, UCI = 0.006).
The pattern observed for gill resorption (Fig. 5D) was markedly different than the patterns observed for upper and lower tail fin resorption. The slope of the left line segment of the piecewise model fit to ΔGL was not statistically different from zero (β1 = −0.004, t = −0.596, P = 0.554) suggesting that gill resorption does not begin immediately. The onset of gill resorption was estimated to occur following 18.82 days of T4 treatment (LCI = 15.26 days, upper UCI = 22.39 days). Upon commencement, gill resorption was rapid and occurred at a rate that was nearly an order of magnitude faster than lower tail fin resorption (m = 0.110 cm/day of T4 treatment).
Analyses of the traits associated with body size (mass, TL, SVL) showed general concordance. For example, piecewise regression models for all three traits had negative β1 estimates indicating that animals continued to grow in length and mass during the first half of the induction interval (Fig. 5E–G). Similarly, the breakpoints for Δmass (Fig. 5E; Ψ = 17.14 days of T4 treatment, LCI = 11.87 days, UCI = 22.42 days), ΔTL (Fig. 5F; Ψ = 12.29 days of T4 treatment, LCI = 8.23 days, UCI = 16.34 days), and ΔSVL (Fig. 5G; Ψ = 16.24 days of T4 treatment, LCI = 11.29 days, UCI = 21.20 days) suggest that a reduction in growth rate occurred following ≈ 12–17 days of exposure to 50 nM T4. The relationship between Δmass and days of T4 administration deteriorated following metamorphic onset (Fig. 5E). However, as expected, TL (Fig. 5F) and SVL (Fig. 5G) showed similar patterns throughout metamorphosis. Before metamorphic onset the slope (β1) for ΔTL was estimated to be −0.124 cm/day of T4 treatment whereas after metamorphic onset, the slope for ΔTL (m) was reduced to −0.035 cm/day of T4 treatment. For ΔSVL β1 = −0.065 cm/day of T4 treatment and m = −0.016 cm/day of T4 treatment. These results suggest that, after a relatively long latency period, growth rates are reduced by 70–75% during T4 induced metamorphosis of A. mexicanum. Fig. 5I summarizes the equations that we estimated for each morphometric trait using linear and piecewise regression.
We recovered four basic temporal gene expression profiles. First, we observed piecewise patterns in which m ≈ zero (e.g., UMOD). Second, we observed piecewise patterns in which m < β1 but m ≠ zero (e.g., CALM2). Third, we observed piecewise patterns that are consistent with the notion of a transient response to T4 that is followed by a return to baseline expression levels (e.g., MMP1). Fourth, we observed linear expression profiles (e.g., TR-α). The shapes (piecewise/quadratic versus linear) of the intercept free models (β0 = 0) that were selected for each gene and their parameter estimates are shown in Table 2. Attempts to fit intercept free piecewise models to TR-α and GLUL were thwarted by convergence issues that are due to the difficulties associated with estimating Ψ for these genes (i.e., several possible breakpoints are essentially equally satisfactory). However, the intercept free quadratic models (β0 = 0) fit to TR-α and GLUL both had larger AIC values than the linear models fit to these genes that were selected. These results support the notion that linear models are sufficient to describe TR-α and GLUL expression.
Three down-regulated genes (UMOD, KRT6A, and CALM2) exhibited piecewise relationships. UMOD was drastically down-regulated in response to T4 (Fig. 6A). Taking the absolute value of the log2(R) predicted for UMOD at time Ψ (≈ 12.78) and back transforming to the raw scale (212.78) results in an estimate of ≈ 7000 fold down-regulation in UMOD transcripts prior to metamorphic offset (Ψ = 14.340 days). This suggests UMOD transcripts are down-regulated at an average rate of approximately (7000/14 = 500) 500 fold per day of T4 treatment prior to metamorphic offset. Analogous calculations for KRT6A (Fig. 6B) suggest that KRT6A is down-regulated by ≈ 65 fold by its time of metamorphic offset (Ψ = 11.540 days). This gives an average rate of a (65/11.540 = 5.6) ≈ 5.6 fold decrease in KRT6A transcript abundance per day of T4 treatment prior to metamorphic offset. The piecewise model fit to CALM2 (Fig. 6C) requires a different interpretation than the models fit to UMOD and KRT6A. The breakpoint in the CALM2 model was estimated to occur following 4.68 days of T4 treatment and corresponds to a transition in the rate of CALM2 down-regulation. The estimate of β1 for CALM2 is over an order of magnitude larger than the estimate of m for CALM2 (Table 2). By time Ψ, CALM2 was down-regulated by approximately 5 fold and by the end of the experimental period (Day32) CALM2 was down-regulated by approximately 15 fold.
We observed three up-regulated genes (KRT14, MMP1, and EMP1) that exhibited piecewise profiles. KRT14, was strongly up-regulated in response to TH (Fig. 6D). Calculations analogous to those computed previously suggest that KRT14 is up-regulated by ≈ 250 fold by time Ψ resulting in an average rate of a (250/10.380 = 24.1) ≈ 24 fold increase in transcript abundance per day of T4 treatment prior to Ψ. By the end of morphological metamorphosis (Day 32) KRT14 transcripts were up-regulated by approximately 650 fold relative to Day 0 expression levels. MMP1 expression was variable throughout metamorphosis and our temporal expression data for MMP1 are equivocal (Fig. 6E). However, model selection (AIC) resulted in a model whose general interpretation (transient expression) is consistent with the presumptive role of MMP1 during amphibian metamorphosis (collagen remodeling). The estimate for Ψ suggests that MMP1 expression levels were maximal following roughly 21 days of T4 treatment. However, there is considerable uncertainty associated with the estimate of Ψ for MMP1 (Table 2). EMP1 was up-regulated (Fig. 6F) by approximately 17 fold by time Ψ resulting in an average rate of a (17/16.330 ≈ 1.0) ≈ one fold increase in EMP1 transcript abundance per day of T4 treatment prior to Ψ. However, in the case of EMP1, the interpretation of Ψ is unclear. Because the estimate of m is negative for EMP1 (Table 2), it is not clear whether EMP1 levels remain elevated outside of our experimental interval (i.e., in animals that have completed morphological metamorphosis).
Finally, we observed three genes (TR-β, TR-β, and GLUL) with linear expression profiles. The linear model selected for TR-β (Fig. 6G) suggests that TR-α is down-regulated by approximately 3.4 fold by the end of morphological metamorphosis (Day 32). This implies that TR-α transcripts decrease in abundance at an average rate of roughly (3.4/32 = 0.11) 0.11 fold per day of T4 treatment. Similarly, the linear model selected for TR-α (Fig. 6H) suggests that TR-α is down-regulated by approximately 5.1 fold by Day 32. Thus TR-α transcripts decreased in abundance at an average rate of approximately (5.1/32 = 0.16) 0.16 fold per day of T4 treatment throughout the experimental period. GLUL expression was extremely variable and its magnitude of change was modest (Fig 6I; back transformation of the predicted value at Day 32 suggests roughly a two fold decrease relative to Day 0). Overall, our Q-RT-PCR data highlight the dramatic changes in gene expression that occur in the skin during TH-induced metamorphosis of A. mexicanum. Fig. 7 summarizes the intercept free equations that we estimated via linear and piecewise regression and shows the timing of these gene expression patterns relative to key morphological events.
Estimating β0 via least squares generally supported the notion that assuming β0 = 0 (i.e., fitting intercept free models) is valid (Fig 6; Table 3). Comparing AIC values from linear models for which β0 was estimated with AIC values from piecewise or quadratic models for which β0 was estimated did not result in the selection of models with different shapes (linear versus piecewise/quadratic) from those presented in Table 2 for any of the genes investigated. The general correspondence between intercept free models and models in which β0 was estimated from the data suggests that our results are robust (Fig 6). We note that seven of the models (UMOD, KRT6A, CALM2, KRT14, EMP1, TRβ, and TR-α) used to estimate β0 account for > 25% of the variation in the data (i.e., R2 > 0.25) with six of these accounting for > 33%, four accounting for > 50%, and two accounting for > 90% (Table 3; see below for KRT6A).
Inspection of Table 3 reveals that the estimates of β1 for UMOD, CALM2, KRT14, EMP1, TR-β, and TR-α are clearly statistically different from zero. Attempts to fit a piecewise model to KRT6A that estimated β0 failed to converge. However, the quadratic model fit to KRT6A (Fig. 6B) had a lower AIC than the analogous linear model and was highly statistically significant (n = 50, df = 47, β0 = −1.11, t = −2.38, P = 0.022, γ1 = −0.485, t = −7.024, P = 7.52 × 10−9, γ2 = 0.012, t = 5.854, P = 4.47 × 10−7, R2 = 0.567). Although the estimate of β1 for GLUL is statistically different from zero at the 0.05 level (Table 3), expression was too variable to allow for inferences on GLUL’s role during metamorphosis, and we do not attach biological significance to this result. The estimate of β1 for MMP1 was not statistically different from zero (Table 3).
Examination of the estimates for m and their 95% confidence intervals reveal that m is not different from zero at the 0.05 level for UMOD and EMP1 (Table 3). Our interpretation of this with respect to EMP1 was discussed in the previous section. However, for UMOD this result suggests that metamorph expression levels are reached at time Ψ (Fig. 6A). Because the estimates of m for CALM2 and KRT14 statistically differed from zero, it seems likely that expression levels for these genes decrease and increase respectively in a more gradual fashion after their respective transition times (Fig. 6C; Fig. 6D). MMP1 represents an interesting case in the sense that piecewise models consistent with transient regulation were selected in favor of simpler linear models irrespective of whether β0 was estimated (see above). When this result is considered in combination with the observation that a value unique to piecewise models (m) is statistically significant but a parameter with a clear analog in simpler linear models (β1) is null, the idea that MMP1 is transiently up-regulated in response to TH is supported. While acknowledging that MMP1 expression is clearly variable, we note that the ISH data presented below are consistent with the view that MMP1 is up-regulated in response to T4 (see below) as are our results from previous expression screens (Page et al., 2008).
In order to verify the temporal patterns described above and associate candidate biomarkers with specific cell populations, we conducted ISH. The spatial expression data obtained via ISH confirm that UMOD (Fig 8A–C), KRT6A (Fig. 8D–F), and CALM2 (Fig. 8G–I) are larval specific. UMOD and KRT6A were expressed exclusively in the apical cell layer that is present only in the larval epidermis where as CALM2 was expressed by basal Leydig cells and dermal fibroblasts (Fig. 8G–I). MMP1 was not detectable in larval skin (Fig. 8J). However, at Day 12 and Day 28, MMP1 was expressed in dermal glands (Fig. 8K, Fig. 8M). ISH confirmed that KRT14 is metamorph specific (Fig. 8N–P) and revealed that KRT14 is strongly and exclusively expressed in the granular layer of the metamorph epidermis (Fig. 8O).
To facilitate further study of metamorphosis in the axolotl, we conducted a detailed time-course study in which we used 50 nM T4 to induce metamorphosis in Mexican axolotls while tracking changes in morphological and transcriptional events. As Fig. 7 makes clear, the earliest events occur at the transcriptional level. In particular, CALM2, KRT6A, KRT14, and UMOD are all substantially differentially regulated between the onset of T4 administration and approximately 14 days of T4 administration. Around the same time that the last of these genes (UMOD) undergoes metamorphic offset, there is a reduction in growth rate, as estimated from SVL and mass, that is accompanied by maximal EMP1 expression levels. This reduction in growth rate is followed by advanced lower tail fin and dorsal ridge resorption, the former of which is completely resorbed following roughly 20 days of T4 treatment. At approximately the same time that the lower tail fin is completely resorbed, MMP1 expression levels are maximal and the gills begin to undergo rapid resorption. Between Day 28 and Day 32, the gills and upper tail fin are completely resorbed, resulting in what has traditionally been considered the end of metamorphosis (Prahlad and DeLanney, 1965; Voss and Smith, 2005). Collectively, these results demonstrate the necessity of transcriptional biomarkers for indexing early metamorphic events and highlight the dynamic changes that occur during induced metamorphosis of the axolotl. Below, we summarize how our gene expression data will enable future studies to index metamorphic progress in the axolotl prior to morphological changes. We then compare and contrast natural and induced metamorphosis in urodeles and discuss the remodeling patterns observed in the epidermis and dermis. Finally, we briefly highlight what appear to be conserved features of urodele and anuran metamorphosis.
Of the genes that we investigated, UMOD, KRT6A, and CALM2 are larval specific. Thus, the responses of these genes to TH are useful for indexing loss of the larval gene expression program in the skin. UMOD and KRT6A are exclusively expressed in the apical cell layer, a larval specific structure that is completely removed during metamorphosis. Hence, UMOD and KRT6A specifically mark the loss of apical cells from the skin during metamorphosis. In addition to these larval specific genes, we have also demonstrated that KRT14 is metamorph specific and expressed exclusively in the granular layer of the metamorph epidermis. For these reasons, KRT14 has considerable utility as a biomarker of differentiation and proliferation of the metamorph epidermis. Of the remaining genes investigated (TR-α, TR-β, EMP1, MMP1, and GLUL), the TRs and EMP1 showed the most potential as transcriptional indices of metamorphosis. Because TH sensitivity is dependent on the presence of TRs (Yoshizato and Frieden, 1975), TR-α and TR-β are useful biomarkers of TH sensitivity. Thus the linear down-regulation that we observed for TR-α and TR-β probably reflects a decrease in TH sensitivity in the skin between metamorphic climax and the end of morphological metamorphosis (see below). With respect to EMP1, in mammals this gene is associated with cell proliferation (Ben-Porath and Benvenisty, 1996) and epithelial differentiation (Marvin et al., 1995). The expression pattern that we recovered for EMP1 is consistent with the axolotl ortholog playing similar roles during TH-induced metamorphosis. Finally, MMP1 and GLUL showed the least potential as metamorphic indices. As mentioned above, MMP1 and GLUL expression were variable; however, MMP1 may have utility as a marker of extra cellular matrix (ECM) remodeling in the dermis (see above; see below). Overall, the genes investigated in this study increase the size of the molecular toolkit for studying TH-induced metamorphosis in the axolotl. These markers identify larval cells (CALM2, KRT6A, UMOD), adult cells (KRT14), factors that influence TH sensitivity (TR-α and TR-β), and putatively, epidermal proliferation and differentiation (EMP1).
In naturally metamorphosing amphibians, serum T3 and T4 gradually rise prior to metamorphic climax (originally proposed by Etkin, 1968 and best documented in anurans; see Shi, 2000 and Denver et al., 2002 for reviews), and the TH competence of cells and tissues increases during the course of larval development (e.g., Rose, 1995b). In contrast, an induced metamorphosis is typically achieved by suddenly exposing larvae to a high dose of one form of TH (T3 or T4); a practice whose effect on serum T3 and T4 levels is not fully understood. Galton (1992) found that juvenile axolotls injected with radiolabeled T4 had undetectable serum T3 at 20 hours post-injection suggesting that T4 is not immediately converted to T3 in large quantities. In agreement with Galton (1992), Alberch et al. (1985) were unable to detect T3 in Eurycea bislineata larvae immersed in 50 nM T4 at 24 hours post-immersion. Moreover, Alberch et al., (1985) observed serum T4 levels as high as 41 nM during the first 48 hours of immersion followed by a leveling off to physiologically relevant T4 concentrations (≈ 1–11 nM; see Alberch et al., 1986). Thus, in addition to discrepancies that may arise between natural and induced metamorphosis due to gradual versus sudden TH exposure, available studies suggest that the serum TH profiles of induced and naturally climaxing urodeles may be quite different, especially early during their respective TH surges.
The dose of T4 that we used to induce metamorphosis is roughly twice the maximum serum T4 level observed in a naturally metamorphosing tiger salamander (28 nM; Larras-Regard, 1981). In a previous study (Page et al., 2008), we found that this relatively high dose of T4 (50 nM) reduces the latency period associated with induced metamorphosis in the axolotl, but does not induce morphological and transcriptional changes that differ from the patterns observed when metamorphosis is induced with a physiologically relevant concentration of T4 (5 nM). Studies that we have conducted on eastern tiger salamanders (Ambystoma tigrinum tigrinum) and tiger salamander-axolotl hybrids show that 120 days post-fertilization is within the range of ages that metamorphic onset occurs spontaneously under laboratory conditions (Voss and Smith, 2005; unpublished data; also see Ducibella, 1974). Thus, the current study presumably reflects a level of TH competence in tissues and cells that is similar to those of naturally metamorphosing members of the tiger salamander complex. While our experimental paradigm is clearly contrived, we have attempted to minimize the disparities between our approach and natural metamorphosis. Despite imperfections, we think that the ability to uniformly administer T4 to animals that are phylogentically close to naturally metamorphosing taxa, but not developing toward metamorphic endpoints, represents a powerful model system that provides information relevant to natural metamorphosis.
The remodeling patterns that we observed show that axolotls rapidly remove the larval epidermis while simultaneously constructing the metamorph epidermis. Thus, while gross morphological changes, such as tissue resorption and weight loss are delayed, gene expression changes for some genes in the epidermis are substantial quite early in the induction period. The fact that UMOD, KRT6A, and KRT14 were drastically down, down, and up-regulated respectively shows that axolotls rapidly down-regulate genes that function as larval barriers while rapidly up-regulating genes with protective function in the metamorph epidermis. The mixture of larval and metamorph cell types that we observed at Day 12 corroborates this account of simultaneous removal of the larval epidermis and construction of the metamorph epidermis. In contrast, the remodeling observed in the dermis was far less dramatic, and a well-formed glandular dermis was present in Day 0 axolotls prior to T4 treatment. Holder and Glade (1984) showed that axolotls naturally develop a well-formed glandular dermis during larval development and that dermal glands form as clusters of basal epidermal cells that pass beneath the basement membrane. Ohmura and Wakahara (1998) observed similar patterns for dermis and gland development, but found that these events were associated with epidermal remodeling in Hynobius retardatus. Interestingly, when Ohmura and Wakahara (1998) inhibited thyroid function in H. retardatus via hyphophysectomy, thyroidectomy, or treatment with sodium perchlorate, they found that dermal thickening and gland formation occurred in the absence of epidermal transformation. Thus, either dermis formation is independent of TH in at least some urodeles, or very low concentrations of TH are required for dermis development. While studying the metamorphic cranial and hyobranchial remodeling that occur in E. bislineata, Rose (1995c) found substantial variation in the TH sensitivities of the tissues that comprise these structures. By immersing Eurycea larvae in a series of T4 concentrations ranging from 0.05 nM to 50 nM, Rose (1995c) was able to show that doses as low as 0.05 nM T4 (well below the sensitivity of RIAs used in studies of urodele metamorphosis) are capable of stimulating metamorphic remodeling in some tissues. Thus, it is possible that, as previously suggested by Ducibella (1974), Rosenkilde and Ussing (1996) and Brown (1997), TH plays important roles in the natural development of axolotls despite the low serum T4 estimates from RIAs (Darras and Kuhn, 1983, 1984; Galton, 1992). Irrespective of whether dermis development is TH dependent in the axolotl, our ISH results for MMP1 suggest that T4 induces remodeling in the ECM of the dermis. In addition, our ISH results suggest that dermal fibroblasts down-regulate CALM2 expression following exposure to T4. Thus, cells within the dermis did undergo some differential transcription and/or mRNA metabolism in response to 50 nM T4 despite modest histological remodeling.
Most of the genes that we have investigated in this study are not known TH response genes in other amphibian taxa. However, TR-α, TR-β, UMOD, and MMP1 expression have been examined in anurans. In X. laevis, TR-α is expressed early in larval development whereas TR-β abundance increases with TH levels (Yaoita and Brown, 1990; Kawahara et al., 1991). Thus, in Xenopus, TR-α is thought to be a repressor of metamorphic gene expression in the absence of TH (i.e., during pre-metamorphosis) whereas TR-β expression is associated with activation of the metamorphic gene expression program (reviewed by Buchholz et al., 2006). In X. laevis, it is known that TR-β is directly targeted by TH (i.e., contains a thyroid hormone response element; Ranjan et al., 1994; Wong and Shi, 1995) and that exogenous TH can induce Xenopus TR-β expression precocially (e.g., Yaoita and Brown, 1990; Kanamori and Brown, 1992). While axolotl TRs have been shown to function in mammalian cells (Safi et al., 2004), very little is known about their expression patterns during natural development. Galton (1992) found that putative T3 receptors were more abundant in juvenile axolotls than in paedomorphic adults suggesting that beyond a certain age TH sensitivity may decrease in the axolotl naturally. Because frogs and salamanders are separated by over 200 million years of independent evolution (Kardong, 1998; Anderson et al., 2008), and the axolotl (Shaffer, 1993) and Xenopus (Bolker, 1995) represent derived modes of development within their respective orders, it is not clear how to interpret our results on TR expression within the broader context of amphibian metamorphosis. While we did not observe any evidence of an increase in the abundance of TR-β, this could be explained by a number of possible scenarios. For example, TR-β may have been up-regulated in response to low concentrations of endogenous TH, or by a TH independent mechanism, during the larval period prior to the beginning of our experiment. If this is the case, TR-β may have already been at maximal levels by 120 days post-fertilization. Irrespective of these considerations, our results from axolotl are consistent with the observation from Xenopus that TRs decrease in abundance subsequent to metamorphic climax (Yaoita and Brown, 1990; Kawahara et al., 1991).
Despite the evolutionary considerations cited above, the expression of UMOD-like genes is strikingly similar between A. mexicanum and X. laevis. In X. laevis, two genes with clear homology to human UMOD are expressed exclusively in apical cells and are dramatically down-regulated in response to TH (Furlow et al., 1997). Furlow et al. (1997) hypothesized that Xenopus UMOD genes are part of a larval-specific complex that has “protective” and “secretory” functions and is homologous to the mammalian periderm. Our current results are extraordinarily similar to those described for X. laevis demonstrating that the expression of UMOD-like genes is strongly conserved between Xenopus and axolotl with respect to spatial location and TH responsiveness.
Similarly, there appear to be parallels between axolotl MMP1 expression and MMP1 expression in Rana catesbeiana (bullfrog). Using immunohistochemistry, Oofusa and Yoshizato (1991) found that bullfrog MMP1 protein increased in back skin in response to T3 treatment. These results are interesting because skin from the tadpole trunk undergoes remodeling that is histologically similar to the remodeling observed in the current study (reviewed by Yoshizato, 2007). Of particular interest is that Oofusa and Yoshizato (1991) reported that bullfrog MMP1 is found in collagen fibers and fibroblasts located in the dermis. Thus, it seems that MMP1 is expressed in the dermis of metamorphosing axolotls and bullfrogs. Collectively, these findings suggest that several aspects of amphibian metamorphosis appear to be conserved between urodeles and anurans at the level of gene expression.
Studying metamorphic remodeling in urodeles at multiple levels of biological organization provides comparative perspective and insights into amphibian metamorphosis. In the axolotl, additional temporal and spatial profiling of TH responsive genes from the skin and other tissues is needed to identify additional biomarkers. Within the skin, transcriptional markers of Leydig cells, basal cells, and various adult cell types are needed so that the fate of different epidermal components at metamorphosis can be better understood. In addition, it will be important to characterize TR expression early during axolotl development, so that the roles of these transcription factors during natural development and induced metamorphosis can be better understood. In the current study, we have identified a sequence of transcriptional and morphological events that will provide a useful benchmark for future studies of induced metamorphosis in the Mexican axolotl. In addition, we have identified biomarkers of larval skin loss and adult skin development. This information will enable future studies to assess metamorphic changes in the skin well before any signs of morphological metamorphosis are evident.
This project was supported by grants R24-RR016344 and P20-RR016741 from the National Center for Research Resources (NCRR), a component of the National Institutes of Health (NIH). Its contents are solely the responsibility of the authors and do not necessarily represent the official views of NCRR or NIH. This project was also supported by the National Science Foundation funded Ambystoma Genetic Stock Center (DBI-0443496). Phil Crowley provided helpful comments on an earlier draft of this manuscript. Arnold Stromberg provided statistical advice. The content of this manuscript was improved by comments from two anonymous reviewers.
Publisher's Disclaimer: This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.