Differential gene expression mediates physiological responses to perceived predation risk in a developmentally plastic vertebrate, the northern leopard frog (Lithobates pipiens)

Abstract

Prey can respond to predation risk through developmental plasticity, generating anti-predator phenotypes. These inducible defenses arise from changes to the stress axis, and neuroendocrine-triggered gene regulation is a likely mechanism influencing such phenotypes. As tadpoles, amphibians improve their escape performance by modifying tail shape in response to perceived predation risk (PPR), and this process should involve tissue and developmentally specific gene regulation. We exposed Lithobates pipiens tadpoles to PPR from Aeshnidae predators and measured tail morphology and transcriptomic response across different tissues (head and tail) and development (pre-metamorphosis to pro-metamorphosis). We found that PPR induced plasticity in tail shape, and this response was suppressed when tadpoles were also exposed to a glucocorticoid synthesis inhibitor. Differential gene expression was associated with predation stress across head and tail tissue, and developmental stage. Predator-exposed tadpoles exhibited up-regulation of genes responsible for muscle tissue and nervous system development, primarily in tail tissue and in pre-metamorphosis. PPR broadly influenced pathways across tissues and metamorphosis, including developmental, endocrine, and immune system pathways. This study provides an important step in understanding transcriptomic responses during predator induced morphological change, and demonstrates that gene expression, as induced by perceived predation risk, is a prominent mechanism of developmental plasticity.

Introduction

Predation is an important force shaping ecosystems, and non-consumptive effects of predators are increasingly recognized as contributing to changes in prey distribution, abundance and evolution (Preisser et al., 2005; Verdolin, 2006; Creel and Christianson, 2008). Prey can detect predators or the risk of predation through visual, olfactory, auditory, or other cues in the environment, and these cues can elicit anti-predator responses such as shifts in development (Gilbert, 2001; Beckerman et al., 2007), morphology (McCollum and Leimberger, 1997; Petrusek et al., 2009), physiology (Mateo, 2010; Clinchy et al., 2013), and behavior (Lima and Dill, 1990; Lind and Cresswell, 2005). Anti-predator responses should be under strong selection if they benefit prey survival, productivity, and fitness (Harvell, 1990; Buskirk and Relyea, 1998). However, the induction of such responses can be costly to prey (Pettersson and Brönmark, 1997; DeWitt et al., 1998; Relyea, 2002), meaning that selection should favor defenses that are phenotypically plastic and responsive to environmental cues (Harvell, 1990). Thus, inducible defenses should be expressed primarily under threat (Relyea, 2004; Kishida et al., 2007) and activated in response to cues such as prey alarm pheromones or predator kairomones (Fraker et al., 2009; Schoeppner and Relyea, 2009; Hettyey et al., 2015).

Cues from predation events, or perceived predation risk (PPR)can trigger the vertebrate neuroendocrine stress axis, and stress hormones can then mediate induction of phenotypically plastic responses to predation risk (Hossie et al., 2010; Middlemis Maher et al., 2013; Vinterstare et al., 2020). This implies a mechanistic pathway between the stress response and corresponding morphological or developmental plasticity (Denver, 2009, 2021), with gene expression potentially playing a key role in this process through stress hormone mediation of gene regulation (Baniahmad and Tsai, 1993; Paul et al., 2022). PPR can stimulate changes in gene expression in adult organisms (Schwarzenberger et al., 2009; Sanogo et al., 2011), and during embryonic development (Mommer and Bell, 2014; Tills et al., 2018). In invertebrates, PPR exposure promotes up-regulation of genes potentially supporting anti-predator defense morphogenesis (Miyakawa et al., 2010; Gu et al., 2022). However, to date, relationships between gene expression, the stress response, and inducible anti-predator defenses are poorly understood in vertebrates. It is therefore crucial to address what transcriptomic pathways might be activated by PPR, and how gene regulation, as a product of the vertebrate stress axis, varies according to tissue-specific needs of different prey phenotypes or at different developmental stages.

Understanding the complete transcriptomic response to PPR requires profiling both relationships between the neuroendocrine stress axis, and tissue-specific plastic responses. The neuroendocrine stress response in brain tissue and endocrine organs induces the subsequent production of stress hormones that can influence gene regulation in relevant body tissues (Baniahmad and Tsai, 1993; Paul et al., 2022). Therefore, gene expression in brain tissue during PPR should involve the neuroendocrine stress axis, and hormone production (Roseboom et al., 2007; Tyler et al., 2020), while genes expressed in tissues that form anti-predator phenotypes should be primarily related to physiological or morphological pathways such as muscle tissue development. For example, major candidate stress-related genes expressed in the brain during PPR are correlated with induction of anti-predator morphologies (Vinterstare et al., 2023). Accordingly, relevant gene expression in morphologically plastic tissues should also be evident during the stress response, but this relationship is not well established in vertebrates. Although PPR can activate key regulatory genes associated with the stress response and morphogenesis (Miyakawa et al., 2010; Vinterstare et al., 2023), these transcriptomic processes have not been described in concert, requiring investigation into broader transcriptomic and functional pathways that move beyond a targeted gene approach.

Animals with high developmental plasticity are especially predisposed toward stress-induced morphogenesis and can reveal differential gene expression as a fundamental molecular mechanism of plasticity (Beldade et al., 2011). Amphibian larvae exhibit high developmental plasticity and are sensitive to a variety of stressors during development, including PPR (Van Buskirk and Schmidt, 2000; Relyea, 2001; Teplitsky et al., 2005; Van Buskirk and Arioli, 2005; Kishida and Nishimura, 2006; Lent and Babbitt, 2020; Melotto et al., 2021). Amphibian tadpoles respond to PPR via morphological plasticity in tail depth (Dayton et al., 2005; Kishida and Nishimura, 2006; Steiner, 2007; Hossie et al., 2017; Sergio et al., 2021), which improves swimming speed and escape performance (Dayton et al., 2005), and directs predator attacks away from more vulnerable body regions (Van Buskirk et al., 2003). The presence of corticosterone (CORT) during development is an important mediator of larval amphibian tail shape plasticity (Hossie et al., 2010; Middlemis Maher et al., 2013), with increased tail depth being generated through exogenous addition of CORT alone, even without PPR (Middlemis Maher et al., 2013; Phuge et al., 2020). This suggests that the neuroendocrine response is developmentally influential and that corticosteroid stress hormones activated from the hypothalamo-pituitary-interrenal (HPI) axis can regulate gene expression in plastic body tissues through glucocorticoid receptor signaling (Baniahmad and Tsai, 1993; Paul et al., 2022).

In addition to the influence of CORT on development, thyroid hormones from the hypothalamo-pituitary-thyroid (HPT) axis have a strong mediation on metamorphosis (White and Nicoll, 1981; Paul et al., 2022) and increase steadily across development, transitioning tadpoles from pre- to pro-metamorphosis (Regard et al., 1978; Kikuyama et al., 1993). As this underlying endocrinological change is time dependent, we can expect stress to be differentially influential across developmental stages, especially considering the known interactions between the HPI and HPT axes (Denver, 2021; Paul et al., 2022). Currently, PPR-induced gene expression in amphibians is known to involve morphologically relevant functions such as blood vessel formation, keratin filament production, and epidermal osmoregulation (Mori et al., 2005, 2009, 2015), but past research has focused on limited targeted gene approaches that leave gaps in our understanding of crucial pathways involved in anti-predator responses. Accordingly, full transcriptome sequencing can provide higher sensitivity to detect a greater diversity and abundance of transcripts in order to reveal informative patterns of gene expression across developmental and neuroendocrine pathways (Mutz et al., 2013).

We examined the effects of predation risk on morphological and transcriptomic responses of northern leopard frog (Lithobates pipiens) tadpoles via experimental manipulation of predator cues and by measuring morphological change and related gene expression in head and tail tissues and across two developmental time periods. Tadpole head tissues contained the brain as well as relevant endocrine organs (pituitary, interrenal and thyroid glands), where neuroendocrine responses occur, and tadpole tail tissue constitutes the site of developmental plasticity. We hypothesize that transcriptional changes mediate the tadpole’s developmentally plastic response to PPR. First, we confirmed presence of anti-predator responses and involvement of neuroendocrine stress pathways by exposing tadpoles to PPR and predicting: 1) increased depth in tail shape following exposure, and 2) suppression of this response when the stress hormone pathway was experimentally blocked. Secondly, we predicted that PPR exposure elicits differential gene expression between: 3) head and tail tissues, and 4) pre- and pro-metaphoric tadpoles, with differentially expressed transcripts being annotated to genes associated with tadpole responses to predation, including: 5) morphological change in the tail tissue, 6) neuroendocrine stress responses, and 7) signatures of PPR exposure over metamorphosis. This study advances our knowledge of the transcriptomic effects of predation stress on a developmentally plastic vertebrate, thereby enhancing our understanding of phenotypic plasticity in response to environmental stressors as well as how chronic stress can influence developmental changes across the highly conserved vertebrate stress axis.

Materials and methods

Morphological responses experiment

