Effect of thyroid hormone concentration on the transcriptional response underlying induced metamorphosis in the Mexican axolotl (Ambystoma)
- Robert B Page†1, 2,
- Stephen R Voss†1, 2Email author,
- Amy K Samuels1, 2,
- Jeramiah J Smith1, 2,
- Srikrishna Putta1, 2 and
- Christopher K Beachy3, 4
© Page et al; licensee BioMed Central Ltd. 2008
Received: 02 November 2007
Accepted: 11 February 2008
Published: 11 February 2008
Thyroid hormones (TH) induce gene expression programs that orchestrate amphibian metamorphosis. In contrast to anurans, many salamanders do not undergo metamorphosis in nature. However, they can be induced to undergo metamorphosis via exposure to thyroxine (T4). We induced metamorphosis in juvenile Mexican axolotls (Ambystoma mexicanum) using 5 and 50 nM T4, collected epidermal tissue from the head at four time points (Days 0, 2, 12, 28), and used microarray analysis to quantify mRNA abundances.
Individuals reared in the higher T4 concentration initiated morphological and transcriptional changes earlier and completed metamorphosis by Day 28. In contrast, initiation of metamorphosis was delayed in the lower T4 concentration and none of the individuals completed metamorphosis by Day 28. We identified 402 genes that were statistically differentially expressed by ≥ two-fold between T4 treatments at one or more non-Day 0 sampling times. To complement this analysis, we used linear and quadratic regression to identify 542 and 709 genes that were differentially expressed by ≥ two-fold in the 5 and 50 nM T4 treatments, respectively.
We found that T4 concentration affected the timing of gene expression and the shape of temporal gene expression profiles. However, essentially all of the identified genes were similarly affected by 5 and 50 nM T4. We discuss genes and biological processes that appear to be common to salamander and anuran metamorphosis, and also highlight clear transcriptional differences. Our results show that gene expression in axolotls is diverse and precise, and that axolotls provide new insights about amphibian metamorphosis.
Amphibian metamorphosis is generally characterized by dramatic and conspicuous developmental changes that are necessary for larvae to function as terrestrial adults. The morphological, behavioral, and physiological changes that occur during metamorphosis are associated with increases in thyroid hormone (triiodothyronine, T3 and thyroxine, T4; TH) [1, 2] and RNA synthesis . These events are interconnected; at metamorphosis, tissue-specific concentrations of TH activate and repress transcriptional networks within target cells that in turn regulate new patterns of development . Many genes that are associated with molecular and morphological events during metamorphosis have been identified from studies of anurans, and in particular Xenopus laevis. In contrast, little is known about patterns of gene expression during salamander metamorphosis.
Although anuran and salamander metamorphosis are regulated by many of the same endocrine factors, there is considerable developmental variation between these groups. Most conspicuously, some salamanders do not undergo a complete metamorphosis in nature. These salamanders are called paedomorphs because they retain larval characteristics into the adult stage, and because genetic and phylogenetic evidence suggests that they evolve from metamorphic ancestors [5, 6]. Paedomorphosis in the Mexican axolotl (Ambystoma mexicanum) is associated with low hypothalamic-pituitary-thyroid (HPT) activity and differential sensitivity of tissues to TH that results in some cryptic biochemical and molecular changes, but not the complete suite of morphological changes seen in related, metamorphic tiger salamanders. Interestingly, A. mexicanum can be induced to undergo anatomical metamorphosis by administering TH and endocrine factors that function upstream of TH synthesis [7, 8]. The axolotl provides an excellent alternative to anuran systems because metamorphosis can be precisely induced and studied in juveniles or adults that are not developing toward a metamorphic outcome.
Functional genomic approaches are beginning to reshape the way transcription is conceptualized during amphibian metamorphosis [9–11]. The transcriptional program for tissue regression, remodeling, and organogenesis is significantly more complicated than was initially predicted for Xenopus [12–15]. Previously, we used microarray technology to identify keratin biomarkers for T4 induced metamorphosis in the integument (epidermis) of the Mexican axolotl . We showed that 50 nM T4 induces a complex transcriptional program and axolotls complete metamorphosis with no mortality. Interestingly, this T4 concentration is known to affect gene expression and mortality in anurans [16–18] and it is higher than T4 concentrations estimated in the serum of spontaneously metamorphosing salamanders. For example, Larras-Regard et al.  reported 28 nM as the maximum serum T4 level in Ambystoma tigrinum, a close relative of the axolotl that typically undergoes metamorphosis. To further investigate the effect of T4 concentration on induced metamorphosis in the Mexican axolotl, we report the results of a second microarray experiment that induced metamorphosis using a much lower concentration of T4 (5 nM). Using 5 and 50 nM T4 microarray datasets, we describe the temporal transcriptional response of T4 induced metamorphosis and specifically address the following question: Does T4 concentration affect morphological metamorphosis and gene expression in the axolotl? We discuss the biological significance of some of the differentially expressed genes (DEGs) that were identified and the relationship between salamander and anuran metamorphic gene expression programs.
Effect of T4 concentration on morphological metamorphosis
During T4 induced metamorphosis, Mexican axolotls progress through developmental stages (0–IV)  that are defined by the resorption of the upper and lower tailfins, dorsal ridge, and gills. We staged all axolotls after 0, 2, 12, and 28 days of T4 treatment. No metamorphic changes were observed after two days of T4 treatment and thus axolotls from both T4 treatments were assigned to Stage 0. At Day 12, morphological changes were observed in 50 nM T4 treated axolotls (Stages I and II) but 5 nM T4 treated axolotls were indistinguishable from control animals (Stage 0). At Day 28, axolotls reared in 50 nM T4 had fully resorbed tailfins and gills, and thus had completed morphological metamorphosis (Stage IV). Between Days 12 and 28, 5 nM T4 treated axolotls initiated metamorphosis but did not complete all morphological changes by Day 28 (Stage III). On average, individuals complete metamorphosis after 35 days in 5 nM T4 (unpublished data). Thus, a low concentration of T4 delays the initiation timing of morphological metamorphosis but not the length of the metamorphic period.
Gene expression in the absence of T4
Our first set of statistical analyses tested control axolotls for temporal changes in mRNA abundance that were independent of T4 treatment [see Additional files 1, 2]. Temporal changes are expected if patterns of transcription (gene expression) change significantly over time as salamanders mature, or if there are uncontrolled sources of experimental variation. After adjusting the false discovery rate (FDR) to 0.05, none of the probe-sets (genes) on the custom Ambystoma GeneChip were identified as significantly differentially abundant as a function of time. Thus, we found no statistical support for differential gene expression among control animals.
Gene expression in the presence of T4: 5 versus 50 nM
Genes identified as differentially abundant between 5 and 50 nM T4 at Day 2
Actin, alpha, cardiac muscle
Solute carrier family 7, member 3
Lipopolysaccharide-induced TNF factor
Regulator of G-protein signaling 5
Cadherin 1, type 1, E-cadherin
ATP-binding cassette, sub-family B, member 4
Solute carrier family 6, member 14
Phosphoenolpyruvate carboxykinase 1
Adipose differentiation-related protein
Chemokine ligand 5
Heterogeneous nuclear ribonucleoprotein L-like
Transmembrane protein 79
Solute carrier family 11, member 1
Mitogen-activated protein kinase kinase 3
Modeling the transcriptional response of genes during induced metamorphosis
The regression analyses identified a number of DEGs that only met both of our criteria (statistically significant and ≥ two-fold change) in one of the T4 concentrations (5 nM, n = 79; 50 nM, n = 246; Figure 1B). Of the 79 genes unique to 5 nM T4, 36 (46%) were statistically significant in 50 nM T4 but did not eclipse our two-fold change criteria. Of the remaining 43 genes unique to 5 nM T4, 25 exhibited ≥ two-fold change in 50 nM T4 at one or more sampling times but were not statistically significant. Inspection of the 50 nM T4 regression patterns and fold change data associated with the 79 genes unique to 5 nM T4 revealed that all of these genes exhibited similar directional trends (up versus down regulation) in 5 and 50 nM T4 [see Additional files 9, 10]. Of the 246 genes unique to 50 nM T4, 96 (39%) were statistically significant in 5 nM T4, but failed to eclipse our two-fold change criteria. An additional 48 genes exhibited ≥ two-fold change for at least one sampling time in 5 nM T4, but were not statistically significant. Inspection of the 5 nM T4 regression patterns and fold change data associated with the 246 genes unique to 50 nM T4 demonstrated that 209 (85%) of these genes exhibited similar directional trends in 5 and 50 nM T4 [see Additional files 11, 12]. An additional 22 genes unique to 50 nM T4 did not exhibit > 1.5 fold changes relative to Day 0 controls until Day 28 (at which time they were differentially regulated by ≥ two-fold), suggesting that they are expressed during the terminal stages of metamorphosis [see Additional files 11, 12]. Presumably, these genes were not detected in 5 nM T4 because we did not sample latter time points for this concentration. These results reiterate the point that essentially all genes identified by our study were similarly, directionally expressed in the 5 and 50 nM T4 treatments.
Genes categorized as QLVD and QLCU in both T4 treatments
Fold Change 5 nM
Fold Change 50 nM
Serpin peptidase inhibitor clade B member 10
Serpin peptidase inhibitor clade B member 2
WAP four-disulfide core domain 5
Bioinformatic comparison: axolotl versus Xenopus
Salamanders and anurans may express similar genes during amphibian metamorphosis. To test this idea, we compared a list of 'core' up regulated metamorphic genes from Xenopus  to DEGs identified from our study of axolotl. Of the 59 genes that were reported as differentially up regulated by ≥ 1.5 fold in limb, brain, tail, and intestine from metamorphosing Xenopus, 23 (39%) are represented by at least one of the 3688 probe-sets analyzed in our study. Of these, only two (FK506 binding protein 2 and glutamate-cysteine ligase modifier subunit) were identified as statistically significant and differentially expressed by ≥ two-fold in our study [see Additional files 13, 14]. FK506 binding protein 2 was up regulated in axolotl and Xenopus. However, glutamate-cysteine ligase modifier subunit was down regulated in axolotl and up regulated in Xenopus. Thus, < 5% of the 'core' DEGs that are commonly expressed among Xenopus tissues during metamorphosis, were identified as DEGs in our study using axolotl.
We also conducted a comprehensive bioinformatics comparison between DEGs from axolotl epidermis and DEGs from T3 induced Xenopus intestine . Gene expression similarities may exist between the Xenopus intestine and axolotl epidermis because both organs undergo extensive extracellular remodeling that is associated with apoptosis of larval epithelial cells and the proliferation and differentiation of adult cell types. The presumptive orthologs of 111 of the 820 non-redundant DEGs from our study correspond to DEGs from Xenopus intestine [see Additional files 13, 14]. Of these 111 genes, 50 (45%) exhibited the same direction of differential expression in axolotl epidermis and Xenopus intestine. This list includes genes that are known to be associated with metamorphic developmental processes in amphibians. For example, two MMPs (MMP9 and MMP13) that are associated with extra cellular matrix turnover were up regulated in Xenopus intestine and axolotl epidermis. However, other genes were regulated in opposite directions. For example, keratin 12 and keratin 15 were down regulated in axolotl epidermis but up regulated in Xenopus intestine. These results show that there are similarities and differences in gene expression between Xenopus and axolotls when comparing tissues that undergo similar remodeling processes.
Biological, technical, and statistical replication
Primer sequences used for Q-RT-PCR
Matrix metallopeptidase 13
Solute carrier family 31 member 1
Genes previously identified using 50 nM T4 that were not identified using 5 nM T4
Actin alpha 2
Actin alpha 2
Chormosome 20 open reading frame 92
Chromosome 21 open reading frame 33
Cytochrome C oxidase subunit VIB polypeptide 1
Cytochrome P450 family 2 subfamily c polypeptide 8
EPH receptor A1
Eukaryotic translation initiation factor 4E binding protein 3
Ferritin heavy polypeptide 1
High-mobility group 20 A
Hypothetical protein FLJ11151
LSM10 U7 small nuclear RNA associated
Mediator of RNA polymerase II transcription subunit 31 homolog
Myosin heavy polypeptide 11
Myosin light polypeptide 9
Nick-associated protein 1
Non-metastatic cells 7
Ras-related C3 botulinum toxin substrate 2
Ubiquitin specific peptidase 10
Paedomorphic Mexican axolotls can be induced to undergo metamorphosis by administering TH. We found that axolotls initiate metamorphosis at least one week earlier in 50 nM versus 5 nM T4 and complete morphological transformations in 28 days. The lower 5 nM T4 concentration was sufficient to induce metamorphosis but the initiation timing was delayed and this proportionally delayed the time to complete metamorphosis. The same sequence of morphological changes was observed between T4 treatments and the majority of DEGs were identified in both T4 treatments, although their expression profiles were temporally shifted. Nearly all DEGs exhibited similar directional trends between treatments, and the subset of genes that were directly compared between the T4 treatments exhibited similar relative abundances. Our results show extremely similar changes in gene expression and morphology in the axolotl when varying T4 by an order of magnitude. This is an interesting result because T4 concentrations within this range are toxic to anurans and are known to affect tissue-specific abundances of transcription factors that regulate metamorphic gene expression programs . Below, we discuss the axolotl's precise transcriptional response to the range of T4 concentrations examined in this study. We then discuss the epidermal gene expression program of the axolotl, noting gene expression similarities and differences between salamander and anuran metamorphosis.
TH levels are known to increase in larval amphibians as they mature and reach maximal levels during metamorphic climax. When the concentration of TH reaches a critical intracellular level, transcriptional changes are initiated that bring about new patterns of development. Because the TH concentration required to alter transcription is cell-specific, tissues are often described as having different sensitivities to TH. The sensitivity of cells to TH involves multiple factors that affect the intracellular concentration of TH and the ability of TH to affect transcription, which is determined in part by the number of nuclear TH binding sites (TRs) . Mexican axolotls have functional TRs , but TH levels are apparently too low to initiate metamorphosis . Direct hypothalamic application of T4, using a dose that is insufficient to initiate metamorphosis via intraperitoneal injection, is sufficient to initiate metamorphosis in the axolotl  and related paedomorphic tiger salamanders . Thus, axolotls are capable of synthesizing TH in sufficient quantities to initiate and complete metamorphosis. However, the pituitary doesn't release a sufficient amount of thyrotropin to trigger the metamorphic process . Axolotl epidermis can be stimulated to initiate metamorphic changes in vitro, in isolation of endogenously synthesized TH . Thus, the metamorphic timing delay that we observed between the 5 and 50 nM T4 treatments probably reflects the time required to autonomously activate gene expression within TH responsive cells of the epidermis and the time required to stimulate the HPT axis. Rosenkilde  showed that this latency period is TH concentration dependent and above a critical dose (37.5 nM T3) there is no variation in latency. After accounting for an estimated one-week difference in the initiation timing of metamorphosis between the T4 treatments, there was not a difference in the length of the metamorphic period. Thus, endogenous TH levels were functionally, if not quantitatively similar between 5 and 50 nM T4 treated axolotls after metamorphosis was initiated. This idea is also supported by the precise gene expression response that we observed between the T4 treatments: essentially all of the genes were expressed in the same direction, and a subset of genes that could be reliably compared showed the same magnitude of gene expression.
The precision of the transcriptional response between T4 treatments indicates that axolotls are surprisingly tolerant to T4 levels that dramatically affect anuran mortality and gene expression. Others before us have also noted the tolerance of axolotls to high levels of T4 [27, 28]. Because anuran metamorphosis involves a more extensive and integrated set of remodeling events that are accomplished over a shorter time frame, there may be greater overlap in the sensitivities of cells to TH among tissues that causes metamorphic remodeling events to occur out of sequence. The fact that salamander metamorphosis encompasses fewer morphological changes and that many of the changes are not as integrated (hindlimb development occurs months before tail metamorphosis) may explain why axolotls are so tolerant to high T4 concentrations. However, failure to observe an increase in thyroid hormone receptor beta transcription in axolotl suggests there may be fundamental regulatory differences between anuran and salamander metamorphosis.
The larval epidermis of axolotls is extensively remodeled during T4 induced metamorphosis [29–31]. Application of TH to paedomorphic axolotls induces many of the same epidermal changes that occur during natural and induced metamorphosis in anurans, including apoptosis of larval cells, proliferation of adult cell types, and epidermal keratinization. Our results show that TH induces a diversity of transcriptional changes that are associated with specific remodeling processes. We observed significant gene expression changes between the T4 treatments at Day 2, prior to observable morphological changes at the whole-organism level. Most of these genes were up regulated in the higher T4 concentration relative to the lower T4 concentration. Day 2 gene expression changes may reflect direct transcriptional activation via the binding of exogenous TH to TRs, which are functional in axolotls . For example, the human ortholog of phosphoenolpyruvate carboxykinase 1, a primary target for transcriptional regulation of gluconeogenesis, is known to have a thyroid hormone response element . Early up regulation of phosphoenolpyruvate carboxykinase 1, as well as fructose 1,6 bisphosphotase, glucose 6 phosphate dehydrogenase, and 70 kD heat shock protein 5 at Day 12, indicates a biochemical response at the cellular level that includes activation of key regulatory enzymes of the gluconeogenic pathway. This is an interesting finding for the epidermis because such responses are generally associated with hepatic cell functions. Several other interesting gene expression changes were detected at Day 2. These include ATP binding cassette, subfamily B, member 4, and transglutaminase 1. ATP binding cassette family genes are up regulated in mammals during epidermal lipid reorganization and keratinocyte differentiation , and transglutaminase 1 encodes an enzyme that functions in the formation of the cross-linked, cornified envelop of keratinocytes . The early expression of these genes is curious because keratinization is assumed to be a terminal differentiation event in the metamorphosis of amphibian epidermis. Our results suggest that the process of keratinization is initiated very early. As a final example, two proteins that are specific to the mammalian inner ear were identified as significantly down regulated: otogelin and otoancorin. The head epidermis of the axolotl contains mechanoreceptors that are homologous to hair cells of the mammalian ear . Our results suggest remodeling of these and other neural components in the axolotl skin at metamorphosis. There are many additional examples that could be highlighted from our gene lists that have not been previously discussed within the context of amphibian metamorphosis.
The most gene expression changes and the greatest changes in transcript abundances were observed at Day 12 in 50 nM T4. For example, keratin 14, a prototypical marker of proliferating keratinocytes in mammals , was up regulated 1146 fold in 50 nM T4. This also marks the time of the greatest morphological remodeling. After this time, gene expression levels of many genes decreased. Thus, as has been described in anurans, many gene expression changes in axolotl are transient, increasing initially and then decreasing. For example, apoptosis is activated and terminated during anuran  and salamander  metamorphosis to regulate the death and replacement of larval epithelial cells. When statistically significant genes were analyzed in the absence of a two-fold change criterion, genes that were transiently up regulated (i.e., exhibited QC profiles) were statistically associated with apoptosis and proteolysis functional ontologies (data not shown). As another example of the similarities between the metamorphic gene expression changes that occur in the epidermis of frogs and salamanders, we identified three distinct probe-sets with established orthologies to human uromodulin that are dramatically down regulated in the epidermis of metamorphosing A. mexicanum. Furlow et al.  have observed analogous results in X. laevis and have shown that Xenopus uromodulin orthologs are exclusively expressed in the apical cells of the larval epidermis. These and other genes that are similarly expressed between urodeles and anurans will provide useful biomarkers for comparative studies of metamorphosis between these two groups.
Our informatics comparison between DEGs identified from axolotl epidermis and Xenopus intestine identified > 100 genes that are commonly expressed in these organs during metamorphic remodeling. However, over half of these genes were differentially expressed in opposite directions in axolotl and Xenopus. For example, several genes associated with immune function (CD74 antigen, chemokine ligand 5, interferon regulatory factor 1, proteasome beta subunit 9, and class I-related major histocompatability complex) were down regulated in axolotl epidermis and up regulated in Xenopus intestine. This is not too surprising because it is well established that the axolotl immune system is fundamentally different from that of other vertebrates, including Xenopus . Additionally, genes expressed in opposite directions in these comparisons may reflect fundamental differences that exist between intestinal and epidermal remodeling. Genes that exhibited similar transcriptional patterns between Xenopus intestine and axolotl epidermis were associated with many different functions. For example, DNA methyltransferase 1 and 17-beta hydroxysteroid dehydrogenase 8 were down regulated in axolotl epidermis and Xenopus intestine during induced metamorphosis. In mammals, DNA methyltransferase 1 functions to maintain DNA methylation patterns that influence gene transcription  and 17-beta hydroxysteroid dehydrogenase 8 preferentially inactivates androgens and estrogens, [42, 43]. This later example suggests a transcriptional response to increase gonadal steroid hormone levels during epithelial remodeling in amphibians. As a final example, a presumptive ortholog to human keratin 24 (SRV_13498_s_at) that was up regulated by > 1000 fold in axolotls exposed to 50 nM T4 was also up regulated in Xenopus intestine, albeit by a comparatively modest four-fold increase. These comparisons emphasize similarities and differences in gene expression during metamorphic epithelial tissue remodeling in anurans and salamanders.
Recent microarray analyses of anurans and salamanders show that amphibian metamorphosis involves thousands of gene expression changes, involving many biological processes that have previously received little attention [9, 10]. Our results show similarities and differences in the metamorphic transcriptional programs of anurans and salamanders. We expected to identify similarly expressed genes because epidermis was included in anuran tissue preparations that were used for microarray analysis, and because tissue remodeling that occurs during metamorphosis appears to involve some evolutionarily conserved biological processes. We also expected to observe transcriptional differences because anuran and salamander lineages diverged > 300 million years ago. Our results suggest that amphibian metamorphosis cannot be fully understood from the study of a few anuran species. We show here that axolotls offer several advantages (inducible metamorphosis, robust transcriptional response, less complex integration of remodeling events) that can be exploited to provide complementary and novel perspectives on amphibian metamorphosis.
Study animals for microarray analyses
All procedures were conducted in accordance with the University of Kentucky Animal Care and Use Guidelines (IACUC #00609L2003). Salamanders (A. mexicanum) were obtained from a single genetic cross, using adults from an inbred strain. Embryos and larvae were reared individually at 20–22 C in 40% Holtfreter's solution. After hatching, larvae were fed freshly hatched brine shrimp (Artemia sp., Brine Shrimp Direct, Ogden, UT) napuli until they were large enough (3 weeks) to eat California blackworms (Lumbriculus sp., J.F. Enterprises, Oakdale, CA). At approximately eight months of age, 30 salamanders were randomly assigned to each of 10 different treatments and reared in 40% Holtfreter's solution with or without T4, (Sigma T2376, St. Louis, MO). One hundred μM stock T4 solutions were made as described in Page et al. . Five and 50 nM T4 were made fresh for each water change by mixing 2.5 or 25 mL of 100 μM stock with 40% Holfreter's solution to a final volume of 50 L. Water was changed every third day.
Skin tissue was collected from salamanders following 0, 2, 12, and 28 days of T4 treatment. These time points were sampled to test for early gene expression changes that might precede morphological metamorphosis, and because 28 days is a sufficient period for complete metamorphosis of 50 nM T4 induced A. mexicanum . To obtain tissue, salamanders were anesthetized in 0.01% benzocaine (Sigma, St. Louis, MO) and ≈ 1 cm2 of skin tissue was removed from the top of the head.
Total RNA was extracted for each tissue sample with TRIzol (Invitrogen, Carlsbad, CA) according to the manufacturer's protocol; additionally, RNA preparations were further purified using a Qiagen RNeasy mini column (Qiagen, Valencia, CA). UV spectrophotometry and a 2100 Agilent Bioanalyzer (Agilent Technologies, Santa Clara, CA) were used to quantify and qualify RNA preparations. Three high quality RNA isolations from each treatment and sampling time combination were used to make individual-specific pools of biotin labeled cRNA probes. Each of the 30 pools was subsequently hybridized to an independent GeneChip. The University of Kentucky Microarray Core Facility generated cRNA probes and performed hybridizations according to standard Affymetrix protocols.
A custom Ambystoma Affymetrix (Santa Clara, CA) GeneChip was designed from curated expressed sequence tag assemblies for A. mexicanum and A. tigrinum [44, 45]. The array contains 4844 probe-sets, 254 of which are controls or replicate probe-sets. Detailed descriptions of this microarray platform can be found in Page et al.  and Monaghan et al. .
Quality control and low-level analyses
We used the Bioconductor package affy  that is available for the statistical programming environment R  to perform a variety of quality control and preprocessing procedures at the individual probe level [49, 50]. These procedures included: 1) generating matrices of M versus A plots for all replicate GeneChips (n = 3 GeneChips for 10 treatment/sampling time combinations), 2) investigating measures of central tendency, measures of dispersion, and the distributions of all 30 GeneChips via boxplots and histograms, 3) viewing images of the log2(intensity) values for each GeneChip to check for spatial artifacts, and 4) viewing an RNA degradation plot  that allows for visualization of the 3' RNA labeling bias across all GeneChips simultaneously. In addition, we used ArrayAssist Lite software (Stratagene, La Jolla, CA) to assess several quality control measures that are recommended by Affymetrix such as average background (minimum = 48.36, maximum = 59.28, n = 30) scale factors (minimum = 0.404, maximum = 0.629, n = 30), and percent present (minimum = 81.5, maximum = 86.5, n = 30). Next, we processed our data by implementing the robust multi-array average (RMA) algorithm of Irizarry et al. .
Assessment of GeneChip precision
To obtain estimates of between GeneChip repeatability, we generated correlation matrices for the hybridization intensities across all probe-sets among replicate GeneChips. Very high and consistent mean r-values were calculated for each of the 10 treatment by sampling time combinations (range of mean r ± standard error = 0.966 ± 0.002 to 0.986 ± 0.001). These results demonstrate that we were able to obtain a high level of repeatability between replicate GeneChips. Our data are MIAME compliant and raw data files can be obtained at Sal-Site [45, 52].
Microarray platforms may not accurately or precisely quantify genes with low intensity values [53, 54]. Because low intensity genes contribute to the multiple testing problem that is inherent to all microarray studies, we filtered probe-sets whose mean expression values across all GeneChips (n = 30 per gene) were smaller than or equal to the mean of the lowest quartiles (25th percentiles) across all GeneChips (n = 30, mean = 6.53, standard deviation = 0.04; data presented on a log2 scale). Upon performing this filtration step, 3688 probe-sets were available for significance testing. We then performed PCA on the centered and scaled data from these probe-sets. This analysis allowed us to visualize the relationships between GeneChips within and across treatments.
Temporal gene expression in the absence of T4
We investigated whether genes exhibited differential expression as a function of time in the absence of T4 (control animals sampled at Days 0, 2, 12, and 28) via linear and quadratic regression . We corrected for multiple testing by evaluating α0 according to the algorithm of Benjamini and Hochberg  with a FDR of 0.05. α1 was set to 0.05.
Detecting and classifying DEGs
We conducted three analyses to investigate the effect of T4 on epidermal gene expression. For the first analysis, we used limma [57, 58] to identify genes that were differentially expressed as a function of T4 concentration. The limma package couples linear models with an empirical Bayes methodology to generate moderated t-statistics for each contrast of interest. This approach has the same effect as shrinking the variance towards a pooled estimate and thus reduces the probability of large test statistics arising due to underestimations of the sample variances. Operating limma requires the specification of two matrices. The first is a design matrix in which the rows represent arrays and the columns represent coefficients in the linear model. The second is a contrast matrix in which the columns represent contrasts of interest and the rows represent coefficients in the linear model. For this analysis, the design matrix specified a coefficient for each unique treatment by sampling time combination (10 coefficients) and the contrast matrix specified the calculation of contrasts between the two T4 concentrations (5 and 50 nM) at each of the non-zero sampling times (Days: 2, 12, and 28). In addition to moderated t-statistics, limma also generates moderated F-statistics. These moderated F-statistics test the null hypothesis that no differences exist among any of the contrasts specified by a given contrast matrix. A FDR correction  of 0.05 was applied to the P-values associated with the moderated F-statistics of the contrast matrix. In order to further reduce the number of false positives, we required that all "identified" DEGs be differentially regulated by ≥ two-fold at one or more of the contrasted time points.
The last two analyses were conducted using the regression-based approach of Liu et al.  to detect genes that exhibit differential expression as a function of days of T4 treatment. This approach also classifies DEGs into nine categories based on their temporal expression profiles as determined via linear and quadratic regression. In the context of our experiment, these profiles have specific biological interpretations. However, exceptions to these interpretations exist (see results for examples). In general, genes that exhibited LD, LU, QLCD, and QLVU expression profiles were still actively undergoing changes in their expression levels when the study was terminated. In contrast, genes that exhibited QLVD and QLCU expression profiles underwent down and up regulation respectively but reached steady state expression levels before the experiment was terminated. Finally, genes that exhibited QC and QV expression profiles were transiently up and down regulated respectively before returning to baseline expression levels. Null results are described by the 'Flat' category. Separate analyses were conducted for the 5 and 50 nM datasets with α0 evaluated at a FDR of 0.05 according to the algorithm of Benjamini and Hochberg  and α1 = 0.05. DEGs were required to exhibit ≥ two-fold changes relative to Day 0 controls at one or more sampling times before they were categorized as "identified".
Over representation analyses for genes with established orthologies
Biological process gene ontology categories that are over represented in our lists of DEGs (statistically significant and ≥ two-fold change) were identified using the Database for Annotation Visualization and Integrated Discovery (DAVID) . In all analyses, the 3085 probe-sets on the Ambystoma GeneChip with established orthologies were used as the background for generating expected values. The EASE threshold was always set to 0.05, and the count threshold was always set to two.
Bioinformatic comparison with Xenopus
In a recent microarray study, Buchholz et al.  presented "a core set of up regulated genes". These genes have been identified as up regulated in response to T3 treatment by ≥ 1.5 fold in every tissue that has been examined in metamorphosing X. laevis via microarray analysis (limb, brain, tail, and intestine) [9, 10]. We determined the orthologies of these genes to human as described in Page et al.  and Monaghan et al. . We then identified genes listed by Buchholz et al.  that were also differentially expressed in our experiment. The same approach was used to compare our gene lists to the 2340 DEGs identified by Buchholz et al.  from the intestine of metamorphosing X. laevis.
Biologically and technically independent verification
We conducted a second experiment using Q-RT-PCR to investigate the temporal expression patterns of five genes identified as differentially expressed by microarray analysis (Table 3). Animals used in our second experiment were raised as described for the animals used in the microarray experiment, with the exception that T4 treatment (50 nM) was initiated at 120 days post-fertilization. Tissue samples from two or three individuals were collected as described above beginning on Day 0 (prior to initiating T4 treatment) and were collected every two days for 32 days (i.e., Day 0, Day 2, Day 4... Day 32).
Total RNA was extracted from integument as described for the microarray experiment with the exception that all samples were treated with RNase-Free DNase Sets (Qiagen, Valencia, CA) according to the manufacturer's protocol. For each sample, the Bio-Rad iScript cDNA synthesis kit (Hercules, CA) was used to synthesize cDNA from 1 μg total RNA. Primers (Table 3) were designed using Primer3 , and design was targeted to the same gene regions that are covered by Affymetrix probe-sets. All PCRs were 25 μL reactions that contained: cDNA template 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 conditions were as follows: 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. Melting curve analysis was used to ensure the amplification of a single product for each reaction. All reactions were run in 96 well plates and blocked by sampling time (i.e., each of the 17 sampling times was equally represented for each gene on each plate). PCRs were performed using a Bio-Rad iCycler iQ Multi-Color Real Time PCR Detection System (Hercules, CA). All plates contained template free controls . Primer efficiencies were estimated via linear regression and relative expression ratios (R) were calculated according to Pfaffl . All expression ratios are relative to the average of the Day 0 animals, and normalized to transcriptional intermediary factor 1 (probe-set ID: L_s_at). This gene was selected as a control because it had an extremely small standard deviation across all treatment regimes in the microarray experiment (n = 30, mean = 14.44, standard deviation = 0.03; data presented on a log2 scale).
Statistical Analysis of the Q-RT-PCR Data
Log2 transformed R-values for each gene were analyzed separately via linear and quadratic regression models in which days of T4 treatment was the predictor variable. These analyses were carried out using JMP, Version 5 (SAS Institute, Cary, NC). We decided whether to use a linear or quadratic model for a given gene via forward selection . In short, quadratic models were accepted when the polynomial terms were significant (P < 0.05) and resulted in an increase in the proportion of variation in the data explained by the model (adjusted R2) of ≥ five percent relative to the linear model. The residuals of all models were inspected graphically. In addition, the residuals of all models were checked for normality. In cases where the assumption of normality was violated, regression analyses were run to obtain equations that describe the response of these genes to T4 as a function of time. However, such analyses were conducted with the understanding that a strict hypothesis testing interpretation could prove problematic.
Aspects of this project were supported by IBN-0242833 from the National Science Foundation (NSF) CAREER Award program and Grant Number 5-R24-RR016344-06 from the National Center for Research Resources (NCRR), a component of the NIH. Additional support for this research came from NIH Grant Numbers P20RR-016741 and P20RR-16481 from the INBRE Program of NCRR. The content of this study is solely the responsibility of the authors and does not necessarily represent the official views of NCRR or NIH or NSF. We thank Donna Walls and the University of Kentucky Microarray Core Facility for processing the GeneChips associated with this study. Arnold Stromberg provided statistical advice during the preparation of this manuscript.
- Leloup J, Buscaglia M: La triiodothyronine: hormone de la metamorphose des amphibiens. CR Acad Sci. 1977, 284: 2261-2263.Google Scholar
- Larras-Regard E, Taurog A, Dorris M: Plasma T4 and T3 levels in Ambystoma tigrinum at various stages of metamorphosis. Gen Comp Endocrinol. 1981, 43: 443-450.PubMedView ArticleGoogle Scholar
- Tata JR: Requirement for RNA and protein synthesis for induced regression of the tadpole tail in organ culture. Dev Biol. 1966, 13: 77-94.PubMedView ArticleGoogle Scholar
- Shi YB: Amphibian Metamorphosis: From Morphology to Molecular Biology. 2000, New York: Wiley-LissGoogle Scholar
- Shaffer HB, Voss SR: Phylogenetic and mechanistic analysis of a developmentally integrated character complex: alternative life history modes in ambystomatid salamanders. Am Zool. 1996, 36: 24-35.View ArticleGoogle Scholar
- Voss SR, Smith JJ: Evolution of salamander life cycles: a major-effect quantitative trait locus contributes to discrete and continuous variation for metamorphic timing. Genetics. 2005, 170: 275-281.PubMedPubMed CentralView ArticleGoogle Scholar
- Rosenkilde P, Ussing AP: What mechanisms control neoteny and regulate induced metamorphosis in urodeles?. Int J Dev Biol. 1996, 40: 665-673.PubMedGoogle Scholar
- Kuhn ER, De Groef B, Van der Geyten S, Darras VM: Corticotropin-releasing hormone-mediated metamorphosis in the neotenic axolotl Ambystoma mexicanum : synergistic involvement of thyroxine and corticoids on brain type II deiodinase. Gen Comp Endocrinol. 2005, 143: 75-81.PubMedView ArticleGoogle Scholar
- Das B, Cai L, Carter MG, Piao YL, Sharov AA, Ko MSH, Brown DD: Gene expression changes at metamorphosis induced by thyroid hormone in Xenopus laevis tadpoles. Dev Biol. 2006, 291: 342-355.PubMedView ArticleGoogle Scholar
- Buchholz DR, Heimeier RA, Das B, Washington T, Shi YB: Pairing morphology with gene expression in thyroid hormone-induced intestinal remodeling and identification of a core set of TH-induced genes across tadpole tissues. Dev Biol. 2007, 303: 576-590.PubMedView ArticleGoogle Scholar
- Page RB, Monaghan JR, Samuels AK, Smith JJ, Beachy CK, Voss SR: Microarray analysis identifies keratin loci as sensitive biomarkers for thyroid hormone disruption in the salamander Ambystoma mexicanum. Comp Biochem Physiol C Toxicol Pharmacol. 2007, 145: 15-27.PubMedView ArticleGoogle Scholar
- Wang Z, Brown DD: A gene expression screen. Proc Natl Acad Sci USA. 1991, 88: 11505-11509.PubMedPubMed CentralView ArticleGoogle Scholar
- Buckbinder L, Brown DD: Thyroid hormone-induced gene expression changes in the developing frog limb. J Biol Chem. 1992, 267: 25786-25791.PubMedGoogle Scholar
- Shi YB, Brown DD: The earliest changes in gene expression in tadpole intestine induced by thyroid hormone. J Biol Chem. 1993, 268: 20312-20317.PubMedGoogle Scholar
- Denver RJ, Pavgi S, Shi YB: Thyroid hormone-dependent gene expression program for Xenopus neural development. J Biol Chem. 1997, 272: 8179-8188.PubMedView ArticleGoogle Scholar
- Zhang F, Degitz SJ, Holcombe GW, Kosian PA, Tietge J, Veldhoen N, Helbing CC: Evaluation of gene expression endpoints in the context of a Xenopus laevis metamorphosis-based bioassay to detect thyroid hormone disruptors. Aquat Toxicol. 2006, 76: 24-36.PubMedView ArticleGoogle Scholar
- Helbing CC, Bailey CM, Ji L, Gunderson MP, Zhang F, Veldhoen N, Skirrow RC, Mu R, Lesperance M, Holcombe GW, Kosian PA, Tietge J, Korte JJ, Degitz SJ: Identification of gene expression indicators for thyroid axis disruption in a Xenopus laevis metamorphosis screening assay part 1. effects on the brain. Aquat Toxicol. 2007, 82: 227-241.PubMedView ArticleGoogle Scholar
- Helbing CC, Ji L, Bailey CM, Veldhoen N, Zhang F, Holcombe GW, Kosian PA, Tietge J, Korte JJ, Degitz SJ: Identification of gene expression indicators for thyroid axis disruption in a Xenopus laevis metamorphosis screening assay part 2. effects on the tail and hindlimb. Aquat Toxicol. 2007, 82: 215-226.PubMedView ArticleGoogle Scholar
- Cano-Martinez A, Vargas-Gonzalez A, Asai M: Metamorphic stages in Ambystoma mexicanum. Axolotl Newsletter. 1994, 23: 64-71.Google Scholar
- Zar JH: Biostatistical Analysis. 1999, Upper Saddle River, NJ: Prentice Hall, 4Google Scholar
- McGowan KM, Coulombe PA: Onset of keratin 17 expression coincides with the definition of major epithelial lineages during skin development. J Cell Biol. 1998, 143: 469-486.PubMedPubMed CentralView ArticleGoogle Scholar
- Yoshizato K, Frieden E: Increase in binding capacity for triiodothyronine in tadpole tail nuclei during metamorphosis. Nature. 1975, 254: 705-707.PubMedView ArticleGoogle Scholar
- Safi R, Bertrand S, Marchand O, Duffraisse M, de Luze A, Vanacker JM, Maraninchi M, Margotat A, Demeneix B, Laudet V: The axolotl (Ambystoma mexicanum), a neotenic amphibian, expresses functional thyroid hormone receptors. Endocrinology. 2004, 145: 760-772.PubMedView ArticleGoogle Scholar
- Kuhn ER, Jacobs GFM: Metamorphosis. Developmental Biology of the Axolotl. Edited by: Armstrong JB, Malacinski GM. 1989, New York: Oxford University Press, 187-197.Google Scholar
- Norris DO, Gern WA: Thyroxine-induced activation of hypothalamo-hyphophysial axis in neotenic salamander larvae. Science. 1976, 194: 525-527.PubMedView ArticleGoogle Scholar
- Takada M, Yai H, Komazaki S: Effect of calcium on development of amiloride-blockable Na+ transport in axolotl in vitro. Am J Physiol. 1998, 275: R69-R75.PubMedGoogle Scholar
- Rosenkilde P: The role of hormones in the regulation of amphibian metamorphosis. Metamorphosis. Edited by: Balls M, Bownes M. 1985, Oxford: Clarendon Press, 221-259.Google Scholar
- Brown DD: The role of thyroid hormone in zebrafish and axolotl development. Proc Natl Acad Sci USA. 1997, 94: 13011-13016.PubMedPubMed CentralView ArticleGoogle Scholar
- Fahrmann W: Die morphodynamik der epidermis des axolotls (Sirendon mexcianum Shaw) unter dem einflub von exogen appliziertem thyroxin I. die epidermis des neotenen axolotls. Z Mikrosk Anat Forsch. 1971, 83: 472-506.PubMedGoogle Scholar
- Fahrmann W: Die morphodynamik der epidermis des axolotls (Sirendon mexcianum Shaw) unter dem einflub von exogen appliziertem thyroxin II. die epidermis wahrend der metamorphose. Z Mikrosk Anat Forsch. 1971, 83: 535-568.PubMedGoogle Scholar
- Fahrmann W: Die morphodynamik der epidermis des axolotls (Sirendon mexcianum Shaw) unter dem einflub von exogen appliziertem thyroxin III. die epidermis des metamorphosierten axolotls. Z Mikrosk Anat Forsch. 1971, 84: 1-25.PubMedGoogle Scholar
- Giralt M, Park EA, Gurney AL, Liu J, Hakimi P, Hanson RW: Identification of a thyroid hormone response element in the phosphoenolpyruvate carboxykinase (GTP) gene. J Biol Chem. 1991, 266: 21991-21996.PubMedGoogle Scholar
- Kielar D, Kaminski WE, Liebisch G, Piehler A, Wenzel JJ, Mohle C, Heimerl S, Langmann T, Friedrich SO, Bottcher A, Barlage S, Drobnik W, Schmitz G: Adenosine triphosphate binding cassette (ABC) transporters are expressed and regulated during terminal keratinocyte differentiation: a potential role for ABCA7 in epidermal lipid reorganization. J Invest Dermatol. 2003, 121: 465-474.PubMedView ArticleGoogle Scholar
- Hitomi K: Transglutaminases in skin epidermis. Eur J Dermatol. 2005, 15: 313-319.PubMedGoogle Scholar
- Northcutt RG, Catania KC, Criley BB: Development of lateral line organs in the axolotl. J Comp Neurol. 1994, 340: 480-514.PubMedView ArticleGoogle Scholar
- Coulombe PA, Kopan R, Fuchs E: Expression of keratin K14 in the epidermis and hair follicle: insights into complex programs of differentiation. J Cell Biol. 1989, 109: 2295-2312.PubMedView ArticleGoogle Scholar
- Yoshizato K: Molecular mechanism and evolutional significance of epithelial-mesenchymal interactions in the body and tail-dependent metamorphic transformation of anuran larval skin. Int Rev Cytol. 2007, 260: 213-260.PubMedView ArticleGoogle Scholar
- Ohmura H, Wakahara M: Transformation of skin from larval to adult types in normally and metamorphosis-arrested salamander, Hynobius retardatus. Differentiation. 1998, 63 (5): 238-246.PubMedGoogle Scholar
- Furlow JD, Berry DL, Wang Z, Brown DD: A set of novel tadpole specific genes expressed only in the epidermis are down-regulated by thyroid hormone during Xenopus laevis metamorphosis. Dev Biol. 1997, 182: 284-298.PubMedView ArticleGoogle Scholar
- Volk H, Charlemagne J, Tournefier A, Ferrone S, Jost R, Parisot R, Kaufman J: Wide tissue distribution of axolotl class II molecules occurs independently of thyroxin. Immunogenetics. 1998, 47: 339-349.PubMedView ArticleGoogle Scholar
- Hermann A, Gowher H, Jeltsch A: Biochemistry and biology of mammalian DNA methyltransferases. Cell Mol Life Sci. 2004, 61: 2571-2587.PubMedView ArticleGoogle Scholar
- Fomitcheva J, Baker ME, Anderson E, Lee GY, Aziz N: Characterization of Ke 6, a new 17β-hydroxysteroid dehydrogenase, and its expression in gonadal tissue. J Biol Chem. 1998, 273: 22664-22671.PubMedView ArticleGoogle Scholar
- Aziz N, Anderson E, Lee GY, Woo DDL: Arrested testis development in the cpk mouse may be the result of abnormal steroid metabolism. Mol Cell Endocrinol. 2001, 171: 83-88.PubMedView ArticleGoogle Scholar
- Putta S, Smith JJ, Walker JA, Rondet M, Weisrock DW, Monaghan J, Samuels AK, Kump K, King DC, Maness NJ, Habermann B, Tanaka E, Bryant SV, Gardiner DM, Parichy DM, Voss SR: From biomedicine to natural history research: EST resources for ambystomatid salamanders. BMC Genomics. 2004, 5: 54-PubMedPubMed CentralView ArticleGoogle Scholar
- Smith JJ, Putta S, Walker JA, Kump DK, Samuels AK, Monaghan JR, Weisrock DW, Staben C, Voss SR: Sal-site: integrating new and existing ambystomatid salamander research and informational resources. BMC Genomics. 2005, 6: 181-PubMedPubMed CentralView ArticleGoogle Scholar
- Monaghan JR, Walker JA, Page RB, Putta S, Beachy CK, Voss SR: Early gene expression during spinal cord regeneration in the salamander Ambystoma mexicanum. J Neurochem. 2007, 101: 27-40.PubMedView ArticleGoogle Scholar
- Bioconductor: Open Source Software for Bioinformatics. [http://www.bioconductor.org]
- The R Project for Statistical Computing. [http://www.r-project.org]
- Bolstad BM, Irizarry RA, Gautier L, Wu Z: Preprocessing high-density oligonucleotide arrays. Bioinformatics and Computational Biology Solutions Using R and Bioconductor. Edited by: Gentleman R, Carey VJ, Huber W, Irizarry RA, Dudoit S. 2005, New York: Springer, 13-32.View ArticleGoogle Scholar
- Bolstad BM, Collin F, Brettschneider J, Simpson K, Cope L, Irizarry RA, Speed TP: Quality assessment of Affymetrix GeneChip data. Bioinformatics and Computational Biology Solutions Using R and Bioconductor. Edited by: Gentleman R, Carey VJ, Huber W, Irizarry RA, Dudoit S. 2005, New York: Springer, 33-47.View ArticleGoogle Scholar
- Irizarry RA, Hobbs B, Collin F, Beazer-Barclay YD, Antonellis KJ, Scherf U, Speed TP: Exploration, normalization, and summaries of high density oligonucleotide array probe level data. Biostatistics. 2003, 4: 249-264.PubMedView ArticleGoogle Scholar
- Sal-Site. [http://www.ambystoma.org]
- Choe SE, Boutros M, Michelson AM, Church GM, Halfon MS: Preferred analysis methods for Affymetrix GeneChips revealed by a wholly defined control dataset. Genome Biol. 2005, 6: R16-PubMedPubMed CentralView ArticleGoogle Scholar
- Draghici S, Khatri P, Eklund AC, Szallasi Z: Reliability and reproducibility issues in DNA microarray measurements. Trends Genet. 2006, 22: 101-109.PubMedPubMed CentralView ArticleGoogle Scholar
- Liu H, Tarima S, Borders AS, Getchell TV, Getchell ML, Stromberg AJ: Quadratic regression analysis for gene discovery and pattern recognition for non-cyclic short time-course microarray experiments. BMC Bioinformatics. 2005, 6: 106-PubMedPubMed CentralView ArticleGoogle Scholar
- Benjamini Y, Hochberg Y: Controlling the false discovery rate: a practical and powerful approach to multiple testing. J R Stat Soc Ser B. 1995, 57: 289-300.Google Scholar
- Smyth GK: Linear models and empirical Bayes methods for assessing differential expression in microarray experiments. Stat Appl Genet Mol Biol. 2004, 3: 3-Google Scholar
- Smyth GK: Limma: linear models for microarray data. Bioinformatics and Computational Biology Solutions Using R and Bioconductor. Edited by: Gentleman R, Carey VJ, Huber W, Irizarry RA, Dudoit S. 2005, New York: Springer, 397-420.View ArticleGoogle Scholar
- Dennis G, Sherman BT, Hosack DA, Yang J, Gao W, Lane HC, Lempicki RA: DAVID: database for annotation, visualization, and integrated discovery. Genome Biol. 2003, 4 (5): P3-Epub 2003 Apr 3.PubMedView ArticleGoogle Scholar
- Rozen S, Skaletsky H: Primer3 on the www for general users and for biologist programmers. Methods Mol Biol. 2000, 132: 365-386.PubMedGoogle Scholar
- Bustin SA, Nolan T: Pitfalls of quantitative real-time reverse-transcription polymerase chain reaction. J Biomol Tech. 2004, 15: 155-166.PubMedPubMed CentralGoogle Scholar
- Pfaffl MW: A new mathematical model for relative quantification in real-time RT-PCR. Nucleic Acids Res. 2001, 29: e45-PubMedPubMed CentralView ArticleGoogle Scholar