- Research article
- Open Access
Multiple stressors produce differential transcriptomic patterns in a stream-dwelling salamander
BMC Genomics volume 20, Article number: 482 (2019)
Global biodiversity is decreasing at an alarming rate and amphibians are at the forefront of this crisis. Understanding the factors that negatively impact amphibian populations and effectively monitoring their health are fundamental to addressing this epidemic. Plasma glucocorticoids are often used to assess stress in amphibians and other vertebrates, but these hormones can be extremely dynamic and impractical to quantify in small organisms. Transcriptomic responses to stress hormones in amphibians have been largely limited to laboratory models, and there have been few studies on vertebrates that have evaluated the impact of multiple stressors on patterns of gene expression. Here we examined the gene expression patterns in tail tissues of stream-dwelling salamanders (Eurycea tynerensis) chronically exposed to the stress hormone corticosterone under different temperature regimes.
We found unique transcriptional signatures for chronic corticosterone exposure that were independent of temperature variation. Several of the corticosterone responsive genes are known to be involved in immune system response (LY-6E), oxidative stress (GSTM2 and TRX), and tissue repair (A2M and FX). We also found many genes to be influenced by temperature (CIRBP, HSC71, HSP40, HSP90, HSP70, ZNF593). Furthermore, the expression patterns of some genes (GSTM2, LY-6E, UMOD, ZNF593, CIRBP, HSP90) show interactive effects of temperature and corticosterone exposure, compared to each treatment alone. Through a series of experiments we also showed that stressor induced patterns of expression were largely consistent across ages, life cycle modes, and tissue regeneration.
Outside of thermal stressors, the application of transcriptomes to monitor the health of non-human vertebrate systems has been vastly underinvestigated. Our study suggests that transcriptomic patterns harbor stressor specific signatures that can be highly informative for monitoring the diverse stressors of amphibian populations.
Environmental stressors such as habitat degradation, climate change, disease, and invasive species are central to the loss of global biodiversity [1, 2]. Effective conservation and management requires understanding when species cannot cope with environmental conditions . Organismal systems potentially express many signs of enduring stress, but the set of metrics commonly analyzed to assess the health of wildlife are typically limited .
In vertebrates, “stress” can activate the hypothalamic-pituitary-interrenal (or adrenal) axis (HPI-axis), leading to the production of elevated plasma glucocorticoids, and these hormones have been widely used for stress assessment [5,6,7,8]. However, glucocorticoid levels can be difficult or impossible to measure in some organisms [9, 10] and may not be indicative of underlying chronic (long-term) stress [11,12,13]. Genes that are regulated directly or indirectly by “stress hormones” or stressors themselves may offer a potentially rich source of informative biomarkers for monitoring population health and understanding adaptation. The transcriptomic responses of climatic variables, specifically temperature, have been well-studied for more than a decade [14, 15], but wildlife are subject to a wide range of stressors [1, 2], some of which may be exacerbated by changing climate [16,17,18]. Nevertheless, only a few studies have examined the transcriptional responses of animals to multiple stressors [18,19,20,21].
Amphibians are at the forefront of the biodiversity crisis and are likely threatened by multiple stressors (habitat loss, disease, climate change) [2, 22,23,24,25]. Many amphibian species appear to be declining in relatively healthy and undisturbed ecosystems, and the causative agents are often unknown [22, 25]. Most amphibians are small, and whole animals or pools of animals may need to be sacrificed in order to extract sufficient quantities of glucocorticoids for analysis via radioimmunoassay (e.g. ). Water-borne [27, 28], urine , fecal , and keratinized tissue  assays have recently been developed as non-invasive means to monitor glucocorticoid levels. However, daily , seasonal [33, 34], and life history variation [8, 35] can present additional challenges to understanding how point measurements of glucocorticoid levels relate to stress response. Gene expression responses in amphibians to “stress” hormones such as glucocorticoids have been studied extensively, but primarily focused on understanding how glucocorticoids regulate cellular, developmental, or physiological mechanisms [6, 36,37,38,39]. Furthermore, such studies have been based on model organisms (e.g. Xenopus). Only a few studies have evaluated the transcriptional responses of amphibians to climatic variables such as temperature [10, 18, 40].
Here we performed transcriptional analyses on adult stream-dwelling salamanders (Eurycea tynerensis) exposed to chronic corticosterone (a primary glucocorticoid) and different temperatures regimes. We tested whether transcriptional patterns provide genetic biomarkers for monitoring stress by simulating chronically elevated corticosterone and fluctuating temperature regimes. We further evaluated the robustness of several “stress response” genes across different life stages (adult and larvae), life histories (biphasic and paedomorphic), and to recent tissue regeneration (Table 1). Many of the differentially regulated genes were specific to temperature variation or corticosterone exposure and are known to be associated with cellular and physiological health in humans and biomedical models. Our analyses suggest that transcriptomic approaches may be key to understanding the diverse types of stressors that can impact amphibians.
Transcriptional responses to corticosterone exposure and temperature
Our reference transcriptome for Eurycea tynerensis included 4348 identified transcripts (median transcript length: 2207 bp; range: 214 to 17,654 bp; total length: 10,944,985 bp; N50: 3109 bp). Thirty-seven percent of 18.1 million 150 bp paired end RNA-Seq reads mapped to the reference and 4082 transcripts had sufficient coverage for analysis (see Methods section). Our RNA-Seq analyses revealed significant transcriptional responses in tail tissue of aquatic adult salamanders chronically exposed (30 days) to corticosterone at high and low temperatures. Out of 4082 genes, 11 were differentially expressed in response to corticosterone (Fig. 1a), which includes 8 upregulated and 3 downregulated genes (Additional file 2: Table S2). Substantially more genes (837) were differentially expressed by long-term exposure to high versus low temperatures (Fig. 1b). Exposure to high temperature (21 °C) upregulated 665 genes and downregulated 172 genes compared to maintenance at a low temperature (11 °C; Additional file 2: Table S2). Pairwise comparisons between low temperature control treatment and the other three treatments revealed that numerous genes share differential patterns of expression among stress conditions, however many are stressor specific (Fig. 1c; Additional file 2: Table S2). At low temperature, 14 genes were differentially regulated by corticosterone, whereas high temperature plus corticosterone differentially regulated 314 genes. High temperature without corticosterone differentially regulated 249 genes compared to low temperature control (Fig. 1c). Heat map of the 100 most differentially expressed genes shows discrete clustering based on treatments (Fig. 2).
Quantitative PCR analyses of six corticosterone and six temperature-regulated genes generally reflected our RNA-Seq analysis (Fig. 3). A2M, FX, GSTM2, and TRX were upregulated, and LY-6E and UMOD were downregulated in response to corticosterone. Temperature specific genes that were upregulated include ZNF593, CIRBP, HSC71, HSP40, and HSP90, whereas HSP70 was downregulated. The effects of corticosterone exposure on some genes (GSTM2, LY-6E, UMOD, ZNF593, CIRBP, HSP90) were temperature dependent (Fig. 3).
Transcriptional responses to corticosterone across life stages and life cycle modes
We found that larvae from both paedomorphic and biphasic populations exposed to corticosterone for 28 days had similar patterns of gene expression to one another (Fig. 4) and to paedomorphic adults (Fig. 3). LY-6E and UMOD were downregulated, and A2M and TRX were upregulated with corticosterone exposure (Fig. 4). However, corticosterone had a more dramatic impact on downregulating LY-6E and UMOD in larvae from biphasic populations compared to larvae from paedomorphic populations (Fig. 4). There was also significant upregulation of FX and GSTM2 in larvae from paedomorphic populations treated with corticosterone, unseen in larvae from biphasic populations. Despite some genes with life cycle dependent responses to chronic corticosterone exposure, other genes (LY-6E, UMOD, A2M, and TRX) had consistent patterns of expression across life cycle modes and life stages.
Transcriptional responses to corticosterone over time and during tissue regeneration
Circulating levels of glucocorticoids can be dynamic [32, 33, 41], and therefore measurements may be highly variable depending on very recent physiological conditions. In order to assess chronic stress, it may be important to identify markers that more consistently reflect signatures of physiological challenges. We found that the corticosterone induced upregulation of A2M and TRX and downregulation of LY-6E and UMOD were constitutively maintained over at least a two-week period (2 weeks to 4 weeks of exposure; Fig. 5). Corticosterone treatment of cultured tails for 48 h demonstrates a similar pattern of upregulation in GSTM2, A2M and FX (Fig. 6), which were also upregulated by chronic corticosterone treatment (30 day) of whole salamanders (Fig. 3). LY-6E was also upregulated by corticosterone in cultured tails, but is notably downregulated under chronic corticosterone treatment (Figs. 3 and 6).
We found that A2M, FX, and TRX were upregulated and LY-6E and UMOD were downregulated in regenerated tissues exposed to corticosterone for 14 days compared to controls (Fig. 5). This demonstrates that chronic exposure to corticosterone constitutively upregulates (A2M, FX, TRX) or downregulates (LY-6E, UMOD) genes across different stress durations and also during the process of regeneration.
Transcriptional responses to temperature variability
Adult paedomorphic and biphasic (metamorphosed) salamanders exposed to a thermally “stressful” environment with a 14 °C daily temperature fluctuation (21 ± 7 °C) for 4 weeks showed significant differential regulation of several genes compared to salamanders maintained at a constant temperature with the same mean (21 ± 0 °C). HSC71 was significantly upregulated under fluctuating temperature regimes in the adults of both paedomorphic and biphasic salamanders. HSPs showed life cycle specific responses to fluctuating temperatures, with HSP40 and HSP70 upregulated in metamorphs, and HSP90 upregulated in paedomorphs. Also, CIRBP was downregulated by temperature fluctuation in biphasic adults compared to those at stable temperatures (Fig. 7). With the exception of GSTM2, genes that were differentially regulated by chronic corticosterone (Fig. 3) did not differ in response to fluctuating temperatures (Fig. 7). In summary, some genes (e.g. GSTM2) are constitutively regulated by chronic corticosterone exposure and by fluctuating temperature. However, we identified almost a dozen genes with responses specific to corticosterone (A2M, FX, LY-6E, TRX, and UMOD) or fluctuations in temperature (CIRBP, HSC71, HSP40, HSP70, HSP90, and ZNF593).
The responses of species to changing climate are difficult to predict, and the environmental variables responsible for stress can be elusive. Compared to humans, there are relatively few metrics to assess the health of wildlife . For amphibians this has been largely limited to body-weight metrics [42, 43], corticosterone levels [27,28,29, 35, 44, 45], and pathogen presence [46, 47]. Therefore it is important to develop more diverse metrics for understanding amphibian stress response and population health. Transcriptional responses of amphibians to stress hormones such as corticosterone have mostly focused on acute stress in the laboratory model Xenopus [36,37,38, 48, 49]. Studies that have investigated gene expression patterns associated with stressors that may not necessarily engage the HPI-axis have been limited to temperature [10, 40] and disease resistance [18, 50]. Here we show that chronic corticosterone exposure and temperature stress instigate unique transcriptomic patterns in the tail tips of a stream-dwelling salamander. Our study demonstrates an example of how transcriptomic data can provide useful information for understanding amphibian responses to diverse stressors in wild or captive populations.
Chronic activation of the HPI-axis and the subsequent production of glucocorticoids can have deleterious effects on immune function and increased vulnerability to disease [51,52,53,54]. We found that chronic corticosterone treatment altered the expression patterns of several genes involved in immune system response, oxidative stress, and tissue repair. In regard to the immune system, the upregulation of LY-6E is thought to be part of a compensatory mechanism against pathological dysfunction following infection [55, 56]. We found that LY-6E was strongly downregulated after chronic corticosterone treatment (Figs. 3 and 5). In comparison, GSTM2 and TRX, which are known to combat cytotoxicity and oxidative stress [57, 58], were upregulated in our chronic corticosterone treatments of adult salamanders (Fig. 3). We also found that chronic corticosterone treatment upregulated FX and A2M, both of which are known to be involved in tissue repair [59,60,61]. Finally, UMOD (encoding Tamm-Horsfall glycoprotein) is best known for its diverse roles in the health of the mammalian kidney including immune functions and osmoregulation . However, this gene is also expressed in diverse amphibian tissues , but the functional role of UMOD like proteins in amphibians is still uncertain. Thus far, UMOD has been shown to be upregulated in tadpole facial tissues in the presence of predators , downregulated in the tails of metamorphosing tadpoles (gene 18 in ), and downregulated in the skin of larval salamanders treated with thyroxin . The expression of UMOD is strongly downregulated in response to corticosterone in the tail tips of adult and larval E. tynerensis in our study (Figs. 3 and 5).
Microarray analyses of Xenopus tails exposed to short-term (18 h) corticosterone showed differential regulation of 1968 genes . We reduced their dataset down to 501 known genes and compared it to a subset of 265 genes from our corticosterone versus control analysis based on uncorrected P-values (α < 0.05). Nine genes were shared between their acute and our chronic corticosterone treatments, but six of these genes displayed opposite patterns of expression (Additional file 2: Table S2). Opposing patterns of gene expression may result from species-specific, durational (“acute” vs. “chronic”), or context dependent responses. For example, we found that LY-6E was strongly upregulated in response to acute corticosterone treatment (Fig. 6), but downregulated after chronic treatment (Figs. 3 and 5). The transcriptional responses of amphibians to elevated glucocorticoids are still a nascent area of research. However, this will likely be a fruitful avenue for future investigation, and has the potential to provide a promising alternative for analyzing corticosterone-based stressors in amphibians.
As previously noted, not all stressors necessarily engage the HPI-axis , and therefore stress assays based on glucocorticoids alone could overlook other important factors that can impact the health of an organism. We found several corticosterone regulated genes lacked differences in expression when adult salamanders were faced with a thermal challenge (Fig. 7). Among the temperature-regulated genes, some HSPs have been previously identified as responsive to temperature stressors [40, 67] and infection . In particular, HSP40, HSP70, and HSP90 are well known for their roles in the cellular response to heat stress , and HSPs have been suggested as potential biomarkers of heat stress . We found that these genes, as well as HSC71, were upregulated in response to fluctuating temperature. Furthermore, HSPs have a variety of functions including minimizing cellular damage and maintaining homeostasis in a thermally fluctuating environment [68, 69].
We found that exposure to the stress hormone corticosterone or a thermally stressful environment produced unique and often consistent transcriptional patterns in salamanders. This suggests that transcriptomic patterns can be useful for monitoring the impact of diverse stressors. However, organisms are often simultaneously subject to multiple stressors , which can have additive or synergistic negative effects on amphibian health and survival [71,72,73,74]. Several of the genes that were assessed showed interactions between corticosterone exposure and heat stress (Table 2). The influence of corticosterone on GSTM2, LY-6E, UMOD, CIRBP, ZNF593, and HSP90 appears to be temperature dependent, whereas the expression of A2M, FX, CIRBP, and HSP40 appears to be influenced only by corticosterone or temperature.
Transcriptomes provide a promising yet under-evaluated source of information for understanding the health of amphibian populations and identifying specific stressors. We show that several genes express consistent patterns of corticosterone or temperature regulation across life stages, life cycle modes, and even during tissue regeneration. We are not suggesting that the genes that we found to be differentially expressed by temperature and corticosterone exposure in E. tynerensis will necessarily exhibit the same patterns in other species. The expression response of a given gene to stress may be species or clades specific, and may also show differences among sexes and tissues (e.g. [37, 75]). Furthermore, just like other biomarkers, there are a number of factors that may need to be resolved for each system before meaningful interpretation can be drawn about wild populations [4, 15]. Initial experiments would ideally be conducted under controlled settings and would entail exposing individuals to one or more stressors and biopsying an accessible tissue for transcriptomic analyses. This is to identify candidate genes that could be further evaluated under a range of contexts (stages, sexes, etc.). Ontogenetic assessment will be particularly important for amphibians that undergo dramatic metamorphosis between life stages, which involves significant changes to endogenous glucocorticoid levels [6, 34, 35] and gene expression patterns [37, 38]. If consistent stress response patterns ultimately emerge across clades, then qPCR of previously identified candidate genes for related species, instead of transcriptomics, could be used to reduce costs of development and assessment.
When facing a stressful environment, the immediate response of a healthy organism is often to compensate for or protect against the stressor . This can involve behavioral, physiological, and molecular responses, which normally subside once the stressor is removed. The impact of the stressor(s) on the health of the organism can vary based on intensity and duration. One of the greatest challenges in deciphering the patterns of any biomarker is to know when a value represents a normal (healthy) response to an acute stressor versus when the system is compromised by severe and/or chronic stress . Furthermore, organisms also exhibit geographic variation in baseline levels of stress hormones  as well as transcriptomes , and can evolve (adapt) to changing conditions . Therefore, geographic variation and adaptation may not necessarily indicate chronic stress. These aspects need to be considered when developing assays for stress, especially across wide ranging species.
Amphibians are negatively responding to ever increasing environmental stressors such as habitat loss, disease, and climate change. Using a transcriptomic approach, we identified a panel of genes that consistently and persistently responded to exposure of the “stress” hormone corticosterone and temperature variation across developmental stages, life cycle modes, and during tissue regeneration. While the specific genes identified here may be only relevant to this species, our study suggests that transcriptomics could be used to identify suites of genes that are indicative of the health of wild amphibian populations. Integrating transcriptomic analyses with other metrics of population stress and health expands the toolkit for conservation and management for understanding the factors that lead to amphibian declines.
The Oklahoma Salamander (Eurycea tynerensis) inhabits small streams in the Ozark Plateau of east-central North America, and exhibits alternative life cycle modes. Most populations of E. tynerensis have aquatic larvae that metamorphose into terrestrial adults (biphasic), whereas others forego metamorphosis and maintain their aquatic larval morphology and ecology into adulthood (paedomorphosis) [78, 79].
Some of our experiments are based on wild-caught E. tynerensis, while other experiments are based on F1 generation larvae and adults raised under controlled conditions in the lab. Prior to each experiment, both wild-caught and captive raised salamanders were acclimated at 18 °C for at least 2 days. For all experiments, larval and paedomorphic salamanders were maintained individually in 500 ml of their assigned solution; biphasic salamanders were kept on wet paper towels. Each experiment was conducted in incubators and experimental solutions were replaced and salamanders were fed bloodworms (chironomid larvae) every other day. All food provided was consumed across all experiments. To avoid disturbance the specimens were intentionally isolated, but this meant that no behavioral data were collected.
At the conclusion of the transcriptome experiment, salamanders were euthanized in a 0.1% solution of tricane methanosulfate (MS-222). For all other experiments salamanders were anesthetized by immersion in a 0.05% solution of MS-222 and awoken with dechlorinated tap water. The tail-tip (< 10% of tail) was dissected off, snap-frozen on dry ice, and stored at − 80 °C. Salamander care, maintenance, and experimentation were approved by the University of Tulsa (IACUC protocol TU-0028), and all experiments were performed in accordance with this protocol.
Organisms experience a variety of stressors and may respond by expressing unique quantifiable symptoms. For example, temperature stress induces the expression of Heat Shock Proteins (HSPs) [14, 80] that may be independent of the HPI-axis. To evaluate whether genes show differential transcriptional responses to specific stressors, our first experiment was designed to test whether chronic corticosterone treatment and different temperature regimes would provide transcriptional signatures in a conveniently biopsied tissue (salamander tail tip ; Table 1). Wild caught adult male paedomorphic E. tynerensis (N = 24) collected from the same locality were randomly split between incubators set at either 11 or 21 °C. Paedomorphic populations of E. tynerensis are adapted to relatively cool streams, and the population used in this experiment is from a groundwater fed stream with an average temperature of 13 °C (range 8 °C to 19 °C; Treglia et al. in prep). These salamanders tend to move to cooler microhabitats, deeper into the streambed when summer temperatures reach upper limits. Over the course of the month long experiment, salamanders kept at 11 °C maintained their body weight, while the body weights of those kept at 21 °C were reduced by ~ 16%. Therefore, 21 °C is above their normal temperature range and was considered a thermal stressor, while 11 °C was not considered stressful.
Within each temperature regime, half were exogenously treated with 100 nM corticosterone for 30 days. This dose is within the upper range or slightly above plasma corticosterone concentrations measured in other salamanders [44, 45]. Ethanol was used as a vehicle for corticosterone, so an equivalent amount of ethanol (< 0.001%) was added to control water. Due to the small size of the animals (average 370 mg), we were unable to obtain sufficient blood plasma to estimate circulating corticosterone levels at the conclusion of the experiment. However, larva and paedomorphic amphibians breathe through their porous skin and gills, and amphibians bathed in exogenous corticosterone solution readily up-take this hormone into their system (reviewed in [26, 81]). We replaced the corticosterone solution every other day across the course of the experiment to ensure a continuous dose.
Variation in physiological processes can be highly dependent upon age; therefore, gene expression patterns of adults may be different from larvae and juveniles [8, 35]. Furthermore, amphibians often exhibit variation in life cycle patterns , including discrete polymorphisms as observed in E. tynerensis [78, 79]. We conducted a series of experiments to validate the efficacy of corticosterone-regulated genes under different stages and life cycle modes using lab-raised larvae from paedomorphic and biphasic populations. Finally, wild amphibians commonly lose and regenerate their tails. Therefore, ideal biomarkers should produce consistent patterns even when tissues have been regenerated. We tested whether corticosterone induced transcriptional patterns were reproducible in newly regenerated tissues by analyzing tail clips that were regenerated while chronically exposed to corticosterone for 2 weeks.
During these experiments larvae from paedomorphic (N = 12) and biphasic (N = 18) populations were borne and raised in the laboratory at 21 °C. These larvae were exposed to either 100 nM of constant corticosterone or control (filtered water) for 28 days. After the first 14 days, 6 larvae from biphasic populations were anesthetized by immersion in MS-222 and their tail tips were biopsied for gene expression, and the salamanders were returned to 100 nM corticosterone treatment. This provided an earlier time point of corticosterone exposure (2 weeks), and also to evaluate the effects of corticosterone on expression patterns during regeneration. After 14 more days (at 28 days from the start of the experiment) tail tips were removed from all salamanders including the regenerated portion of the tails that were previously biopsied at 14 days.
We also conducted a tissue culture experiment on excised tail tips from 12 lab-raised adult, but non-reproductive (18 month old) paedomorphic E. tynerensis to test the effects of corticosterone on transcription when tissues are isolated from the rest of the endocrine system. The distal portions of tails (< 25% of total lengths) were cultured at 21 °C in 6-well plates and bathed in Leibovitz L-15 solution (diluted 2:1) with penicillin/streptomycin (100 units per ml). Tails were treated with either 100 nM corticosterone diluted in EtOH or an equivalent amount of EtOH as a control. Treatment solutions were replaced every 24 h. After 48 h of exposure to treatment conditions, tail tips were rinsed with 1x PBS and snap frozen on dry ice and stored at − 80 °C until RNA extraction.
We used paedomorphic (N = 9) and biphasic (N = 8) adult wild-caught salamanders to test for gene expression differences when exposed to a “stressful” thermal regime that involved dramatic daily changes in temperature. Salamanders were randomly split between 30 day temperature treatments, constant 21 °C or a thermally stressful fluctuating temperature regime with an average of 21 °C and a cyclical daily range of 14–28 °C.
RNA extraction, transcriptome sequencing, and quantitative PCR
RNA was isolated from tail tip tissue using Trizol Reagent (Invitrogen, Carlsbad, CA) following the manufacturer’s protocol. RNA concentrations were determined using either a QuBit fluorometer 2.0 (Thermofisher Scientific) for RNA-Seq samples, or a NanoDrop 8000 for samples that would be analyzed by quantitative PCR (qPCR). RNA-Seq libraries were prepared using the TruSeq RNA Library Preparation Kit (Illumina) and sequenced using 300 or 500 cycle V2 paired end read kits on an Illumina MiSeq at the University of Tulsa. All reads with a Q score less than 30 were discarded and adaptors were trimmed using MiSeq Reporter prior to analyses.
We have been iteratively building a partial Eurycea tynerensis transcriptome based on diverse tissues from several larval and adult individuals. These tissues included adult tail tips (N = 24), adult skin (N = 30), larval brains (N = 2), adult brains (N = 2), larval livers (N = 2), adult kidneys (N = 1), oviducts (N = 1), and testes (N = 1). We performed de novo assemblies of each tissue type and of individuals using CLC Genomics Workbench version 7.5.1 (Qiagen). These assemblies included a total of more than 100 million 150 to 250 bp paired-end Illumina reads, and were assembled with a similarity fraction of 0.95 or higher. Consensus sequences were extracted using a minimum coverage of 5x and ambiguity threshold of 0.25. We primarily identified genes by individually BLASTx searching transcripts against NCBI’s non-redundant protein database (parameters: organism = Vertebrata or Amphibian; max target sequence = 100; expected threshold = 10; max word size = 6; matrix = BLOSUM62; filter = low complexity regions). Groups of similarly identified transcripts were aligned using Clustal Omega and their uniqueness was evaluated by visually inspecting alignments. We identified 4348 transcripts with unique coding sequences (presumably non-redundant genes) totaling ~ 10.9 million bp, which was used as the reference for transcriptomic analyses.
We used the RNA-Seq function in CLC Genomic Workbench (95% similarity, 50% length fraction) to map 18.1 million, 150 bp paired-end, pass-filtered reads (~ 754 K reads per sample) to the partial E. tynerensis transcriptome (4348 genes). Our number of reference transcripts and depth of sequencing was sufficient to identify a large number of “highly expressed” genes that show large disparities among stress treatments. This made them readily quantifiable via qPCR. We used EdgeR  in the statistical platform R version 3.4.0  to identify differentially expressed genes between treatment groups (corticosterone or temperature) based on total read counts. To determine differentially expressed genes, we first reduced our initial 4348 genes to 4082 based on a minimum of 1 count per million across at least 6 of our 24 RNA-Seq libraries. The ‘calcNormFactors’ function was used to normalize each sample library based on scaling factors that minimize the log-fold changes between each sample. We used the ‘estimateDisp’ function to fit negative binomial models based on weighted likelihood empirical Bayes method to determine dispersion estimates for each sample. The ‘decideTests’ function was used to assess differential expression with a Benjamini-Hochberg adjusted P-value of 0.05 to minimize false discovery rates. A heat map was plotted using ‘hclust’  and ‘heatmap.2’ in R  to assess the degree of clustering among treatment groups.
TaqMan BHQ1a-6FAM qPCR assays were developed for 12 differently regulated genes (Additional file 1: Table S1; see Results section). This included six corticosterone regulated genes: Alpha-2 Macroglobulin (A2M), Coagulation Factor X-like (FX), Glutathione-S Transferase Mu 2 (GSTM2), Lymphocyte Antigen 6E (LY-6E), Thioredoxin (THIO), Uromodulin-like (UMOD), and six temperature regulated genes: Cold Inducible RNA Binding Protein (CIRBP), Zinc Finger 593 (ZNF593), Heat Shock Cognate 71 (HSC71), and the Heat Shock Proteins 40, 70, and 90 (HSP40, HSP70, and HSP90). cDNA was synthesized using SuperScript II (Invitrogen) and random hexamer. Reactions for qPCR were run with ABI TaqMan Gene Expression Master Mix on an ABI StepOne Plus (Thermofisher Scientific). Samples for a given gene were run simultaneously with a five-point standard curve, negative RT reactions, and negative controls. Expression quantity values were interpolated from CT values (number of cycles) based on the standard curves for each gene. Expression values were normalized with ribosomal protein L8 (rpL8), which is commonly used for normalization in amphibian gene expression studies [86, 87]. Relative gene expression values were log transformed and significant differences among groups were determined using ANOVA and a multiple-tests adjusted P-value (Benjamini and Hochberg method) in the R statistical platform .
Availability of data and materials
Public access to databases is open. RNA-Seq reads are available on Genbank (BioProject PRJNA531501), and log fold changes for each experiment and each gene are available as supplemental files.
Analysis of variance
Complimentary deoxyribonucleic acid
Hypothalamic pituitary interrenal-axis
Heat shock protein
quantitative Polymerase Chain Reaction
Ribonucleic acid sequencing
Pimm SL, Gareth JR, Gittleman JL, Brooks TM. The future of biodiversity. Science. 1995;269:347–50.
Hayes TB, Falso P, Gallipeau S, Stice M. The cause of global amphibian declines: a developmental endocrinologist’s perspective. J Exp Biol. 2010;213:921–33.
Busch DS, Hayward LS. Stress in a conservation context: a discussion of glucocorticoid actions and how levels change with conservation-relevant variables. Biol Conserv. 2009;142:2844–53.
Romero LM, Platts SH, Schoech SJ, Wada H, Crespi E, Martin LB, et al. Understanding stress in the healthy animal – potential paths for progress. Stress. 2015;18:491–7.
Sapolsky RM, Romero LM, Munck AU. How do glucocorticoids influence stress responses? Integrating permissive, suppressive, stimulatory, and preparative actions. Endocr Rev. 2000;21:55–89.
Denver RJ. Stress hormones mediate environment-genotype interactions during amphibian development. Gen Comp Endocrinol. 2009;164:20–31.
Carr JA. Stress and reproduction in amphibians. In: Norris DO, Lopez KH, editors. Hormones and reproduction of vertebrates, vol. 2. Cambridge: Academic Press; 2010. p. 99–116.
Crespi EJ, Williams TD, Jessop TS, Delehanty B. Life history and the ecology of stress: how do glucocorticoid hormones influence life-history variation in animals? Funct Ecol. 2013;27:93–106.
Romero LM, Fairhurst GD. Measuring corticosterone in feathers: strengths, limitations, and suggestions for the future. Comp Biochem Physiol A Mol Integr Physiol. 2016;202:112–22.
Cyzpionka T, Krugman T, Altmüller J, Blaustein L, Steinfartz S, Templeton AR, Nolte AW. Ecological transcriptomics – a non-lethal sampling approach for endangered fire salamanders. Methods Ecol Evol. 2015;6:1417–25.
Moore IT, Jessop TS. Stress, reproduction, and adrenocortical modulation in amphibians and reptiles. Horm Behav. 2003;43:39–47.
Dickens MJ, Romero LM. A consensus endocrine profile for chronically stressed wild animals does not exist. Gen Comp Endocrinol. 2013;191:177–89.
Woodley SK, Freeman P, Ricciardella LF. Environmental acidification is not associated with altered plasma corticosterone levels in the stream-side salamander, Desmognathus ochrophaeus. Gen Comp Endocrinol. 2014;201:8–15.
Sørenson JG. Application of heat shock protein expression for detecting natural adaptation and exposure to stress in natural populations. Curr Zool. 2010;56:703–13.
Alvarez M, Schrey AW, Richards CL. Ten years of transcriptomics in wild populations: what have we learned about their ecology and evolution? Mol Ecol. 2015;24:710–25.
Katzenberger M, Hammond J, Duarte H, Tejedo M, Calabuig C, Relyea RA. Swimming with predators and pesticides: how environmental stressors affect the thermal physiology of tadpoles. PLoS One. 2014;9:e98265.
Folt CL, Chen CY, Moore MV, Burnaford J. Synergism and antagonism among multiple stressors. Limnol Oceanogr. 1999;44:864–77.
Ribas L, Li MS, Doddington BJ, Robert J, Seidel JA, Kroll S, et al. Expression profiling the temperature-dependent amphibian response to infection by Batrachochytrium dendrobatis. PLoS One. 2009;4:e8408.
Chapman RW, Mancia A, Beal M, Veloso A, Rathburn C, Blair A, et al. The transcriptomic responses of the eastern oyster, Crassostrea virginica, to environmental conditions. Mol Ecol. 2011;20:1431–49.
Colburne JK, Pfrender ME, Gilbert D, Thomas WK, Tucker A, Oakley TH, et al. The ecoresponsive genome of Daphnia pulex. Science. 2011;331:555–61.
Maor-Landaw K, Ben-Asher HW, Karako-Lampert S, Salmon-Divon M, Prada F, Caroselli E, et al. Mediterranean versus Red Sea corals facing climate change, a transcriptome analysis. Sci Rep. 2017;7:42405.
Pechmann JHK, Wake DB. Enigmatic declines and disappearances of amphibian populations. In: Groom MJ, Meffe GK, Carroll CR, editors. Principles of conservation biology. 3rd ed. Sunderland: Sinauer Associates; 2006. p. 93–8.
Whitfield SM, Lips KR, Donnelly MA. Decline and conservation of amphibians in Central America. Copeia. 2016;104:351–79.
Grant EHC, Miller DAW, Schmidt BR, Adams MJ, Amburgey SM, Chambert T, et al. Quantitative evidence for the effects of multiple drivers on continental-scale amphibian declines. Sci Rep. 2016;6:25625.
Stuart SN, Chanson JS, Cox NA, Young BE, Rodrigues ASL, Fischman DL, et al. Status and trends of amphibian declines and extinctions worldwide. Science. 2004;306:1783–6.
Burraco P, Arribas R, Kulkarni SS, Buchholz DR, Gomez-Mestre I. Comparing techniques for measuring corticosterone in tadpoles. Curr Zool. 2015;61:835–45.
Gabor CR, Bosch J, Fries JN, Davis DR. A non-invasive water-borne hormone assay for amphibians. Amphib-Reptil. 2013;34:151–62.
Gabor CR, Zabierek KC, Kim DS, Barbiano LA, Mondelli M, Bendik N, et al. A non-invasive water-borne assay of stress hormones in aquatic salamanders. Copeia. 2016;104:172–81.
Narayan EJ. Non-invasive reproductive and stress endocrinology in amphibian conservation physiology. Conserv Physiol. 2013;1:cot11.
Woodruff JA, Lacey EA, Bentley G. Contrasting fecal corticosterone metabolite levels in captive and free-living colonial Tuco-Tucos (Ctenomys sociabilis). J Exp Zool. 2010;313:498–507.
Lattin CR, Reed M, DesRochers DW, Romero LM. Elevated corticosterone in feathers correlates with corticosterone-induce decreased feather quality: a validation study. J Avian Biol. 2011;42:247–52.
Bruener CW, Wingfield JC, Romero LM. Diel rhythms of basal and stress-induced corticosterone in a wild, seasonal vertebrate, Gambel’s white-crowned sparrow. J Exp Biol. 1999;284:334–42.
Romero LM. Seasonal changes in plasma glucocorticoid concentrations in free-living vertebrates. Gen Comp Endocrinol. 2002;128:1–24.
Glennemeier KA, Denver RJ. Developmental changes in interrenal responsiveness in anuran amphibians. Integr Comp Biol. 2002;42:565–72.
Chambers DL, Wojdak JM, Du P, Belden LK. Corticosterone level changes throughout larval development in the amphibians Rana sylvatica and Ambystoma jeffersonianum reared under laboratory, mesocosm, or free-living conditions. Copeia. 2011;2011:530–8.
Bonett RM, Hu F, Bagamasbad P, Denver RJ. Stressor and glucocorticoid induction of the immediate early gene Krüppel-like factor 9: implications for neural development and plasticity. Endocrinology. 2009;150:1757–65.
Bonett RM, Hoopfer ED, Denver RJ. Molecular mechanisms of corticosteroid synergy with thyroid hormone during tadpole metamorphosis. Gen Comp Endocrinol. 2010;168:209–19.
Kulkarni SS, Buchholz DR. Beyond synergy: corticosterone and thyroid hormone have numerous interaction effects of gene regulation in Xenopus tropicalis tadpoles. Endocrinology. 2012;153:5309–24.
Bagamasbad PD, Bonett RM, Sachs L, Buisine N, Raj S, Knoedler JR, et al. Deciphering the regulatory logic of an ancient, ultraconserved nuclear receptor enhancer module. Mol Endocrincol. 2015;29:856–72.
Sørenson JG, Pekkonen M, Lindgren B, Loeschcke V, Laurila A, Merilä J. Complex patterns of geographic variation in heat tolerance and Hsp70 expression levels in the common frog Rana temporaria. J Therm Biol. 2009;34:49–54.
Romero LM, Reed JM. Collecting baseline corticosterone samples in the field: is under 3 min good enough? Comp Biochem Physiol A Mol Integr Physiol. 2005;140:73–9.
Huntsman BM, Venarsky MP, Benstead JP, Huryn AD. Organic matter availability on the life history and production of a top vertebrate predator (Plethodontidae: Gyrinophilus palleucus) in two cave streams. Freshw Biol. 2011;56:1746–60.
Janin A, Léna JP, Joly P. Beyond occurrence: body condition and stress hormone as integrative indicators of habitat availability and fragmentation in the common toad. Biol Conserv. 2011;144:1008–16.
Homan RN, Reed JM, Romero LM. Corticosterone concentrations in free-living spotted salamanders (Ambystoma maculatum). Gen Comp Endocrinol. 2003;130:165–71.
Wack CL, DuRant SE, Hopkins WA, Lovern MB, Feldhoff RC, Woodley SK. Elevation of plasma corticosterone to physiologically relevant levels increased metabolic rate in a terrestrial salamander. Comp Biochem Physiol Part A: Mol Integr Physiol. 2012;161:153–8.
Greer AL, Collins JP. Sensitivity of a diagnostic test for amphibian ranavirus varies with sampling protocol. J Wildl Dis. 2007;43:525–32.
DiRenzo GV, Grant EHC, Longo AV, Che-Castaldo ZKR, Lips KR. Imperfect pathogen detection from non-invasive skin swabs biases disease inference. Methods. 2018;9:380–9.
Chen X, Li Y, Zhu M, Qin Z. An ex vivo assay for screening glucocorticoid signaling disruption based on glucocorticoid-response transcription in Xenopus tails. J Environ Sci. 2018;66:104–12.
Schneider KA, Shewade LH, Buisine N, Sachs LM, Buchholz DR. A novel stress hormone response gene in tadpoles of Xenopus tropicalis. Gen Comp Endocrinol. 2018;260:107–14.
Rosenblum EB, Poorten TJ, Settles M, Murdoch GK. Only skin deep: shared genetic responses to chytrid fungus in susceptible frog species. Mol Ecol. 2012;21:3110–20.
Marketon JIW, Glaser R. Stress hormones and immune function. Cell Immunol. 2008;252:16–26.
Cohen S, et al. Chronic stress, glucocorticoid receptor resistance, inflammation, and disease risk. Proc Natl Acad Sci U S A. 2012;109:5995–9.
Thomas JR, Woodley SK. Treatment with corticosterone delays wound healing in male and female salamanders. Gen Comp Endocrinol. 2015;216:33–8.
Fonner C, Patel S, Boord S, Venesky MD, Woodley SK. Effects of corticosterone on infection and disease in salamanders exposed to the amphibian fungal pathogen Batrachocytrium dendrobatidis. Dis Aquat Org. 2017;123:159–71.
Flood PM, Murphy DB, Horowitz M, LeClair KP, Smith FR, Stockert E, et al. A monoclonal antibody that recognizes an Ly-6-linked antigen inhibits the generation of functionally active T cell subsets. J Immunol. 1985;135:63–72.
Xu X, Qiu C, Zhu L, Huang J, Li L, Fu W, et al. IFN-stimulated gene LY6E in monocytes regulates the CD14/TLR4 pathway but inadequately restrains the hyperactivation of monocytes during chronic HIV-1 infection. J Immunol. 2014;193:4125–36.
Hayes J, Strange RC. Glutathione S-transferase polymorphisms and their biological consequences. Pharmacology. 2000;61:154–66.
Benhar M, Shytaj IL, Stamler JS, Savarino A. Dual targeting of the thioredoxin and glutathione systems in cancer and HIV. J Clin Invest. 2016;126:1630–9.
Armstrong PB, Quigley JP. α2-macroglobulin: an evolutionarily conserved arm of the innate immune system. Dev Comp Immunol. 1999;23:375–90.
Scotton CJ, Krupiczojc KM, Mercer PF, Lee GYC, Kaminski N, et al. Increased local expression of coagulation factor X contributes to the fibrotic response in human and murine lung injury. J Clin Invest. 2009;119:2250–563.
Rehman AA, Ahsan H, Khan FH. Alpha-2-macroglobulin: a physiological guardian. J Cell Physiol. 2013;228:1665–75.
Devuyst O, Olinger E, Rampoldi L. Uromodulin: from physiology to rare and complex kidney disorders. Nat Rev Nephrol. 2017;13:525–44.
Howie AJ, Lote CJ, Cunningham AA, Zaccone G, Fasulo S. Distribution of immunoreactive Tamm-Horsfall protein in various species in the vertebrate classes. Cell Tissue Res. 1993;274:15–9.
Mori T, Kawachi H, Imai C, Sugiyama M, Kurata Y, Kishida O, et al. Identification of a novel uromodulin-like gene related to predator-induced bulgy morph in anuran tadpoles by functional microarray analysis. PLoS One. 2009;4:e5936.
Brown DD, Wang Z, Furlow JD, Kanamori A, Schwartzman RA, Remo BF, et al. The thyroid hormone-induced tail resorption program during Xenopus laevis metamorphosis. Proc Natl Acad Sci U S A. 1996;93:1924–9.
Page RB, Monaghan JR, Walker JA, Voss SR. A model of transcriptional and morphological changes during thyroid hormone-induced metamorphosis of the axolotl. Gen Comp Endocrinol. 2009;162:219–32.
Heikkila JJ. Heat shock protein gene expression and function in amphibian model systems. Comp Biochem Physiol A Mol Integr Physiol. 2010;156:19–33.
Lindquist S. The heat-shock response. Annu Rev Biochem. 1986;55:1151–91.
Todgham AE, Iwama GK, Schulte PM. Effects of the natural tidal cycle and artificial temperature cycling on HSP levels in the tidepool sculpin Oligocottus maculosus. Physiol Biochem Zool. 2006;79:1033–45.
Semlitsch RD, Scott DE, Pechmann JHK, Gibbons W. Structure and dynamics of an amphibian community: evidence from a 16-year study of a natural pond. In: Cody ML, Smallwood JA, editors. Long-term studies of vertebrate communities. San Diego: Academic Press; 1996. p. 217–48.
Kiesecker JM. Synergism between trematode infection and pesticide exposure: a link to amphibian limb deformities in nature? Proc Natl Acad Sci U S A. 2002;99:9900–4.
Boone MD, Semlitsch RD, Little EE, Doyle MC. Multiple stressors in amphibian communities: effects of chemical contamination, bullfrogs, and public. Ecol Appl. 2007;17:291–301.
Chen CY, Hathaway KM, Thompson DG, Folt CL. Multiple stressor effects of herbicide, pH, and food on wetland zooplankton and a larval amphibian. Ecotoxicol Envrion Saf. 2008;71:209–18.
Hanlon SM, Parris MJ. The interactive effects of chytrid fungus, pesticides, and exposure timing on gray treefrogs (Hyla versicolor) larvae. Environ Toxicol Chem. 2014;33:216–22.
Calisi RM, Austin S, Lang A, MacManes M. Sex-biased transcriptomic response of the reproductive axis to stress. Horm Behav. 2018;100:56–68.
Eikenaar C, Husak J, Escallón C, Moore IT. Variation in testosterone and corticosterone in amphibians and reptiles: relationships with latitude, elevation, and breeding season length. Am Nat. 2012;180:642–54.
Schoville SD, Barreto FS, Moy GW, Wolff A, Burton RS. Investigating the molecular basis of local adaptation to thermal stress: population differences in gene expression across the transcriptome of the copepod Tigriopus californicus. BMC Evol Biol. 2012;12:170.
Bonett RM, Chippindale PT. Speciation, phylogeography and evolution of life history and morphology in the salamanders of the Eurycea multiplicate complex. Mol Ecol. 2004;13:1189–203.
Emel SL, Bonett RM. Considering alternative life history modes and genetic divergence in conservation: a case study of the Oklahoma salamander. Conserv Genet. 2011;12:1243–59.
Chen B, Feder ME, Kang L. Evolution of heat-shock protein expression underlying adaptive responses to environmental stress. Mol Ecol. 2018;27:3040–3054.
Belden LK, Moore IT, Wingfield JC, Blaustein AR. Corticosterone and growth in Pacific treefrog (Hyla regilla) tadpoles. Copeia. 2005;2005:324–30.
Duellman WJ, Trueb L. Biology of amphibians. 1st ed. Baltimore: Johns Hopkins University Press; 1994.
Robinson MD, McCarthy DJ, Smith GK. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Biostats. 2010;26:139–40.
R Core Team. R: a language and environment for statistical computing. Vienna: R Foundation for Statistical Computing; 2017. Available online at: URL https://www.r-project.org
Warnes GR, Bolker B, Bonebakker L, Gentleman R, Huber W, Liaw A, et al. Gplots: various R programming tools for plotting data. R package version 3.0.1. 2016.
Mazon RG, Denver RJ. Regulation of pituitary thyrotropin gene expression during Xenopus metamorphosis: negative feedback is functional throughout metamorphosis. J Endocrinol. 2004;182:273–85.
Aran RP, Steffen MA, Martin SD, Lopez OI, Bonett RM. Reduced effects of thyroid hormone on gene expression and metamorphosis in a paedomorphic plethodontid salamander. J Exp Zool B Mol Dev Evol. 2014;322:294–303.
We thank A. Boardman, M. Gibson, and L. Lehman for help in maintaining salamanders, and M. Herrboldt, A. Hess, N. Ledbetter, and R. Tovar for assistance with labwork. We thank K. Irwin of the Arkansas Game of Fish and M. Howery of the Oklahoma Department of Wildlife for granting permits for salamanders in their respective states.
This research was funded by Oklahoma NSF-EPSCoR (IIA-1301789 to RMB). Our funding agency did not play a role in the study design, data collection, analysis, interpretation of the data, or preparation of the manuscript.
Ethics approval and consent to participate
This study and experiments within were approved by the University of Tulsa’s Institutional Animal Care and Use Committee (IACUC) with permit #TU-0028.
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
About this article
Cite this article
Clay, T.A., Steffen, M.A., Treglia, M.L. et al. Multiple stressors produce differential transcriptomic patterns in a stream-dwelling salamander. BMC Genomics 20, 482 (2019) doi:10.1186/s12864-019-5814-y