This study extends our previous work on northern leopard frog tadpole transcriptomics (Row et al., 2016) and focuses on tadpole morphological and transcriptomic responses to predation risk. All animal care procedures were approved and conducted in accordance with the Trent University Animal Care Committee, under protocol #13005. Two experiments were conducted in parallel, to avoid potential cross-over effects of handling stress during morphological measurements on transcriptomic responses. In Spring 2013, we collected three northern leopard frog egg masses near Peterborough, ON (44 ° 36’N -78 ° 26’W). Egg masses were hatched, and tadpoles were reared in large outdoor mesocosms until they reached Gosner Stage (GS) 24 (~ 3 weeks), when they were transferred to an indoor laboratory and housed in 8L aquaria. Tadpoles (n = 480) were randomly assigned to one of three experimental treatments at GS 24, and each treatment consisted of eight replicates of 20 tadpoles per replicate. A PPR treatment, with tadpoles exposed to 50 mL of water used to rear tadpole predators, late-instar dragonfly nymphs (Aeshnidae), that were fed two tadpoles three times per week (Bennett et al., 2013, 2016). A Control group, with tadpoles exposed to 50 mL of reverse osmosis water as a sham; and A PPR + MTP treatment, involving application of the same predator chemical cue as in the PPR treatment along with metyrapone (MTP) dissolved in a water and ethanol solution to a concentration of 110 µM (Glennemeier and Denver, 2002b). Final ethanol concentration was < 0.005% of total water volume. MTP inhibits corticosteroid synthesis in vertebrates (Glennemeier and Denver, 2002b) and is widely used in experiments assessing the physiology of stress responses (Glennemeier and Denver, 2002c; Hossie et al., 2010; Middlemis Maher et al., 2013; Fraker et al., 2021). The PPR + MTP treatment was included to evaluate the role of the stress axis in eliciting tadpole morphological responses to perceived risk; adding MTP should counteract PPR application by blocking CORT, with tadpoles in the PPR + MTP treatment presumably being unresponsive to PPR and thus morphologically similar to control tadpoles. Such a response would support that gene expression related to the neuroendocrine response should be present in tadpoles exposed to PPR in the transcriptomics experiment. Previous work shows that in the doses we used for this experiment, MTP effectively reduces whole-body CORT by > 50% in L. pipiens tadpoles (Glennemeier and Denver, 2002b). Furthermore, no toxic effects have been detected when applying this concentration of MTP in tadpoles (Glennemeier and Denver, 2002b, 2002c).We also did not include an ethanol control because no influences on tadpole development or the stress axis have been detected when applying the associated concentration of ethanol with MTP treatments (Glennemeier and Denver, 2002c). PPR exposure began at GS 24, and the experiment was divided into early (6 weeks) and late (13 weeks) termination dates, at which point tadpoles were measured for their morphometric responses to treatment, corresponding to ranges of development at GS 26-31 (6 weeks – pre-metamorphosis) and GS 36-41 (13 weeks – pro-metamorphosis), respectively (Glennemeier and Denver, 2002a). We applied the PPR, sham, and PPR + MTP treatments three times per week.

Transcriptomic responses experiment

We measured tadpole gene expression in response to PPR exposure using the identical aforementioned study design including PPR and Control (sham application) treatments only. After 6 weeks (pre-metamorphosis), three tadpoles from different tank replicates of each treatment were randomly selected, euthanized and dissected, with head and tail tissue being isolated to assess gene expression in neuroendocrine-relevant (brain and endocrine organs) and morphology-relevant (tail) tissues. At 13 weeks (pro-metamorphosis) the process was repeated with remaining tadpoles so that differential gene expression could be compared between both developmental periods. Thus, 24 tissue samples were selected, which included six PPR-exposed samples and six control samples, three at each developmental period (6 weeks and 13 weeks), and samples from head and tail for each tadpole.

Morphometric analysis

The morphometric experiment was completed during week 6 and week 13, with tadpoles being removed from aquaria and photographed laterally in a standardized plexiglass frame (n = 303). Tadpole morphology was determined via landmarked images on the body and tail (Ferland-Raymond and Murray, 2008), and we defined tadpole shape with the coordinate data via a generalized Procrustes analysis (GPA) using the ‘Geomorph’ R package (Adams and Otárola-Castillo, 2013). Tadpole morphology was compared between treatments using principal component analysis of a covariance matrix of the GPA shape variables (Adams and Otárola-Castillo, 2013). We excluded PC1 because it usually represents positional and size variation rather than shape variation (McCollum and Van Buskirk, 1996; McCollum and Leimberger, 1997; Touchon and Warkentin, 2008; Collins and Gazley, 2017), and we generated schematic deformation diagrams of tadpole tails to reflect morphological variation represented by each principal component (Adams and Otárola-Castillo, 2013). PC2 was selected as the shape variable of interest, reflecting tail depth, a relevant anti-predator phenotype across many anurans (McCollum and Leimberger, 1997; Dayton et al., 2005; Hossie et al., 2017). We used a Linear Mixed-Effects model to assess individual variation in tail depth morphology (PC2) across treatments and developmental stages, with individual tank serving as a random effect. We confirmed normality and homoscedasticity of the tail depth PC2 scores with Shapiro-Wilk (Shapiro and Wilk, 1965) and Breusch-Pagan (Breusch and Pagan, 1979) tests, respectively. Our categorical variables were coded as deviations from the grand mean, using sum-to-zero contrast coding (Kaufman and Sweet, 1974). We present p-values generated from the model using a Type III sum of squares and a Kenward-Roger approximation (Kenward and Roger, 1997) (Supplementary Table 1).

RNA-sequencing and de novo transcriptome assembly

After total RNA extraction, we sequenced samples that passed initial quality control estimates (RNA Integrity number > 6.0), with 2 tadpoles yielding low quality tail tissue, leaving 22 total transcriptomes (12 head, and 10 tail) for RNA-sequencing (Supplementary Table 2). We sequenced samples on two lanes of an Illumina HiSeq 2000 using 100 bp paired-end reads. We quality-checked raw Illumina reads using FastQC v0.11.9 (Andrews, 2010) and filtered poor quality reads using the Trimmomatic v0.36 (Bolger et al., 2014) module in Trinity v2.12.0 (parameters: TruSeq3-PE.fa:2:30:10, slidingwindow:4:5, leading:5, trailing:5, minlength:25) (Grabherr et al., 2011). We performed a de novo assembly of transcript fragments (transfrags) of quality trimmed reads into transcriptomes, using Trinity v2.12.0 (parameters: default settings in strand-specific mode, RF; Grabherr et al., 2011). We estimated the quality of the assembled transcriptomes with the N50 and E90N50 metrics, using Trinity utilities. We evaluated completeness of the transfrag assembly with BUSCO v5.2.2, using transcriptome mode and auto-lineage eukaryotic search (Waterhouse et al., 2018).

Transcript abundance

We aligned quality trimmed reads from each library to the de novo transfrag assembly with the Trinity “align_and_estimate_abundance” perl script, where Bowtie2 v2.4.4 (Langmead and Salzberg, 2012) mapped the paired end reads and RSEM v1.3.3 (Li and Dewey, 2011) and estimated transfrag abundance in each library. The resulting expression matrices used weighted trimmed mean (TMM) of the log expression ratios (Grabherr et al., 2011). To assess similarities between different libraries based on transcript abundance, we used the Trinity “PtR” perl script on the RSEM-generated isoform counts matrix, where hierarchical clustering is used to generate a Pearson correlation matrix. This script was also used to generate a principal component analysis across all sample replicates (Grabherr et al., 2011).

Differential gene expression

We ran EdgeR v3.34 (Robinson et al., 2010) and DESeq2 v1.32 (Love et al., 2014) using the SARools v1.7.4 package (Varet et al., 2016) with rounded RSEM gene counts as input. These differential gene expression analyses compared expression of PPR vs. Control group samples from our different sample groups. We performed pairwise comparisons individually by groupings of tissue type and exposure duration (Supplementary Table 3). For all analyses, the unexposed Control samples were set as the reference condition and PPR samples were set to the factor of interest. Blocking factor was set to “Null”. EdgeR was run with settings alpha = 0.05, pAdjustMethod = BH, cpmCutoff = 1, and normalizationMethod = TMM, while DESeq2 was run with settings cooksCutoff = True, independantFiltering = True, alpha = 0.05, pAdjustMethod = BH, and locfunc = median (Varet et al., 2016). For both analyses we corrected p-values for multiple comparisons using the Benjamini-Hochberg method to control for false discovery rate (Benjamini and Hochberg, 1995). The DESeq2 analyses revealed more differentially expressed transcripts than EdgeR. Additionally, of the total set of differentially expressed transcripts identified by EdgeR, all but three were also identified by the DESeq2 analysis. We have thus restricted further analyses and interpretations to the DESeq2 results only. To address the use of wild-caught animals for this study (Todd et al., 2016), we estimated the statistical power of our tests of differential gene expression using RnaSeqSampleSize v3.19 within each of our treatment groups (Zhao et al., 2018).

Functional annotation

We performed functional annotation on the differentially expressed transfrags having a fold change > 2 for up-regulation, and < 0.5 for down-regulation (Schurch et al., 2016; Todd et al., 2016). For each, the longest open reading frame (ORF) was determined using TransDecoder v5.50 (Haas, 2018). The longest ORFs were queried using HMMER v 3.2.1 (Finn et al., 2011) to identify protein domains, TMHMM v2.0 (Krogh et al., 2001) to predict transmembrane helices, and SignalIP v4.1 (Petersen et al., 2011) to predict signal peptides. Blast+ v2.11.0 (Johnson et al., 2008) was used to query the Swissprot (Boeckmann et al., 2003), non-redundant (nr) proteomes (Pruitt et al., 2005), and the amphibian model organism Xenopus tropicalis proteome UP000008143 (Hellsten et al., 2010). We ran blastp queries of these databases on the longest ORFs and blastx queries on significantly differentially expressed transfrags. Functional annotation results were loaded into a Sqlite3 database (Gaffney et al., 2022) in order to generate a combined annotation report with Trinotate v3.2.2 (Bryant et al., 2017). We then extracted gene ontology (GO) assignments at the gene level from the Trinotate annotation report. To perform GO enrichment analysis, GOseq was run using the Trinity “run_GOseq” perl script (Grabherr et al., 2011). We ran GOseq individually by group (6 week tail, 6 week head, 13 week tail, 13 week head) using the GO assignments generated by Trinotate. The GOseq gene lengths input file was generated by the RSEM TMM normalized expression matrix, and the background input file was the full set of assembled transcripts at the gene level (Grabherr et al., 2011). Gene products may be associated with multiple biological processes, and many annotated transcripts may show redundant functional categories. Therefore, we used REVIGO (Supek et al., 2011) to help summarize the full list of enriched GO terms and generate a reduced set of the most relevant functions of gene products differentially expressed following exposure to PPR. We have further restricted our analysis to only examine GO terms which represent relevant biological processes. We also used REVIGO to generate semantic clustering plots which use Multidimensional scaling to cluster GO terms based on the SimRel (Schlicker et al., 2006) functional similarity measurement (Supek et al., 2011).

Results

Tadpole tail morphology

Our morphometric analysis revealed that PPR treatment affected tadpole tail shape, with tail depth (PC2) morphology varying across developmental stages (p=0.029), and treatment groups (p= 0.001), compared to controls. At 6 weeks, PPR-exposed tadpoles had deeper tail shape compared to both Control and PPR + MTP tadpoles (Figure 1), and tail shape for the PPR+MTP group was similar to the Control group (Figure 1). At 13 weeks, PPR-exposed tadpoles continued to have deeper tail morphology than Control tadpoles, although Control tadpoles developed a deeper tail shape over time. PPR + MTP tadpoles did not develop a deeper tail shape over time, even at 13 weeks (Figure 1). Therefore, morphological measurements revealed that PPR exposure: 1) elicited inducible defense morphology during pre-metamorphosis, but 2) in pro-metamorphosis, morphological responses to PPR treatment were dampened relative to Controls; and 3) PPR + MTP elicited responses that more closely align with Controls during pre-metamorphosis, thereby confirming the involvement of the stress axis in tadpole anti-predator responses in our experimental system.

Figure 1

Sequence reads, transcriptome quality assessment and abundance

We generated a total of 408,355,745 paired-end reads across 12 tadpoles, with separate head and tail transcriptomes for 22 individual/tissue combinations. RNA-seq libraries averaged 18,561,624 reads (range: 15,476,753 – 22,193,962), and the de novo transfrag Trinity assembly contained 333,449 genes and 502,104 transcripts. Quality checks of transfrag length revealed 50% of assembled sequences contained in transfrags were >1,296 bp in length (N50). More importantly, of the most highly expressed genes (90% of total normalized expression data representing 27,446 total genes), the mean gene length was 1,938 bp (E90N50). BUSCO analysis comparing the assembled L. pipiens transcriptomes to 255 core eukaryote genes revealed that 98.1% were complete, 0.8% fragmented, and 1.1% missing. The average alignment of individual RNA-seq libraries to the assembled transcriptome was 99.2% (Supplementary Table 4). The correlation matrix (Supplementary Figure 1) resolved sample type through hierarchical clustering, with the highest correlation within groups of head and tail tissue, indicating tissue-specific correspondence in transcriptomes. Similarly, principal component analysis revealed that tissue type explained most variation among samples (PC1 = 33.88%). Interestingly, transcriptomes also showed some clustering based on duration of PPR exposure, indicating tadpole transcriptomes had higher correlation within similar periods of metamorphosis (PC2 = 7.71%; Figure 2). Clustering associated with treatment exposure time and developmental stage were more pronounced in tail tissue than in head tissue (Figure 2). These results indicate the variance between tissue types and between developmental stages is greater than the biological variance within these groups, and thus these are appropriate within-group comparisons for differential gene expression tests between PPR-treatments and control.

Figure 2

Differential gene expression

Tadpoles exposed to PPR treatment showed differential gene expression compared to Controls. Differentially expressed transcripts due to PPR treatment differed according to tissue type and exposure duration, with transcripts being both up- and down-regulated in each group (Table 1). A total of 661 transcripts were differentially expressed across all groupings, with tail tissue showing more differentially expressed transcripts than head tissue (Tail - 366, Head - 295). Longer exposure to PPR led to more differential gene expression (6 weeks - 234, 13 weeks - 427), and within groupings, 13 week tails showed the most differential expression (264), followed by 13 week heads (163), 6 week heads (132), and 6 week tails (102). Of the 661 differentially expressed transcripts, 634 were unique to their tissue and duration of exposure grouping, and 27 transcripts were expressed in more than one grouping. Statistical power within our treatments was greater than might be expected from wild animals (Todd et al., 2016) with fold changes between 3.3 and 4.2 for upregulation, or 0.30 to 0.23 for downregulation, providing a minimum power of 0.8 (Supplementary Figure 2).

Functional annotation

Functional information was obtained from GO terms enriched to the corresponding gene products that were differentially expressed during PPR treatment. Annotation information was available for 32.4% of differentially expressed genes, and this rate is standard when using Xenopus to identify gene homologies in other amphibian species (Price et al., 2015; Yang et al., 2016; Liedtke et al., 2019). The GO terms showed 2,628 functional categories across differentially expressed genes of all groupings (6 week head= 398, 6 week tail= 440, 13 week head= 662, 13 week tail= 1,128). We present here a subset of functional pathways associated with differential expression responding to PPR treatment, reduced by removing many redundant or uninformative GO terms.

Overall, PPR exposed tadpoles displayed both up- and down-regulation of genes which were enriched to a variety of functional categories, including developmental, endocrine, and immune system processes. Enriched GO terms, reduced by REVIGO, were grouped into their broader “parent term” functional categories to show common enriched pathways across tissue groupings (Figure 3; Supplementary Table 5). We have also highlighted a selection of the REVIGO output focusing only on enrichment from genes associated with biological processes, including genes supporting muscle growth in tail tissue, glucocorticoid and thyroid hormone functions, and a diverse set of immune and stress response functions (Table 2). Developmental and immune system processes were common enriched pathways across all four groupings (Figure 3). We found developmental functions indicative of morphological change enriched in 6 week tail tissue (Figure 3; Table 2). Hormone-related pathways were enriched in 6 week heads, 13 week heads, and in 6 week tail tissue (Figure 3; Table 2). Immune system processes showed a diverse response across all tissue groupings (Figure 3), and a number of immune system genes were down-regulated in 13 week tail tissue (Table 2).

Figure 3

Pre-metamorphosis (GS 26-31): 6-week head

Head tissue of tadpoles exposed to 6 weeks of PPR showed up-regulation of the ZN205 gene, which is associated with programed cell death (GO:0043068; Table 2). Down-regulated genes were associated with such pathways as endocrine process, developmental processes, immune system processes, and others (Figure 3A). Specifically, down-regulation of HVM13 was associated with enriched immune system processes such as regulation of B cell activation (GO:0050864; Table 2). Genes RETNB and RETNG, which code for Resistin-like molecules, were down-regulated and associated with enriched endocrine pathways such as positive regulation of steroid hormone secretion (GO:2000833) and positive regulation of progesterone secretion (GO:2000872; Table 2).

Pre-metamorphosis (GS 26-31): 6-week tail

Tail tissue at 6 weeks showed enrichment to a variety of functional categories, including response to hormone, response to stress, development processes, immune system processes, and others. Developmental processes included pathways associated with anatomical growth, such as nervous system development, vasculature development, and muscle tissue development (Figure 3B). Specifically, up-regulation of the c-Fos gene was observed and associated with enriched functional pathways such as skeletal muscle cell differentiation (GO:0035914) and nervous system development (GO:0007399; Table 2). Likewise, myelin expression factor 2 was up-regulated and is associated with the enriched pathways of myotube differentiation (GO:0014902), neuron differentiation (GO:0030182), and striated muscles cell differentiation (GO:0051146; Table 2). 6-week tail tissue also showed up-regulation of the transcription factor JunB, which is associated with enriched endocrine functions such as response to progesterone (GO:0032570) and response to corticosterone (GO:0051412; Table 2). 6-week tail tissue of tadpoles exposed to PPR also displayed changes in transcripts related to immune and stress responses. For example, up-regulation of genes PXN1 and FUCL6 were associated with the enriched pathways regulation of humoral immune response (GO:0002920) and regulation of acute inflammatory response (GO:0002673). 6-week tail tissue only showed down-regulation of gene MYH8, which was associated with muscle filament sliding (GO:0030049; Table 2).

Pro-metamorphosis (GS 36-41): 13-week head

Head tissue of tadpoles exposed to 13 weeks of PPR showed enrichment to pathways such as programed cell death, immune system processes, response to stress, and developmental processes (Figure 3C). Up-regulation of pulmonary surfactant-associated protein (SFTPA) was observed and is associated with surfactant homeostasis (GO:0043129) from the gene SFTPA (Table 2). A broad range of up-regulated genes were associated with enriched immune functions, including innate immune response (GO:0045824), humoral immune response (GO:0006959), and positive regulation of B cell activation (GO:0050871; Table 2). 13-week head tissue also showed down-regulation of genes related to both thyroid hormones and the thyroid gland. Genes responsible for the production of thyroglobulin were down-regulated and associated with enriched functional categories such as thyroid hormone generation (GO:0006590) and thyroid gland development (GO:0030878; Table 2).

Pro-metamorphosis (GS 36-41): 13-week tail

13-week tail tissue of tadpoles exposed to PPR showed a variety of enriched functional categories, including response to stress, programmed cell death, immune system processes, and several others (Figure 3D). Of note, the gene PR15A was observed to be up-regulated and is associated with the enriched pathways of anatomical structure size (GO:0090066) and regulation of angiogenesis (GO:0045765). Up-regulation of the gene CATIN was also observed and is associated with multiple enriched immune functions, such as negative regulation of toll-like receptor signaling pathway (GO:0034122), negative regulation of tumor necrosis factor production (GO:0032720) and negative regulation of innate immune response (GO:0045824). ALKB5 was up-regulated and is associated with the enriched pathway response to hypoxia (GO:0001666). A broad set of down-regulated genes that were observed were associated with enriched functional categories such as response to stress (GO:0006950) and immune system processes (GO:0002376; Table 2).

Discussion

Amphibian tadpoles respond to perceived predation risk through developmental plasticity, and stress physiology mediates this response (Denver, 2021). Our experiments in the Lithobates – Aeshnidae system reinforce this relationship and also demonstrate that such developmental plasticity corresponds to functionally informative differential gene expression across tissues and tadpole development. PPR triggers the up-regulation of genes associated with developmental pathways, which may play a role in deeper tail growth, connections between the HPI stress axis and tail development, and other effects of chronic stress in a developmentally plastic vertebrate undergoing metamorphosis. By revealing significant differences in transcript abundance and tissue-specific functional pathways during exposure to a stressor, we provide an important foundation for discerning transcriptome-wide pathways of gene expression and relevant gene products, in support of a deeper understanding of how stress-mediated gene regulation may broadly influence developmental plasticity.

Our morphological assessment confirmed that PPR generates an inducible defense that is reliant on HPI axis functionality during tadpole pre-metamorphosis, where the window of developmental plasticity is strongest (Denver, 2009). During pro-metamorphosis, differences between the tail depth of PPR-exposed tadpoles and control tadpoles were minimal. This is not unexpected, as pro-metamorphic tadpoles often reach a large enough size they are no longer subject to predation by ambush predators (Jara and Perotti, 2010), and thus the development of deeper tails into late-stage metamorphosis does not incur a fitness benefit. Furthermore, similar changes in tail shape arising in pro-metamorphic control tadpoles, which has previously been observed in the Lithobates – Aeshnidae system (Hossie and Murray, 2012), implies development of deeper tail morphology may be an underlying metamorphic process that can be accelerated to occur in pre-metamorphosis via stress-induced gene expression (Denver, 2021). During pre-metamorphosis, where developmental plasticity is evident, up-regulation of c-Fos in pre-metamorphic tail tissues is associated with pathways likely responsible for tail growth, including skeletal muscle cell differentiation (GO:0035914) and nervous system development (GO:0007399). C-Fos is commonly induced by growth factors (Ong et al., 1987) and can play a role in angiogenesis (Marconcini et al., 1999) and myogenesis (Kami et al., 1995). Similarly, myelin expression factor 2 was up-regulated and is associated with pathways implying pre-metamorphic tail growth, including development of striated muscle cells (GO:0051146) and neuron differentiation (GO:0030182). More broadly, additional indications of morphological or developmental processes in head and tail tissue included pathway enrichment for programmed cell death and regulation of apoptotic processes (GO:0043068, GO:0043065, GO:0012501). Programmed cell death is an important morphogenic process for amphibians during metamorphosis (Estabel et al., 2003; Nakajima et al., 2005) and can be mediated by hormonal regulation (Tata, 1994). Broad differential regulation of apoptotic genes may, therefore, correspond to chronic stimulation of the stress axis in organisms that are exposed to perceived risk. Finally, up-regulation of protein Phosphatase 1 Regulatory Subunit 15A (also referred to as Gadd34) in tail tissue of pro-metamorphic tadpoles is associated with enriched pathways that regulate anatomical structure size (GO:0090066). Gadd34 can be up-regulated following stress-induced growth arrest conditions (Fornace et al., 1989) and this gene family broadly suppresses cell growth and stimulates apoptosis (Zhan et al., 1994). Up-regulation of cellular growth inhibition, associated with enriched pathways regulating anatomical structure size, may indicate developmental consequences of chronic stress exposure, as evidenced by reduced body mass post-metamorphosis in amphibians exposed to long-term chronic stress (Hu et al., 2008).

During exposure to PPR, the HPI stress axis, in particular corticosterone (CORT) production, mediates developmental plasticity, producing inducible defense morphologies in tadpoles (Denver, 2021). Accordingly, we predicted that differential gene expression would be informative regarding connections between the stress axis and observed morphological plasticity, especially in pre-metamorphic tissues. While our findings did not reveal the complete regulatory pathways connecting exactly how, for example, CORT may interact with glucocorticoid receptors (GRs) to activate gene expression at relevant tissues stimulating tail growth, as a crucial starting point, we did reveal genes which may be implicated in these pathways. For example, PPR caused up-regulation of the transcription factor JunB in pre-metamorphic tail tissues, and this gene is associated with the enriched pathways of response to corticosterone (GO:0051412) and progesterone (GO:0032570). JunB is a potential indicator of the presence of CORT because it can increase following prolonged CORT treatment (Hansson and Fuxe, 2008), which is likely an indication of crosstalk between the AP-1 transcription factor (of which JunB is a subunit), and GRs (Karin and Chang, 2001; Kassel and Herrlich, 2007). Furthermore, previous studies have indicated JunB may be a prominent regulator of tail development in amphibians (Yoshida et al., 2016) as well as playing a role in cell proliferation during tail regeneration (Nakamura et al., 2020). Due to the connections between GRs and the AP-1 transcription factor, it is likely that aspects of AP-1’s induction of tissue re-modeling genes are at play in stress axis mediated developmental plasticity (Karin et al., 1997; Yamashita and McCauley, 2006) and these yet unexplored pathways are promising candidates for connecting the HPI axis to inducible defense morphologies in the tail.

Thyroglobulin (Tg) was also differentially expressed in the head tissues of pro-metamorphic tadpoles, and as an important precursor, likely indicates thyroid hormone (TH) production and thyroid gland development. Given the importance of THs for mediating amphibian metamorphosis (Paul et al., 2022) and the interactive effects of the HPI and HPT axes, it is not surprising to observe Tg differentially regulated by prolonged exposure to PPR (Paul et al., 2022). Specifically, we observed lower Tg transcript levels in head tissue from pro-metamorphic tadpoles exposed to PPR, suggesting that TH levels should also be lower. In pro-metamorphic tadpoles, levels of endogenous TH should increase steadily (Regard et al., 1978; Galton, 1992) and only decrease after the point of metaphoric climax (Thambirajah et al., 2019). Given known inhibitory effects of CORT on the HPT axis (Gutiérrez-Mariscal et al., 2012; Castillo-Campos et al., 2021), it is possible that the down-regulation of Tg is a result of stress eliciting a premature decline in TH production or slower up-regulation of TH. This may be a product of feedback inhibition, as may arise from prolonged stress exposure by CORT on the hypothalamus (Kulkarni and Buchholz, 2014; Kim et al., 2019). This would elicit a decrease in CRF, a peptide that stimulates the HPT axis in tadpoles (Denver, 2021) and thereby leads to TH release. Alternatively, it is also possible that synergistic interactions between TH and CORT influenced tadpole physiology during stress (Kulkarni and Buchholz, 2012, 2014; Shewade et al., 2020; Sterner et al., 2020; Sterner and Buchholz, 2022), especially during pro-metamorphosis. For example, the overall timing and tissue-specific effects of TH on metamorphosis can be regulated by CORT (Sachs and Buchholz, 2019; Paul et al., 2022), through one mechanism of altering tissue sensitivity to TH (Galton, 1990; Bonett et al., 2010). Therefore, the lower abundance of Tg producing transcripts in pro-metamorphosis may also reflect an alteration in developmental timing, such as an accelerated decline in endogenous TH mediated by the stress response in PPR exposed tadpoles. Further potential indications of altered metamorphic processes occurring in PPR-exposed tadpoles include up-regulation of SFTPA in pro-metamorphic head tissue, an important respiratory protein for lung development (Chang et al., 2022), and up-regulation of ALKB5, an indicator of response to hypoxia, which often occurs in amphibian tissues during metamorphosis (Hastings and Burggren, 1995).

PPR exposure elicited differential gene expression in a variety of functional pathways unrelated to inducible defenses. For example, immune system genes were both up and down-regulated across all tissue types and metamorphic stages. This observed differential regulation may be a product of the vertebrate neuroendocrine system’s multiple interactive connections to immune system processes (Dunn, 2007). Furthermore, chronic stress exposure can have immune impacts in amphibians (Gervasi and Foufopoulos, 2008; Falso et al., 2015), and predation risk is a known stressor and immune system trigger in other vertebrate systems (Meuthen et al., 2020; Roncalli et al., 2020). In the case of the amphibian immune system, dynamic re-modeling occurs during pro-metamorphosis, resulting in immune impairment (Robert and Ohta, 2009). Our observed differential regulation of immune system processes, especially during pro-metamorphosis, may indicate heightened amphibian immune system susceptibility as a result of chronic stress exposure.

We conclude that this study is an important first step in understanding transcriptome-wide responses to stress in developing amphibians and highlights how observed phenotypic plasticity can be accompanied by changes in gene expression related to developmental plasticity. The inducible defenses of anuran tadpoles are among the most impressive examples of developmental and phenotypic plasticity triggered by environmental stressors. The tadpole model system is not only useful for studying complex ecological interactions, but also as a convenient method to assess the sensitivity of neuroendocrine pathways to environmental stressors during vertebrate development. Future studies may build on our foundation of transcription-wide analysis to focus on more specific tissues or endocrine level responses. This approach should be complemented by improved data resolution. For example, a high quality assembled and annotated Ranidae genome could support better gene annotation and thereby provide a more comprehensive understanding of complete regulatory pathways connecting the neuroendocrine response to developmental plasticity. Furthermore, we sought to limit the bias of handling stress by measuring transcriptomic responses separately from morphological responses, but future studies should seek to measure these responses together, to better elicit the effect of individual variation on developmental plasticity. The results we obtained from de novo assembly highlight the utility of this study system, and that this group of animals are well deserving of more extensive genomic resources. This is not only an investment for the increasing ecological understanding and conservation of these declining taxa (Green et al., 2020), but also support for future research that is relevant as a model system for developmental plasticity and the neuroendocrinology of stress.

References

• 1

  AdamsD. C.Otárola-CastilloE. (2013). geomorph: an r package for the collection and analysis of geometric morphometric shape data. Methods Ecol. Evol.4, 393–399. doi: 10.1111/2041-210X.12035

  • CrossRef
  • Google Scholar

• 2

  AndrewsS. (2010). FastQC A Quality Control tool for High Throughput Sequence Data. Available online at: https://www.bioinformatics.babraham.ac.uk/projects/fastqc/ (Accessed February 21, 2022).

  • Google Scholar

• 3

  BaniahmadA.TsaiM. J. (1993). Mechanisms of transcriptional activation by steroid hormone receptors. J. Cell. Biochem.51, 151–156. doi: 10.1002/jcb.240510206

  • CrossRef
  • Google Scholar

• 4

  BeckermanA. P.WieskiK.BairdD. J. (2007). Behavioural versus physiological mediation of life history under predation risk. Oecologia152, 335–343. doi: 10.1007/s00442-006-0642-6

  • CrossRef
  • Google Scholar

• 5

  BeldadeP.MateusA. R. A.KellerR. A. (2011). Evolution and molecular mechanisms of adaptive developmental plasticity. Mol. Ecol.20, 1347–1363. doi: 10.1111/j.1365-294X.2011.05016.x

  • CrossRef
  • Google Scholar

• 6

  BenjaminiY.HochbergY. (1995). Controlling the false discovery rate: A practical and powerful approach to multiple testing. J. R. Stat. Soc Ser. B Methodol.57, 289–300. doi: 10.1111/j.2517-6161.1995.tb02031.x

  • CrossRef
  • Google Scholar

• 7

  BennettA. M.LonghiJ. N.ChinE. H.BurnessG.KerrL. R.MurrayD. L. (2016). Acute changes in whole body corticosterone in response to perceived predation risk: A mechanism for anti-predator behavior in anurans? Gen. Comp. Endocrinol.229, 62–66. doi: 10.1016/j.ygcen.2016.02.024

  • CrossRef
  • Google Scholar

• 8

  BennettA. M.PereiraD.MurrayD. L. (2013). Investment into Defensive Traits by Anuran Prey (Lithobates pipiens) Is Mediated by the Starvation-Predation Risk Trade-Off. PloS One8, e82344. doi: 10.1371/journal.pone.0082344

  • CrossRef
  • Google Scholar

• 9

  BoeckmannB.BairochA.ApweilerR.BlatterM.-C.EstreicherA.GasteigerE.et al. (2003). The SWISS-PROT protein knowledgebase and its supplement TrEMBL in 2003. Nucleic Acids Res.31, 365–370. doi: 10.1093/nar/gkg095

  • CrossRef
  • Google Scholar

• 10

  BolgerA. M.LohseM.UsadelB. (2014). Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics30, 2114–2120. doi: 10.1093/bioinformatics/btu170

  • CrossRef
  • Google Scholar

• 11

  BonettR. M.HoopferE. D.DenverR. J. (2010). Molecular mechanisms of corticosteroid synergy with thyroid hormone during tadpole metamorphosis. Gen. Comp. Endocrinol.168, 209–219. doi: 10.1016/j.ygcen.2010.03.014

  • CrossRef
  • Google Scholar

• 12

  BreuschT. S.PaganA. R. (1979). A simple test for heteroscedasticity and random coefficient variation. Econometrica47, 1287–1294. doi: 10.2307/1911963

  • CrossRef
  • Google Scholar

• 13

  BryantD. M.JohnsonK.DiTommasoT.TickleT.CougerM. B.Payzin-DogruD.et al. (2017). A tissue-mapped axolotl de novo transcriptome enables identification of limb regeneration factors. Cell Rep.18, 762–776. doi: 10.1016/j.celrep.2016.12.063

  • CrossRef
  • Google Scholar

• 14

  BuskirkJ.RelyeaR. A. (1998). Selection for phenotypic plasticity in Rana sylvatica tadpoles. Biol. J. Linn. Soc65, 301–328. doi: 10.1111/j.1095-8312.1998.tb01144.x

  • CrossRef
  • Google Scholar

• 15

  Castillo-CamposA.Gutiérrez-MataA.CharliJ.-L.Joseph-BravoP. (2021). Chronic stress inhibits hypothalamus–pituitary–thyroid axis and brown adipose tissue responses to acute cold exposure in male rats. J. Endocrinol. Invest.44, 713–723. doi: 10.1007/s40618-020-01328-z

  • CrossRef
  • Google Scholar

• 16

  ChangL.ZhangM.ChenQ.LiuJ.ZhuW.JiangJ. (2022). From Water to Land: The Structural Construction and Molecular Switches in Lungs during Metamorphosis of Microhyla fissipes. Biology11, 528. doi: 10.3390/biology11040528

  • CrossRef
  • Google Scholar

• 17

  ClinchyM.SheriffM. J.ZanetteL. Y. (2013). Predator-induced stress and the ecology of fear. Funct. Ecol.27, 56–65. doi: 10.1111/1365-2435.12007

  • CrossRef
  • Google Scholar

• 18

  CollinsK. S.GazleyM. F. (2017). Does my posterior look big in this? The effect of photographic distortion on morphometric analyses. Paleobiology43, 508–520. doi: 10.1017/pab.2016.48

  • CrossRef
  • Google Scholar

• 19

  CreelS.ChristiansonD. (2008). Relationships between direct predation and risk effects. Trends Ecol. Evol.23, 194–201. doi: 10.1016/j.tree.2007.12.004

  • CrossRef
  • Google Scholar

• 20

  DaytonG. H.SaenzD.BaumK. A.LangerhansR. B.DeWittT. J.LindströmJ. (2005). Body shape, burst speed and escape behavior of larval anurans. Oikos111, 582–591. doi: 10.1111/j.1600-0706.2005.14340

  • CrossRef
  • Google Scholar

• 21

  DenverR. J. (2009). Stress hormones mediate environment-genotype interactions during amphibian development. Gen. Comp. Endocrinol.164, 20–31. doi: 10.1016/j.ygcen.2009.04.016

  • CrossRef
  • Google Scholar

• 22

  DenverR. (2021). Stress hormones mediate developmental plasticity in vertebrates with complex life cycles. Neurobiol. Stress14, 100301. doi: 10.1016/j.ynstr.2021.100301

  • CrossRef
  • Google Scholar

• 23

  DeWittT. J.SihA.WilsonD. S. (1998). Costs and limits of phenotypic plasticity. Trends Ecol. Evol.13, 77–81. doi: 10.1016/S0169-5347(97)01274-3

  • CrossRef
  • Google Scholar

• 24

  DunnA. J. (2007). “The HPA axis and the immune system: A perspective,” in NeuroImmune Biology (Elsevier), 3–15s. doi: 10.1016/S1567-7443(07)00201-3

  • CrossRef
  • Google Scholar

• 25

  EstabelJ.MercerA.KönigN.ExbrayatJ.-M. (2003). Programmed cell death in Xenopus laevis spinal cord, tail and other tissues, prior to, and during, metamorphosis. Life Sci.73, 3297–3306. doi: 10.1016/j.lfs.2003.06.015

  • CrossRef
  • Google Scholar

• 26

  FalsoP. G.NobleC. A.DiazJ. M.HayesT. B. (2015). The effect of long-term corticosterone treatment on blood cell differentials and function in laboratory and wild-caught amphibian models. Gen. Comp. Endocrinol.212, 73–83. doi: 10.1016/j.ygcen.2015.01.003

  • CrossRef
  • Google Scholar

• 27

  Ferland-RaymondB.MurrayD. L. (2008). Predator diet and prey adaptive responses: Can tadpoles distinguish between predators feeding on congeneric vs. conspecific prey? Can. J. Zool.86, 1329–1336. doi: 10.1139/Z08-117

  • CrossRef
  • Google Scholar

• 28

  FinnR. D.ClementsJ.EddyS. R. (2011). HMMER web server: interactive sequence similarity searching. Nucleic Acids Res.39, W29–W37. doi: 10.1093/nar/gkr367

  • CrossRef
  • Google Scholar

• 29

  FornaceA. J.NebertD. W.HollanderM. C.LuethyJ. D.PapathanasiouM.FargnoliJ.et al. (1989). Mammalian genes coordinately regulated by growth arrest signals and DNA-damaging agents. Mol. Cell. Biol.9, 4196–4203. doi: 10.1128/mcb.9.10.4196-4203.1989

  • CrossRef
  • Google Scholar

• 30

  FrakerM. E.HuF.CuddapahV.McCollumS. A.RelyeaR. A.HempelJ.et al. (2009). Characterization of an alarm pheromone secreted by amphibian tadpoles that induces behavioral inhibition and suppression of the neuroendocrine stress axis. Horm. Behav.55, 520–529. doi: 10.1016/j.yhbeh.2009.01.007

  • CrossRef
  • Google Scholar

• 31

  FrakerM. E.LudsinS. A.LuttbegB.DenverR. J. (2021). Stress hormone-mediated antipredator morphology improves escape performance in amphibian tadpoles. Sci. Rep.11, 4427. doi: 10.1038/s41598-021-84052-9

  • CrossRef
  • Google Scholar

• 32

  GaffneyK. P.PrammerM.BrasfieldL.HippD. R.KennedyD.PatelJ. M. (2022). SQLite: past, present, and future. Proc. VLDB Endow.15, 3535–3547. doi: 10.14778/3554821.3554842

  • CrossRef
  • Google Scholar

• 33

  GaltonV. A. (1990). Mechanisms underlying the acceleration of thyroid hormone-induced tadpole metamorphosis by corticosterone. Endocrinology127, 2997–3002. doi: 10.1210/endo-127-6-2997

  • CrossRef
  • Google Scholar

• 34

  GaltonV. A. (1992). The role of thyroid hormone in amphibian metamorphosis. Trends Endocrinol. Metab.3, 96–100. doi: 10.1016/1043-2760(92)90020-2

  • CrossRef
  • Google Scholar

• 35

  GervasiS. S.FoufopoulosJ. (2008). Costs of plasticity: responses to desiccation decrease post-metamorphic immune function in a pond-breeding amphibian. Funct. Ecol.22, 100–108. doi: 10.1111/j.1365-2435.2007.01340.x

  • CrossRef
  • Google Scholar

• 36

  GilbertS. F. (2001). Ecological developmental biology: developmental biology meets the real world. Dev. Biol.233, 1–12. doi: 10.1006/dbio.2001.0210

  • CrossRef
  • Google Scholar

• 37

  GlennemeierK. A.DenverR. J. (2002a). Developmental changes in interrenal responsiveness in anuran amphibians1. Integr. Comp. Biol.42, 565–573. doi: 10.1093/icb/42.3.565

  • CrossRef
  • Google Scholar

• 38

  GlennemeierK. A.DenverR. J. (2002b). Role for corticoids in mediating the response of Rana pipiens tadpoles to intraspecific competition. J. Exp. Zool.292, 32–40. doi: 10.1002/jez.1140

  • CrossRef
  • Google Scholar

• 39

  GlennemeierK. A.DenverR. J. (2002c). Small changes in whole-body corticosterone content affect larval Rana pipiens fitness components. Gen. Comp. Endocrinol.127, 16–25. doi: 10.1016/S0016-6480(02)00015-1

  • CrossRef
  • Google Scholar

• 40

  GrabherrM. G.HaasB. J.YassourM.LevinJ. Z.ThompsonD. A.AmitI.et al. (2011). Trinity: reconstructing a full-length transcriptome without a genome from RNA-Seq data. Nat. Biotechnol.29, 644–652. doi: 10.1038/nbt.1883

  • CrossRef
  • Google Scholar

• 41

  GreenD. M.LannooM. J.LesbarrèresD.MuthsE. (2020). Amphibian population declines: 30 years of progress in confronting a complex problem. Herpetologica76, 97–100. doi: 10.1655/0018-0831-76.2.97

  • CrossRef
  • Google Scholar

• 42

  GuL.QinS.SunY.HuangJ.AkbarS.ZhangL.et al. (2022). Coping with antagonistic predation risks: Predator-dependent unique responses are dominant in Ceriodaphnia cornuta. Mol. Ecol.31, 3951–3962. doi: 10.1111/mec.16550

  • CrossRef
  • Google Scholar

• 43

  Gutiérrez-MariscalM.SánchezE.García-VázquezA.Rebolledo-SolleiroD.CharliJ.-L.Joseph-BravoP. (2012). Acute response of hypophysiotropic thyrotropin releasing hormone neurons and thyrotropin release to behavioral paradigms producing varying intensities of stress and physical activity. Regul. Pept.179, 61–70. doi: 10.1016/j.regpep.2012.08.010

  • CrossRef
  • Google Scholar

• 44

  HaasB. J. (2018). TransDecoder Release v5.5.0. Available online at: https://github.com/TransDecoder/TransDecoder/releases/tag/TransDecoder-v5.5.0 (Accessed December 2, 2024).

  • Google Scholar

• 45

  HanssonA. C.FuxeK. (2008). Time-course of immediate early gene expression in hippocampal subregions of adrenalectomized rats after acute corticosterone challenge. Brain Res.1215, 1–10. doi: 10.1016/j.brainres.2008.03.080

  • CrossRef
  • Google Scholar

• 46

  HarvellC. D. (1990). The ecology and evolution of inducible defenses. Q. Rev. Biol.65, 323–340. doi: 10.1086/416841

  • CrossRef
  • Google Scholar

• 47

  HastingsD.BurggrenW. (1995). Developmental changes in oxygen consumption regulation in larvae of the South African clawed frog Xenopus Laevis. J. Exp. Biol.198, 2465–2475. doi: 10.1242/jeb.198.12.2465

  • CrossRef
  • Google Scholar

• 48

  HellstenU.HarlandR. M.GilchristM. J.HendrixD.JurkaJ.KapitonovV.et al. (2010). The genome of the western clawed frog xenopus tropicalis. Science328, 633–636. doi: 10.1126/science.1183670

  • CrossRef
  • Google Scholar

• 49

  HettyeyA.TóthZ.ThonhauserK. E.FrommenJ. G.PennD. J.Van BuskirkJ. (2015). The relative importance of prey-borne and predator-borne chemical cues for inducible antipredator responses in tadpoles. Oecologia179, 699–710. doi: 10.1007/s00442-015-3382-7

  • CrossRef
  • Google Scholar

• 50

  HossieT. J.Ferland-RaymondB.BurnessG.MurrayD. L. (2010). Morphological and behavioural responses of frog tadpoles to perceived predation risk: A possible role for corticosterone mediation? Écoscience17, 100–108. doi: 10.2980/17-1-3312

  • CrossRef
  • Google Scholar

• 51

  HossieT.LandoltK.MurrayD. L. (2017). Determinants and co-expression of anti-predator responses in amphibian tadpoles: a meta-analysis. Oikos126. doi: 10.1111/oik.03305

  • CrossRef
  • Google Scholar

• 52

  HossieT. J.MurrayD. L. (2012). Assessing behavioural and morphological responses of frog tadpoles to temporal variability in predation risk. J. Zool.288, 275–282. doi: 10.1111/j.1469-7998.2012.00955.x

  • CrossRef
  • Google Scholar

• 53

  HuF.CrespiE. J.DenverR. J. (2008). Programming neuroendocrine stress axis activity by exposure to glucocorticoids during postembryonic development of the frog, Xenopus laevis. Endocrinology149, 5470–5481. doi: 10.1210/en.2008-0767

  • CrossRef
  • Google Scholar

• 54

  JaraF. G.PerottiM. G. (2010). Risk of predation and behavioural response in three anuran species: influence of tadpole size and predator type. Hydrobiologia644, 313–324. doi: 10.1007/s10750-010-0196-9

  • CrossRef
  • Google Scholar

• 55

  JohnsonM.ZaretskayaI.RaytselisY.MerezhukY.McGinnisS.MaddenT. L. (2008). NCBI BLAST: a better web interface. Nucleic Acids Res.36, W5–W9. doi: 10.1093/nar/gkn201

  • CrossRef
  • Google Scholar

• 56

  KamiK.NoguchiK.SenbaE. (1995). Localization of myogenin, c-fos, c-jun, and muscle-specific gene mRNAs in regenerating rat skeletal muscle. Cell Tissue Res.280, 11–19. doi: 10.1007/BF00304506

  • CrossRef
  • Google Scholar

• 57

  KarinM.ChangL. (2001). AP-1–glucocorticoid receptor crosstalk taken to a higher level. J. Endocrinol.169, 447–451. doi: 10.1677/joe.0.1690447

  • CrossRef
  • Google Scholar

• 58

  KarinM.LiuZ.ZandiE. (1997). AP-1 function and regulation. Curr. Opin. Cell Biol.9, 240–246. doi: 10.1016/S0955-0674(97)80068-3

  • CrossRef
  • Google Scholar

• 59

  KasselO.HerrlichP. (2007). Crosstalk between the glucocorticoid receptor and other transcription factors: molecular aspects. Mol. Cell. Endocrinol.275, 13–29. doi: 10.1016/j.mce.2007.07.003

  • CrossRef
  • Google Scholar

• 60

  KaufmanD.SweetR. (1974). Contrast coding in least squares regression analysis. Am. Educ. Res. J.11, 359–377. doi: 10.2307/1162790

  • CrossRef
  • Google Scholar

• 61

  KenwardM. G.RogerJ. H. (1997). Small sample inference for fixed effects from restricted maximum likelihood. Biometrics53, 983–997. doi: 10.2307/2533558

  • CrossRef
  • Google Scholar

• 62

  KikuyamaS.KawamuraK.TanakaS.YamamotoK. (1993). Aspects of amphibian metamorphosis: hormonal control. Int. Rev. Cytol.145, 105–148. doi: 10.1016/s0074-7696(08)60426-x

  • CrossRef
  • Google Scholar

• 63

  KimJ. S.HanS. Y.IremongerK. J. (2019). Stress experience and hormone feedback tune distinct components of hypothalamic CRH neuron activity. Nat. Commun.10, 5696. doi: 10.1038/s41467-019-13639-8

  • CrossRef
  • Google Scholar

• 64

  KishidaO.NishimuraK. (2006). Flexible architecture of inducible morphological plasticity. J. Anim. Ecol.75, 705–712. doi: 10.1111/j.1365-2656.2006.01091.x

  • CrossRef
  • Google Scholar

• 65

  KishidaO.TrussellG. C.NishimuraK. (2007). Geographic variation in a predator-induced defense and its genetic basis. Ecology88, 1948–1954. doi: 10.1890/07-0132.1

  • CrossRef
  • Google Scholar

• 66

  KroghA.LarssonB.von HeijneG.SonnhammerE. L. (2001). Predicting transmembrane protein topology with a hidden Markov model: application to complete genomes. J. Mol. Biol.305, 567–580. doi: 10.1006/jmbi.2000.4315

  • CrossRef
  • Google Scholar

• 67

  KulkarniS. S.BuchholzD. R. (2012). Beyond synergy: corticosterone and thyroid hormone have numerous interaction effects on gene regulation in Xenopus tropicalis tadpoles. Endocrinology153, 5309–5324. doi: 10.1210/en.2012-1432

  • CrossRef
  • Google Scholar

• 68

  KulkarniS. S.BuchholzD. R. (2014). Corticosteroid signaling in frog metamorphosis. Gen. Comp. Endocrinol.203, 225–231. doi: 10.1016/j.ygcen.2014.03.036

  • CrossRef
  • Google Scholar

• 69

  LangmeadB.SalzbergS. L. (2012). Fast gapped-read alignment with Bowtie 2. Nat. Methods9, 357–359. doi: 10.1038/nmeth.1923

  • CrossRef
  • Google Scholar

• 70

  LentE. M.BabbittK. J. (2020). The effects of hydroperiod and predator density on growth, development, and morphology of wood frogs (Rana sylvatica). Aquat. Ecol.54, 369–386. doi: 10.1007/s10452-020-09748-y

  • CrossRef
  • Google Scholar

• 71

  LiB.DeweyC. N. (2011). RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinf.12, 323. doi: 10.1186/1471-2105-12-323

  • CrossRef
  • Google Scholar

• 72

  LiedtkeH. C.GarridoJ. G.Esteve-CodinaA.GutM.AliotoT.Gomez-MestreI. (2019). De novo assembly and annotation of the larval transcriptome of two spadefoot toads widely divergent in developmental rate. G3 GenesGenomesGenetics9, 2647–2655. doi: 10.1534/g3.119.400389

  • CrossRef
  • Google Scholar

• 73

  LimaS. L.DillL. M. (1990). Behavioral decisions made under the risk of predation: a review and prospectus. Can. J. Zool.68, 619–640. doi: 10.1139/z90-092

  • CrossRef
  • Google Scholar

• 74

  LindJ.CresswellW. (2005). Determining the fitness consequences of antipredation behavior. Behav. Ecol.16, 945–956. doi: 10.1093/beheco/ari075

  • CrossRef
  • Google Scholar

• 75

  LoveM. I.HuberW.AndersS. (2014). Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol.15, 550. doi: 10.1186/s13059-014-0550-8

  • CrossRef
  • Google Scholar

• 76

  MarconciniL.MarchiòS.MorbidelliL.CartocciE.AlbiniA.ZicheM.et al. (1999). c-fos-induced growth factor/vascular endothelial growth factor D induces angiogenesis in vivo and in vitro. Proc. Natl. Acad. Sci.96, 9671–9676. doi: 10.1073/pnas.96.17.9671

  • CrossRef
  • Google Scholar

• 77

  MateoJ. M. (2010). Alarm calls elicit predator-specific physiological responses. Biol. Lett.6, 623–625. doi: 10.1098/rsbl.2010.0118

  • CrossRef
  • Google Scholar

• 78

  McCollumS. A.LeimbergerJ. D. (1997). Predator-induced morphological changes in an amphibian: predation by dragonflies affects tadpole shape and color. Oecologia109, 615–621. doi: 10.1007/s004420050124

  • CrossRef
  • Google Scholar

• 79

  McCollumS. A.Van BuskirkJ. (1996). Costs and benefits of a predator-induced polyphenism in the gray treefrog hyla chrysoscelis. Evolution50, 583–593. doi: 10.1111/j.1558-5646.1996.tb03870.x

  • CrossRef
  • Google Scholar

• 80

  MelottoA.FicetolaG. F.PennatiR.AnconaN.ManentiR. (2021). Raised by aliens: constant exposure to an invasive predator triggers morphological but not behavioural plasticity in a threatened species tadpoles. Biol. Invasions23, 3777–3793. doi: 10.1007/s10530-021-02603-7

  • CrossRef
  • Google Scholar

• 81

  MeuthenD.MeuthenI.BakkerT. C. M.ThünkenT. (2020). Anticipatory plastic response of the cellular immune system in the face of future injury: chronic high perceived predation risk induces lymphocytosis in a cichlid fish. Oecologia194, 597–607. doi: 10.1007/s00442-020-04781-y

  • CrossRef
  • Google Scholar

• 82

  Middlemis MaherJ.WernerE. E.DenverR. J. (2013). Stress hormones mediate predator-induced phenotypic plasticity in amphibian tadpoles. Proc. R. Soc B Biol. Sci.280, 20123075. doi: 10.1098/rspb.2012.3075

  • CrossRef
  • Google Scholar

• 83

  MiyakawaH.ImaiM.SugimotoN.IshikawaY.IshikawaA.IshigakiH.et al. (2010). Gene up-regulation in response to predator kairomones in the water flea, Daphnia pulex. BMC Dev. Biol.10, 45. doi: 10.1186/1471-213X-10-45

  • CrossRef
  • Google Scholar

• 84

  MommerB. C.BellA. M. (2014). Maternal experience with predation risk influences genome-wide embryonic gene expression in threespined sticklebacks (Gasterosteus aculeatus). PloS One9, e98564. doi: 10.1371/journal.pone.0098564

  • CrossRef
  • Google Scholar

• 85

  MoriT.HirakaI.KurataY.KawachiH.KishidaO.NishimuraK. (2005). Genetic basis of phenotypic plasticity for predator-induced morphological defenses in anuran tadpole, Rana pirica, using cDNA subtraction and microarray analysis. Biochem. Biophys. Res. Commun.330, 1138–1145. doi: 10.1016/j.bbrc.2005.03.091

  • CrossRef
  • Google Scholar

• 86

  MoriT.KawachiH.ImaiC.SugiyamaM.KurataY.KishidaO.et al. (2009). Identification of a novel uromodulin-like gene related to predator-induced bulgy morph in anuran tadpoles by functional microarray analysis. PloS One4, e5936. doi: 10.1371/journal.pone.0005936

  • CrossRef
  • Google Scholar

• 87

  MoriT.YanagisawaY.KitaniY.SugiyamaM.KishidaO.NishimuraK. (2015). Gene expression profiles in Rana pirica tadpoles following exposure to a predation threat. BMC Genomics16, 258. doi: 10.1186/s12864-015-1389-4

  • CrossRef
  • Google Scholar

• 88

  MutzK.-O.HeilkenbrinkerA.LönneM.WalterJ.-G.StahlF. (2013). Transcriptome analysis using next-generation sequencing. Curr. Opin. Biotechnol.24, 22–30. doi: 10.1016/j.copbio.2012.09.004

  • CrossRef
  • Google Scholar

• 89

  NakajimaK.FujimotoK.YaoitaY. (2005). Programmed cell death during amphibian metamorphosis. Semin. Cell Dev. Biol.16, 271–280. doi: 10.1016/j.semcdb.2004.12.006

  • CrossRef
  • Google Scholar

• 90

  NakamuraM.YoshidaH.TakahashiE.WlizlaM.Takebayashi-SuzukiK.HorbM. E.et al. (2020). The AP-1 transcription factor JunB functions in Xenopus tail regeneration by positively regulating cell proliferation. Biochem. Biophys. Res. Commun.522, 990–995. doi: 10.1016/j.bbrc.2019.11.060

  • CrossRef
  • Google Scholar

• 91

  OngJ.YamashitaS.MelmedS. (1987). Insulin-like growth factor I induces c-fos messenger ribonucleic acid in L6 rat skeletal muscle cells*. Endocrinology120, 353–357. doi: 10.1210/endo-120-1-353

  • CrossRef
  • Google Scholar

• 92

  PaulB.SternerZ. R.BuchholzD. R.ShiY.-B.SachsL. M. (2022). Thyroid and corticosteroid signaling in amphibian metamorphosis. Cells11, 1595. doi: 10.3390/cells11101595

  • CrossRef
  • Google Scholar

• 93

  PetersenT. N.BrunakS.von HeijneG.NielsenH. (2011). SignalP 4.0: discriminating signal peptides from transmembrane regions. Nat. Methods8, 785–786. doi: 10.1038/nmeth.1701

  • CrossRef
  • Google Scholar

• 94

  PetrusekA.TollrianR.SchwenkK.HaasA.LaforschC. (2009). A “crown of thorns” is an inducible defense that protects Daphnia against an ancient predator. Proc. Natl. Acad. Sci.106, 2248–2252. doi: 10.1073/pnas.0808075106

  • CrossRef
  • Google Scholar

• 95

  PetterssonL. B.BrönmarkC. (1997). Density-dependent costs of an inducible morphological defense in crucian carp. Ecology78, 1805–1815. doi: 10.1890/0012-9658(1997)078[1805:DDCOAI]2.0.CO;2

  • CrossRef
  • Google Scholar

• 96

  PhugeS.TapkirS.BhandV.KourG.PanditR. (2020). Comparative analysis of anti-predator behaviour and life history traits of the tadpoles exposed to predation risk and corticosterone. Proc. Zool. Soc73, 220–226. doi: 10.1007/s12595-019-00320-7

  • CrossRef
  • Google Scholar

• 97

  PreisserE. L.BolnickD. I.BenardM. F. (2005). Scared to death? The effects of intimidation and consumption in predator–prey interactions. Ecology86, 501–509. doi: 10.1890/04-0719

  • CrossRef
  • Google Scholar

• 98

  PriceS. J.GarnerT. W. J.BallouxF.RuisC.PaszkiewiczK. H.MooreK.et al. (2015). A de novo Assembly of the Common Frog (Rana temporaria) Transcriptome and Comparison of Transcription Following Exposure to Ranavirus and Batrachochytrium dendrobatidis. PloS One10, e0130500. doi: 10.1371/journal.pone.0130500

  • CrossRef
  • Google Scholar

• 99

  PruittK. D.TatusovaT.MaglottD. R. (2005). NCBI Reference Sequence (RefSeq): a curated non-redundant sequence database of genomes, transcripts and proteins. Nucleic Acids Res.33, D501–D504. doi: 10.1093/nar/gki025

  • CrossRef
  • Google Scholar

• 100

  RegardE.TaurogA.NakashimaT. (1978). Plasma thyroxine and triiodothyronine levels in spontaneously metamorphosing Rana catesbeiana tadpoles and in adult anuran amphibia. Endocrinology102, 674–684. doi: 10.1210/endo-102-3-674

  • CrossRef
  • Google Scholar

• 101

  RelyeaR. A. (2001). Morphological and behavioral plasticity of larval anurans in response to different predators. Ecology82, 523–540. doi: 10.1890/0012-9658(2001)082[0523:MABPOL]2.0.CO;2

  • CrossRef
  • Google Scholar

• 102

  RelyeaR. A. (2002). Costs of phenotypic plasticity. Am. Nat.159, 272–282. doi: 10.1086/338540

  • CrossRef
  • Google Scholar

• 103

  RelyeaR. A. (2004). Fine-tuned phenotypes: tadpole plasticity under 16 combinations of predators and competitors. Ecology85, 172–179. doi: 10.1890/03-0169

  • CrossRef
  • Google Scholar

• 104

  RobertJ.OhtaY. (2009). Comparative and developmental study of the immune system in Xenopus. Dev. Dyn.238, 1249–1270. doi: 10.1002/dvdy.21891

  • CrossRef
  • Google Scholar

• 105

  RobinsonM. D.McCarthyD. J.SmythG. K. (2010). edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinforma. Oxf. Engl.26, 139–140. doi: 10.1093/bioinformatics/btp616

  • CrossRef
  • Google Scholar

• 106

  RoncalliG.SolerM.TielemanB. I.VersteeghM. A.Ruiz-RayaF.ColomboE.et al. (2020). Immunological changes in nestlings growing under predation risk. J. Avian Biol.51. doi: 10.1111/jav.02271

  • CrossRef
  • Google Scholar

• 107

  RoseboomP. H.NandaS. A.BakshiV. P.TrentaniA.NewmanS. M.KalinN. H. (2007). Predator threat induces behavioral inhibition, pituitary-adrenal activation and changes in amygdala CRF-binding protein gene expression. Psychoneuroendocrinology32, 44–55. doi: 10.1016/j.psyneuen.2006.10.002

  • CrossRef
  • Google Scholar

• 108

  RowJ. R.DonaldsonM. E.LonghiJ. N.SavilleB. J.MurrayD. L. (2016). Tissue-specific transcriptome characterization for developing tadpoles of the northern leopard frog (Lithobates pipiens). Genomics108, 232–240. doi: 10.1016/j.ygeno.2016.10.002

  • CrossRef
  • Google Scholar

• 109

  SachsL. M.BuchholzD. R. (2019). Insufficiency of thyroid hormone in frog metamorphosis and the role of glucocorticoids. Front. Endocrinol.10. doi: 10.3389/fendo.2019.00287

  • CrossRef
  • Google Scholar

• 110

  SanogoY.HankisonS.BandM.ObregonA.BellA. (2011). Brain transcriptomic response of threespine sticklebacks to cues of a predator. Brain. Behav. Evol.77, 270–285. doi: 10.1159/000328221

  • CrossRef
  • Google Scholar

• 111

  SchlickerA.DominguesF. S.RahnenführerJ.LengauerT. (2006). A new measure for functional similarity of gene products based on Gene Ontology. BMC Bioinf.7, 302. doi: 10.1186/1471-2105-7-302

  • CrossRef
  • Google Scholar

• 112

  SchoeppnerN. M.RelyeaR. A. (2009). Interpreting the smells of predation: how alarm cues and kairomones induce different prey defences. Funct. Ecol.23, 1114–1121. doi: 10.1111/j.1365-2435.2009.01578.x

  • CrossRef
  • Google Scholar

• 113

  SchurchN. J.SchofieldP.GierlińskiM.ColeC.SherstnevA.SinghV.et al. (2016). How many biological replicates are needed in an RNA-seq experiment and which differential expression tool should you use? RNA22, 839–851. doi: 10.1261/rna.053959.115

  • CrossRef
  • Google Scholar

• 114

  SchwarzenbergerA.CourtsC.von ElertE. (2009). Target gene approaches: Gene expression in Daphnia magna exposed to predator-borne kairomones or to microcystin-producing and microcystin-free Microcystis aeruginosa. BMC Genomics10, 527. doi: 10.1186/1471-2164-10-527

  • CrossRef
  • Google Scholar

• 115

  SergioC.LucaR.OlivierF. (2021). Plasticity and flexibility in the anti-predator responses of treefrog tadpoles. Behav. Ecol. Sociobiol.75, 142. doi: 10.1007/s00265-021-03078-1

  • CrossRef
  • Google Scholar

• 116

  ShapiroS. S.WilkM. B. (1965). An analysis of variance test for normality (Complete samples). Biometrika52, 591–611. doi: 10.2307/2333709

  • CrossRef
  • Google Scholar

• 117

  ShewadeL. H.SchoephoersterJ. A.PatmannM. D.KulkarniS. S.BuchholzD. R. (2020). Corticosterone is essential for survival through frog metamorphosis. Endocrinology161, bqaa193. doi: 10.1210/endocr/bqaa193

  • CrossRef
  • Google Scholar

• 118

  SteinerU. K. (2007). Investment in defense and cost of predator-induced defense along a resource gradient. Oecologia152, 201–210. doi: 10.1007/s00442-006-0645-3

  • CrossRef
  • Google Scholar

• 119

  SternerZ. R.BuchholzD. R. (2022). Glucocorticoid receptor mediates corticosterone-thyroid hormone synergy essential for metamorphosis in Xenopus tropicalis tadpoles. Gen. Comp. Endocrinol.315, 113942. doi: 10.1016/j.ygcen.2021.113942

  • CrossRef
  • Google Scholar

• 120

  SternerZ. R.ShewadeL. H.MertzK. M.SturgeonS. M.BuchholzD. R. (2020). Glucocorticoid receptor is required for survival through metamorphosis in the frog Xenopus tropicalis. Gen. Comp. Endocrinol.291, 113419. doi: 10.1016/j.ygcen.2020.113419

  • CrossRef
  • Google Scholar

• 121

  SupekF.BošnjakM.ŠkuncaN.ŠmucT. (2011). REVIGO summarizes and visualizes long lists of gene ontology terms. PloS One6, e21800. doi: 10.1371/journal.pone.0021800

  • CrossRef
  • Google Scholar

• 122

  TataJ. R. (1994). Hormonal regulation of programmed cell death during amphibian metamorphosis. Biochem. Cell Biol.72, 581–588. doi: 10.1139/o94-077

  • CrossRef
  • Google Scholar

• 123

  TeplitskyC.PlénetS.JolyP. (2005). Costs and limits of dosage response to predation risk: to what extent can tadpoles invest in anti-predator morphology? Oecologia145, 364–370. doi: 10.1007/s00442-005-0132-2

  • CrossRef
  • Google Scholar

• 124

  ThambirajahA. A.KoideE. M.ImberyJ. J.HelbingC. C. (2019). Contaminant and environmental influences on thyroid hormone action in amphibian metamorphosis. Front. Endocrinol.10. doi: 10.3389/fendo.2019.00276

  • CrossRef
  • Google Scholar

• 125

  TillsO.TruebanoM.FeldmeyerB.PfenningerM.MorgenrothH.SchellT.et al. (2018). Transcriptomic responses to predator kairomones in embryos of the aquatic snail Radix balthica. Ecol. Evol.8, 11071–11082. doi: 10.1002/ece3.4574

  • CrossRef
  • Google Scholar

• 126

  ToddE. V.BlackM. A.GemmellN. J. (2016). The power and promise of RNA-seq in ecology and evolution. Mol. Ecol.25, 1224–1241. doi: 10.1111/mec.13526

  • CrossRef
  • Google Scholar

• 127

  TouchonJ. C.WarkentinK. M. (2008). Fish and dragonfly nymph predators induce opposite shifts in color and morphology of tadpoles. Oikos117, 634–640. doi: 10.1111/j.0030-1299.2008.16354.x

  • CrossRef
  • Google Scholar

• 128

  TylerR. E.WeinbergB. Z. S.LovelockD. F.OrnelasL. C.BesheerJ. (2020). Exposure to the predator odor TMT induces early and late differential gene expression related to stress and excitatory synaptic function throughout the brain in male rats. Genes Brain Behav.19, e12684. doi: 10.1111/gbb.12684

  • CrossRef
  • Google Scholar

• 129

  Van BuskirkJ.AnderwaldP.LüpoldS.ReinhardtL.SchulerH. (2003). The lure effect, tadpole tail shape, and the target of dragonfly strikes. J. Herpetol.37, 420–424. doi: 10.1670/0022-1511(2003)037[0420:TLETTS]2.0.CO;2

  • CrossRef
  • Google Scholar

• 130

  Van BuskirkJ.ArioliM. (2005). Habitat specialization and adaptive phenotypic divergence of anuran populations. J. Evol. Biol.18, 596–608. doi: 10.1111/j.1420-9101.2004.00869.x

  • CrossRef
  • Google Scholar

• 131

  Van BuskirkJ.SchmidtB. R. (2000). Predator-induced phenotypic plasticity in larval newts: trade-offs, selection, and variation in nature. Ecology81, 3009–3028. doi: 10.1890/0012-9658(2000)081[3009:PIPPIL]2.0.CO;2

  • CrossRef
  • Google Scholar

• 132

  VaretH.Brillet-GuéguenL.CoppéeJ.-Y.DilliesM.-A. (2016). SARTools: A DESeq2- and EdgeR-based R pipeline for comprehensive differential analysis of RNA-Seq data. PloS One11, e0157022. doi: 10.1371/journal.pone.0157022

  • CrossRef
  • Google Scholar

• 133

  VerdolinJ. L. (2006). Meta-analysis of foraging and predation risk trade-offs in terrestrial systems. Behav. Ecol. Sociobiol.60, 457–464. doi: 10.1007/s00265-006-0172-6

  • CrossRef
  • Google Scholar

• 134

  VinterstareJ.HulthénK.NilssonP. A.LangerhansR. B.ChauhanP.HanssonB.et al. (2023). Sex matters: predator presence induces sexual dimorphism in a monomorphic prey, from stress genes to morphological defences. Evolution77, 304–317. doi: 10.1093/evolut/qpac030

  • CrossRef
  • Google Scholar

• 135

  VinterstareJ.HulthénK.NilssonP. A.Nilsson SköldH.BrönmarkC. (2020). Experimental manipulation of perceived predation risk and cortisol generates contrasting trait trajectories in plastic crucian carp. J. Exp. Biol.223, jeb213611. doi: 10.1242/jeb.213611

  • CrossRef
  • Google Scholar

• 136

  WaterhouseR. M.SeppeyM.SimãoF. A.ManniM.IoannidisP.KlioutchnikovG.et al. (2018). BUSCO applications from quality assessments to gene prediction and phylogenomics. Mol. Biol. Evol.35, 543–548. doi: 10.1093/molbev/msx319

  • CrossRef
  • Google Scholar

• 137

  WhiteB. A.NicollC. S. (1981). “Hormonal control of amphibian metamorphosis,” in Metamorphosis: A Problem in Developmental Biology. Eds. GilbertL. I.FriedenE. (Springer US, Boston, MA), 363–396. doi: 10.1007/978-1-4613-3246-6_11

  • CrossRef
  • Google Scholar

• 138

  YamashitaJ.McCauleyL. K. (2006). The activating protein-1 transcriptional complex. Clin. Rev. Bone Miner. Metab.4, 107–122. doi: 10.1385/BMM:4:2:107

  • CrossRef
  • Google Scholar

• 139

  YangW.QiY.FuJ. (2016). Genetic signals of high-altitude adaptation in amphibians: a comparative transcriptome analysis. BMC Genet.17, 134. doi: 10.1186/s12863-016-0440-z

  • CrossRef
  • Google Scholar

• 140

  YoshidaH.OkadaM.Takebayashi-SuzukiK.UenoN.SuzukiA. (2016). Involvement of JunB proto-oncogene in tail formation during early xenopus embryogenesis. Zoolog. Sci.33, 282–289. doi: 10.2108/zs150136

  • CrossRef
  • Google Scholar

• 141

  ZhanQ.LordK. A.AlamoI.HollanderM. C.CarrierF.RonD.et al. (1994). The gadd and MyD genes define a novel set of mammalian genes encoding acidic proteins that synergistically suppress cell growth. Mol. Cell. Biol.14, 2361–2371. doi: 10.1128/mcb.14.4.2361-2371.1994

  • CrossRef
  • Google Scholar

• 142

  ZhaoS.LiC.-I.GuoY.ShengQ.ShyrY. (2018). RnaSeqSampleSize: real data based sample size estimation for RNA sequencing. BMC Bioinf.19, 191. doi: 10.1186/s12859-018-2191-5

  • CrossRef
  • Google Scholar