Temporal Profile of Brain Gene Expression After Prey Catching Conditioning in an Anuran Amphibian

A key goal in modern neurobiology is to understand the mechanisms underlying learning and memory. To that end, it is essential to identify the patterns of gene expression and the temporal sequence of molecular events associated with learning and memory processes. It is also important to ascertain if and how these molecular events vary between organisms. In vertebrates, learning and memory processes are characterized by distinct phases of molecular activity involving gene transcription, structural change, and long-term maintenance of such structural change in the nervous system. Utilizing next generation sequencing techniques, we profiled the temporal expression patterns of genes in the brain of the fire-bellied toad Bombina orientalis after prey catching conditioning. The fire-bellied toad is a basal tetrapod whose neural architecture and molecular pathways may help us understand the ancestral state of learning and memory mechanisms in tetrapods. Differential gene expression following conditioning revealed activity in molecular pathways related to immediate early genes (IEG), cytoskeletal modification, axon guidance activity, and apoptotic processes. Conditioning induced early IEG activity coinciding with transcriptional activity and neuron structural modification, followed by axon guidance and cell adhesion activity, and late neuronal pruning. While some of these gene expression patterns are similar to those found in mammals submitted to conditioning, some interesting divergent expression profiles were seen, and differential expression of some well-known learning-related mammalian genes is missing altogether. These results highlight the importance of using a comparative approach in the study of the mechanisms of leaning and memory and provide molecular resources for a novel vertebrate model in the relatively poorly studied Amphibia.


INTRODUCTION
Broadly defined, learning is the change in behavior that comes as a result of experience (Mackintosh, 1974;Rudy, 2014). Information about prior experiences is established and stored as memories in brain neuronal networks and allows for the subsequent modification of behavior. At the level of the neuron, information encoding is mediated by changes in the efficacy of synaptic transmission (Castillo, 2012;Squire et al., 2012). The persistent increase in synaptic efficacy is termed long-term potentiation (LTP), while the weakening of such efficacy is termed long-term depression; both presynaptic and post-synaptic processes mediate these phenomena. Among the molecular mechanisms of long-term memory and learning, our knowledge of LTP is most advanced. LTP involves structural modification of the synapse as well as growth of new synaptic connections. While many mechanisms underlie LTP, these mechanisms share distinct temporal phases: induction is the initial event that triggers LTP, expression is an established change in synaptic transmission, and maintenance is the process of stabilizing the change over time (Lamprecht, 2014;Schaefer et al., 2017). Induction of long-term memory formation is independent of protein synthesis and mediated by the depolymerization and lengthening of actin cytoskeletal elements in axons and dendrites (Lynch et al., 2007;Rudy, 2015a). Induction also involves activity-dependent Ca +2 second messenger signaling cascades which lead to the recruitment of AMPA receptors to the post-synaptic density, as well as activation of the MAPK-CREB pathway which results in learning-related gene transcription (Kandel, 2012). These early memory processes begin the events leading to long-term changes in the strength of synaptic transmission. Learningrelated changes in neuroarchitecture also involve modification of larger neuronal structures by neurite growth or pruning, and possibly the birth and differentiation of new neurons in brain regions that display adult neurogenesis (Cameron and Glover, 2015;Augusto-Oliveira et al., 2019). After induction, expression, and maintenance of the change in synaptic transmission are needed for the consolidation of long-term memory. Expression depends on the establishment of mechanisms that potentiate synaptic efficacy. For example, AMPA receptors recruited to the post-synaptic density during induction are phosphorylated, leading to increased conductance post-synaptically. Additionally, the availability of neurotransmitters is increased pre-synaptically by increasing the number of synaptic vesicles available (Abraham and Williams, 2003). Maintenance of long-term synaptic change is the last phase of LTP, and it involves protein synthesis. Waves of delayed expression of many genes relating to signal transduction (camk2a), cytoskeletal organization (arc) and transcriptional regulation (c-fos and egr-1) have been implicated in the persistence of synaptic changes and memories (Katche et al., 2010;Lisman et al., 2012). Functionally, the maintenance phase involves stabilizing the increase in AMPA receptors at synapses, as well as structural stabilization of synaptic modifications, such as synaptic enlargement (Desmond and Levy, 1986;Stewart et al., 2000;Blitzer, 2005;Rudy, 2015b). While much is known about the molecular pathways underlying vertebrate LTP (Malenka and Bear, 2004), there is still much to be learned especially in terms of the mechanisms underlying maintenance of stored memory.
Invertebrates have long been popular models to investigate learning and memory mechanisms due to the simplicity and accessibility of their nervous systems and vast behavioral repertoires. Many arthropods and cephalopods are capable of complex forms of learning like contextual, spatial, and concept learning (Mather and Kuba, 2013;Perry et al., 2013;Byrne and Hawkins, 2015). Our understanding of the molecular correlates of learning and memory comes from a large body of literature that started with work on Aplysia californica, Caenorhabditis elegans, and Drosophila melanogaster (Kandel and Schwartz, 1982;Sokolowski, 2001;Kaletta and Hengartner, 2006). Many molecular mechanisms underlying long-term memory formation are conserved between invertebrates and mammals (for reviews see Bailey et al., 1996;Cayre et al., 2002;Kandel et al., 2014). Long-term facilitation, the invertebrate analog of LTP, engages many similar molecular signaling pathways as mammalian LTP, moreover, both invertebrates and mammals show synaptic potentiation and depression as parallel processes underling memory formation and maintenance. However, while invertebrates share some learning mechanisms with vertebrates, robust information about the transition from invertebrate to vertebrate learning mechanisms is lacking in the literature. Invertebrates are incredibly diverse and important differences in central nervous system organization could limit their usefulness as comparative models for vertebrate learning. Furthermore, molecular differences in invertebrate and vertebrate learning mechanisms have not been studied in detail, and while some fundamental similarities exist, important differences have been noted, such as the absence of retrograde signaling pathways in invertebrates and different mechanisms of modulation of synaptic plasticity between groups (Glanzman, 2010;Shomrat et al., 2015). Lastly, it has not been determined if similarities in learning mechanisms represent homology or homoplasy, thus compounding the limited usefulness of invertebrates as comparative models. Such differences may suggest evolutionary divergence in some learning and memory mechanisms. Hence, the study of the molecular mechanisms involved learning could benefit from a broader selection of animals, especially representatives of basal vertebrate lineages, which are understudied. A basal vertebrate comparative model would help us to determine if learning mechanisms have been conserved or lost during the evolution of the complex behaviors demonstrated by amniote vertebrates.
Anuran amphibians, while generally understudied, are capable of several forms of learning (e.g., Schmajuk et al., 1980;Daneri et al., 2007;Mitchell and McCormick, 2013;Ramsay et al., 2013;Liu et al., 2016;Miller et al., 2018). The fire-bellied toad Bombina orientalis is a member of the Bombinatoridae, a group which occupies a basal position in anuran phylogeny (Yu et al., 2007) and follows a typical anuran life history with aquatic larvae and terrestrial adults, unlike the more established anuran models of the Xenopus genus which are derived aquatic specialists. B. orientalis offers a closer comparison with amniotes than phylogenetically more distant models, such as invertebrates or teleost fishes and shows extensive homology in brain structure with mammals (Laberge and Roth, 2007). Their phylogenetic position as a basal tetrapod, and brain homology, may help establish the ancestral tetrapod condition for brain structure and behavioral capacity, allowing us to address questions about the evolution of learning and memory.
Here we investigate the temporal sequence of gene expression after conditioning in the brain of the fire-bellied toad. We used a previously developed prey-catching conditioning task to initiate gene expression and assessed whole-brain gene expression by transcriptomic analysis at different time points following conditioning. The initial steps of short-term memory formation occur immediately after a learning event, and as such, transcriptomic analysis is not amenable to such short-term changes mediating early plasticity. Given that immediate early gene (IEG) expression is consistently associated with learning in vertebrates (Jarvis et al., 1995;Guzowski et al., 2001;Mokin and Keifer, 2005;Lau et al., 2011), and is activated very quickly and transiently, time-points must be considered which capture both IEG expression as well as related downstream genes whose expression depends on IEGs. Many IEGs are known to show expression as early as 30 min following a learning event and can persist for up to 4-6 h. Furthermore, there is evidence that structural changes of the neuron are still taking place for long periods after a learning event in mammals and late phase IEG activity appears to be a critical component for long-term memory persistence, for instance, c-fos shows late expression linked to the persistence of fear memory (Bekinschtein et al., 2008;Katche et al., 2010;Hertler et al., 2016). Therefore, early (2 h), intermediate (4 h), and late timepoints were included (24 h) in our study. If the underlying molecular pathways are conserved from anurans to mammals, we expect to see similar changes in transcription over time following conditioning. This study serves as a foundational step for investigating the molecular correlates of learning in B. orientalis and provides a new vertebrate model for the comparative study of the molecular mechanisms of learning and memory.

Animals
To quantify gene expression following conditioning, we trained 14 fire-bellied toads on a prey catching task in which the modification of a snapping response toward a visual prey stimulus by the administration of food reward was previously shown to involve learning (Ramsay et al., 2013). Animals were purchased from National Reptile Supply (Mississauga, ON) and housed in glass terrariums of 37 × 22 × 25 cm in dimension. Temperature was held constant at 21 • C, and a 12/12 h lightdark photoperiod was maintained with light onset at 7 h. Terrariums were furnished with a rock-gravel substrate, broken clay pot pieces and large flat rocks for cover. The substrate was kept moist and the animals had constant access to water in a bowl. Prior to behavioral assays, toads were fed juvenile crickets (Acheta domesticus) dusted with vitamin and calcium supplements weekly, ad libitum. Individual toads were identified based on dorsal patterns of coloration. All procedures were approved by the University of Guelph animal care committee (animal utilization protocol #3590) under guidance from the Canadian Council on Animal Care.

Testing Arena
Behavioral assays were performed in a Bussey-Saksida rodent touch screen chamber running the ABET II software (Lafayette Instrument, Lafayette, IN, United States) (Figure 1). Each chamber wall and the touch screen were covered with matte vinyl FIGURE 1 | Training area for prey-catching conditioning. Toads were required to touch the prey stimulus five times, at which point the video was turned off and the toad fed a cricket (Photo by V. Lewis).
coverings adhering to the interior surfaces to prevent the toads from reacting to their own reflections. At the end of every testing session, the Plexiglas floor and screen surfaces were cleaned with a 5% Tergazyme TM solution and rinsed with distilled water to remove mucus build-up left by the toads. The visual stimulus used to stimulate toad prey catching was an 8 × 3 cm video of moving crickets against a white background previously used for similar conditioning experiments (Ramsay et al., 2013), located in the center-bottom portion of the monitor flush with the floor of the enclosure (Figure 1)

Conditioning
Prior to conditioning, each toad was gently handled for 1 min once per day over 5 days, followed by one shaping trial, and 2 days of rest. Shaping consisted of placing a toad in the touch screen chamber for 1 min, with hand feeding of one cricket after 30 s. No video stimulus was playing during shaping trials. All toads successfully completed shaping (consuming the cricket presented) on their first attempt. Two days after the shaping trial, toads underwent prey catching conditioning for a period of about 1.5 h. Conditioning involved two sessions of six trials. The trials within each session were separated by an interval of 3 min, while the two sessions were separated by 30 min. A preliminary experiment showed that this conditioning protocol was sufficient to see a significant change in performance between the two sessions (Supplementary Figure S2). In each trial, toads were required to snap five times at the cricket video stimulus before manual administration of a cricket reward by the experimenter. During each conditioning trial, a toad was placed in the arena facing the video stimulus from a distance of 15 cm and given 2 min to complete the task. The time it took toads to execute the five snaps (i.e., latency to task completion) was measured. Toads were placed temporarily in a separate container between trials and sessions. To insure high prey catching activity, toads were not fed for 7 days before shaping and the size of cricket rewards administered during conditioning was kept small to minimize satiation.
Following completion of the two conditioning sessions, only toads whose performance improved by a pre-defined amount were selected for genomic analysis. This performance criterion was based on the average change in performance (i.e., decreased latency to task completion) between sessions 1 and 2 (Figure 2A). Toads which showed less than 20% decrease in average latency between sessions were not used for genomic analysis. Two toads out of 14 were rejected based on this criterion.

Time-Course
We chose to assess brain gene expression at 2, 4, and 24 h post-conditioning to evaluate gene expression associated with the induction, expression and maintenance phases of long-term memory. The temporal pattern of gene expression following conditioning was studied by sacrificing groups of toads (n = 3 per group) at the three time points post-training. Sacrifice timepoints were based on time elapsed from the start of the first conditioning trial in the first training session. Each toad was returned to its home terrarium during the interval between the end of conditioning and sacrifice. Control toads (n = 3) were housed with the conditioned toads, and handled daily, but never experienced the touch-screen chamber and were sacrificed immediately after being removed from the terrarium. We faced important constraints in our choice of controls. First, we could not use unpaired controls within the time period of 2 h before sampling of the first group of toads because Ramsay et al. (2013) showed that unpairing the prey catching response and cricket administration by 5 min was insufficient to prevent an improvement in prey caching performance with training. In fact, longer delays of 1 h and 1 day had to be used to demonstrate that learning was involved in modification of the innate prey catching response. Second, because prey catching toward a moving prey stimulus is an innate response, a control procedure where toads exposed to the cricket stimulus would not be rewarded could produce aversive learning and thus not be a good control. Third, similar to exposure to the cricket stimulus on its own, controls that would simply be put into the experimental context without the video stimulus playing and without cricket administration could also produce aversive learning because we have observed regularly that fire-bellied toads put into a restricted space try to escape if nothing is there to capture their attention (Laberge, personal observations). In the end, we chose to use unexposed controls as the best method to assess gene expression induced by the prey catching conditioning procedure. This method allows only the assessment of gene candidates that could be associated with learning. We verified the likelihood that the expression of these gene candidates was involved in learning by comparison with other studies (see section Discussion).

Brain Dissection
Before sacrifice, toads were anesthetized by immersion in a 0.05% solution of buffered tricaine methanesulfonate (Argent Chemical Laboratories, Redmond, WA, United States) until loss of reaction to a leg pinch (about 5 min). Thereafter, the head was cut, and the brain was quickly dissected (∼10 min) from a ventral approach in a physiological Ringer's solution consisting of Na + 129 mM, K + 4 mM, Ca 2+ 2.4 mM, Mg 2+ 1.4 mM, Cl − 115 mM, HCO 3− 25 mM, and glucose 10 mM. Cranial nerves, the pituitary gland and the dura mater were removed before transfer of the brain in 700 µl of TRIzol TM (Life Technologies Inc., Carlsbad, CA, United States) followed by freezing on dry ice. The whole brain was used from the olfactory bulbs to the caudal hindbrain, excluding the spinal cord. Therefore, the pathways responsible for control of prey catching behavior (e.g., midbrain: Ewert, 1967;medulla: Takei et al., 1987;Matesz et al., 2014) and their potential modification by learning or other processes (e.g., Carr et al., 2002) were all included in the analysis of molecular responses.

RNA Extraction and Sequencing
To isolate brain tissue RNA, brains were thawed and ground up using a pestle in a 1.5 ml Eppendorf tube before a modified TRIzol TM extraction protocol was used (Supplementary Material, Section 2). Extracted RNA was cleaned up using the RNeasy mini kit (Qiagen, Venlo, Netherlands) before samples were stored at −80 o C. RNA concentration and quality were respectively assessed using the NanoDrop1500 (Thermo Fisher Scientific) and the Agilent 2100 Bioanalyzer (Agilent Technologies; RIN ≥ 1.9; Supplementary Figure S4 and Supplementary Table S8).
Library preparation and sequencing was done by The Hospital for Sick Children's Center for Applied Genomics (Toronto, ON, Canada). Libraries were poly-a filtered and 125 bp strand specific, paired-end sequenced reads were produced at a depth of 25 million reads with a HiSeq2500 system (Illumina Inc., CA, United States). All raw RNA-seq data has been uploaded to the Gene Expression Omnibus (GSE135054).

Quality Control
Removal of indexing base pairs and trimming of low-quality segments was performed using the FastQC/Trimmomatic tools (Babraham Bioinformatics, Cambridgeshire, United Kingdom). Approximately 250 million reads across the 12 samples were trimmed for quality and length resulting in ∼230 million clean reads (Supplementary Table S1).

De novo Assembly-Alignment -Annotation
De novo assembly of the B. orientalis brain transcriptome was performed using custom scripts within the Trinity assembly suite, which reconstructs transcripts by partitioning RNA-seq data into multiple De Bruijn graphs (Haas et al., 2013). Read alignment was performed using the Bowtie2 − > eXpress pipeline, and differential gene expression was analyzed using DESeq2 (Langmead and Salzberg, 2012;Roberts and Pachter, 2013;Love et al., 2014). Transcript abundance was estimated in a genomefree manner using the alignment-based eXpress tool. eXpress aligns raw reads back against the de novo transcriptome and outputs expression-count estimates in Transcripts Per Million (TPM). DESeq2 uses these count estimates to analyze differential gene expression. Analysis was done in a pair-wise fashion, with each post-training time point being compared against the control. FIGURE 2 | Prey catching conditioning results. (A) Percent decrease in latency to five snaps from session one to session two for each toad. The dashed line represents the rejection threshold, which was set to 1 standard deviation below the mean of two previous 2-session acquisition experiments (n = 20; 37% ± 15%; Lewis, unpublished). Dark bars are two rejected toads that were not used for transcriptomic analysis. Only toads showing a ≥20% decrease in latency between sessions 1 and 2 were used for transcriptomic analysis. (B) Mean prey-catching performance over 2 sessions of 6 trials for trained toads included in the study. Error bars represent ± Standard Deviation. The decrease in latency between sessions is statistically significant (Paired t-test: t = 9.07, df = 8, p < 0.01).
Differentially expressed genes (DEG) were sorted by a strict ≥ 1 log 2 fold-change and the P-values reported by DESeq2 were used to calculate False Discovery Rates (FDR) for all transcripts using the Benjamini-Hochberg transformation (0.1 ≥ FDR, Benjamini and Hochberg, 1995). Annotation of the de novo transcriptome was performed using the Trinotate annotation suite (Haas et al., 2013), which is part of the Trinity platform. Trinotate takes the de novo assembly, uses the program Transdecoder (Haas et al., 2013) to estimate the longest open reading frames, and performs multiple BLAST homology searches against local databases (NCBI_nr, Swissprot, KEGG, GO).
Analysis of possible taxonomically-restricted anuran gene orthologs in the B. orientalis transcriptome was accomplished by homology search against model Xenopus spp. genomes. Two current Xenopus genomes (X. laevis, v2; X. tropicalis, v9.1) were combined to create a local database and a homology search was performed with the B. orientalis transcriptome using NCBI's BLASTx tool (Altschul et al., 1990;Camacho et al., 2009). The resulting orthologs were then subjected to a protein motif survey to determine coding vs. non-coding transcripts. NCBI's ORFinder website was used to find all open reading frames (ORF), then a motif search was performed for every ORF with an amino acid length greater or equal to 100. The motif survey was performed using the MyHits online portal (Pagni et al., 2007).

Hierarchical Clustering and Enrichment Analyses
Hierarchical clustering of co-expressed genes was performed using the Java-based desktop application Multi-Experiment Viewer (MeV; Chu et al., 2008) using K-means clustering with test parameters set to default (Pearson Correlation distance measure, Absolute distance used, and a Maximum Iteration of 50). Unsorted DEG in the dataset (122981 transcripts; Gene Omnibus: GSE135054) were sorted into 11 co-expression clusters. Gene Ontology enrichment analysis was performed on three clusters (2, 4, and 24 h absolute expression peaks; Figure 3) using the DAVID Bioinformatics Resources website and the PANTHER classification system (Thomas, 2003;Huang et al., 2009), with each cluster individually compared to a background gene list consisting of all gene annotations (n = 19,290 nonredundant annotated transcripts). The ratio of functional group genes to all genes in the set was compared between the individual clusters and the background, and functional groups were considered enriched if the number of genes in a cluster was significantly different compared to the background (Fisher's exact test p ≥ 0.05). Finally, while Gene Ontology enrichment analysis provides an estimate of what functions may be affected by the treatment, they tell us little about possible molecular pathways that may be involved. In order to investigate the pathways involved, we searched the Kyoto Encyclopedia of Genes and Genomes (KEGG) database using the DAVID Bioinformatics resource and performed a similar enrichment analysis.

Behavior
The 12 toads included in the study significantly improved prey catching performance between the first and second training sessions as demonstrated by a reduction in mean latency to complete 5 snaps between the training sessions ( Figure 2B; paired t-test: t = 9.07, df = 8, p < 0.01; Supplementary Figure S3).

Transcriptome Details
RNA-sequencing resulted in 268,925,312 raw reads across the 12 samples and yielded a total of 259,530,365 clean reads (96.5% of raw reads) following quality control procedures (Supplementary Table S1; Gene Omnibus: GSE135054). Transcriptome assembly resulted in 1,236,332 total transcripts consisting of 656,704,220 BP and ultimately 714,091 non-redundant "Trinity" genes (non-isoform transcripts; Supplementary Table S2). Replicate variability was tested using a correlation matrix comparing TPM values across treatments and clearly showed greater between-than within-treatment variability (Supplementary Figure S1). Alignment of raw reads back to the transcriptome resulted in an average of 23% uniquely mapped reads, with an average of 93% overall alignment quality across all samples (Supplementary Table S2).
As there is currently no publicly available genome for B. orientalis, 714,091 aligned transcripts were annotated using sequence Blast homology searches against the NCBI_nr and Swissprot protein databases, resulting in 137,305 transcripts annotated to 24,830 non-redundant genes (19.2% of transcriptome). All BLAST searches reported here were based on a cut-off of E ≥ 0.005 with a minimum Percent Similarity of 30%.
Differential expression was analyzed across three time points by comparing gene expression against the controls. A total of 122,981 contigs were found to be differentially expressed before sorting for fold change and FDR. Ultimately, a total of 4,692 transcripts were found to be significantly differentially expressed overall (p ≤ 0.1 FDR, | Fold change| ≥ 2; Figure 3). Of these, 1,344 annotated DEG were of use (557 genes with Uncharacterized Protein accessions were ignored for subsequent enrichment analysis). All subsequent expression profile and enrichment analyses were performed on these 1,344 annotated DEG.
Finally, a homology search of the 557 uncharacterized DEG against the Xenopus spp. database found 39 orthologs. Of these 39, 22 orthologs were upregulated and 17 were downregulated, with the majority (34) showing expression at a single timepoint (2, 4, or 24 h; Supplementary Table S6). Finally, 23 out of the 39 uncharacterized orthologs contained a protein motif, 15 of which showed a strong motif match (Sigrist, 2002), suggesting that most of these taxonomically-restricted transcripts were effectively coding for proteins (Supplementary Table S6).

Identification of Distinct Temporal Co-expression Patterns
To analyze gene expression patterns over time (both up and down regulation of genes), we adapted the time points of our experiment based on a priori expectations informed by the literature on molecular responses associated with learning. These time points are: Early (2 h), Intermediate (4 h), Late (24 h) in comparison to the beginning of the experiment (0 h). We then performed an unbiased clustering analysis on the entire unsorted DEG dataset (122981 DEG), which resulted in 11 coexpression patterns based on absolute expression levels. Of the 11 co-expression patterns, the majority of DEG (>90%) were found in the Early, Intermediate, and Late co-expression peak profiles (Figure 3). Hence, all subsequent analyses were performed on these three co-expression profiles.

Enrichment Analysis and Annotation of Co-expression Groups
High level gene ontology annotation (GO) was performed on the three co-expression profiles identified above. GO annotations are separated into three broad categories: biological processes (BP), cellular component (CC), and molecular function (MF). All genes were sorted into first level GO sub-groups (herein referred to as functional groups) within these categories (Figure 4). In the biological process category, 989 genes were assigned to eight functional groups with the largest number of genes within the cellular and metabolic processes groups (Figure 4A), indicating that at the highest GO functional group level most DEG are involved in processes related to cell metabolism, growth, and maintenance. The cellular process functional group was unpacked to give a better idea of the specific processes involved. Interestingly, most genes in this group fall into the cytoskeleton organization (GO:0007010), and gene expression (GO:0010467) functional groups by GO level 4 ( Figure 4D). In the cellular component category, 620 genes were assigned to eight functional groups with cell part and organelles showing the greatest gene assignment ( Figure 4B). Lastly, in the molecular function category, 441 genes were assigned to six functional groups, with catalytic activity and binding being most represented ( Figure 4C).
Following GO annotation, a functional group enrichment analysis was performed. Significant GO groups consisting of 2 or less genes were omitted. In total, the Early co-expression profile consisted of 10 BP, 16 CC, and 6 MF enriched functional groups, the Intermediate profile consisted of 18 BP, 8 CC, and 3 MF enriched groups, and the Late profile consisted of 11 BP, 7 CC, and 12 MF enriched groups (Figure 5 and Supplementary Tables S3A-C).
At the Early time point most upregulated genes are related to transcriptional regulation, followed by structural plasticityrelated genes (actin filament and cytoskeleton organization; axon guidance; positive regulation of synapse assembly), and intracellular transport genes (transport, protein transport; Supplementary Table S3A). This activity is localized mainly to the nucleus, intracellularly, and to neuronal cell bodies, with some activity in the axon and growth cone. The MF group shows some minor protein dimerization activity but is dominated by DNA binding activity. Conversely, other groups related to structural plasticity and genesis of neuronal extensions are downregulated at 2 h (neuron projection development, microtubule cytoskeleton organization).
The Intermediate time point shows up-regulation in several genes related to structural plasticity GO groups (axon guidance; single organismal cell-cell adhesion; homophilic cell adhesion via plasma membrane adhesion molecules; heterophilic cellcell adhesion via plasma membrane cell adhesion molecules; cell-matrix adhesion; Figure 5 and Supplementary Table S3B) as well as endocytosis, apoptosis, and metabolic activity. The MF category shows nucleic acid binding activity, while activity localizes mainly to the mitochondrion, transcription factor complex, and filopodium. The most suppressed biological process group is regulation of transcription followed by protein modification, axonogenesis, and metabolic activity.
Finally, the Late time point also shows most up-regulated genes that are related to transcriptional regulation (regulation of transcription, DNA-templated; transcription, DNA-templated). Interestingly, while showing some structural plasticity gene enrichment (actin filament organization, focal adhesion assembly, cytoskeleton organization), and glucose metabolic gene enrichment (positive regulation of glucose import), the  Table S3C). Conversely, the most suppressed groups include cell cycle, actin filament bundle assembly, and cell-cell adhesion. The late profile only shows molecular function activity in the Metabolic Processes and Binding groups. The binding group involves activity related to peptide binding (metal ion binding), DNA/RNA binding, and actin binding. While the metabolic group includes nuclease activity (endonuclease activity), proteolysis activity (aspartic-type endopeptidase activity) and polymerase activity (RNA-directed DNA polymerase activity). Furthermore, the most suppressed MF group is cadherin binding involved in cell-cell adhesion. Lastly, most enriched activity is localized to the spindle and nuclear matrix.
Enriched functional groups within the Biological Process category were used to group DEG under the following general, high level, subjectively selected functional headings: transcriptional regulation, structural plasticity, transport, apoptosis, and neurite genesis (Supplementary Tables S4A-C). In this way genes not represented in the enrichment analysis could be included and discussed. A total of 239 additional DEG from this study were ultimately included.

Pathway Analysis
We searched the KEGG database using all 1,344 annotated DEG to identify putative biological pathways represented in our transcriptome (Table 1). Interestingly, pathways related to virus infection were enriched (hsa05169: Epstein-Barr, hsa05168: herpes simplex), as well as two structural plasticity pathways (hsa04360: Axon guidance, hsa04514: Cell adhesion molecules), and a protein processing pathway (hsa04141: Protein processing in endoplasmic reticulum).

Temporal Expression Patterns After Conditioning
There are currently over 100 described immediate early genes (IEG), and a subset of those are known to be induced during learning tasks (Clayton, 2000;Minatohara et al., 2016). Multiple learning related IEG's were positively expressed at 2 h postconditioning in B. orientalis, with a few at 4 h but none at 24 h ( Table 2). Members of the Fos family c-fos and fosB were found to peak either early or intermediately and return to baseline by 24 h. The genes arc/arg3.1, nr4a1, and egr1 were found to peak at 2 h, showed reduced expression by 4 h and returned to baseline by 24 h. The genes jund and ier2 were both differentially expressed only in the 2 h treatment. Interestingly, fosb and nr4a3, which are known to be expressed early in response to acute stress or other environmental cues (30 min to 2 h; Perrotti et al., 2004;Hawk et al., 2012) were only differentially expressed at 4 h (73.5 and 21.1-fold-change) in our transcriptome. Further, bdnf, homer1, and slc2a3, which are known IEG's associated with learning (Lu et al., 2008;Zhao et al., 2010;Chowdhury and Hell, 2018), were not differentially expressed even though they are present in the transcriptome.
Based on the results of the enrichment and KEGG pathway analyses, the temporal expression of a number of important learning-related genes (based on the mammalian literature) was compared to mammalian expression patterns. These genes, including the aforementioned IEG's, were grouped into generalized functional groups based roughly on gene ontology functional groups: IEG's, structural genes, axon guidance genes, and apoptotic genes as well as a general learning-related category ( Table 2). The structural group includes all genes related to cytoskeletal, microtubule, adhesion, actin, and tubulin activity, and consists of 13 DEG. Structural activity peaks at 2 h with 9 of 13 genes showing differential expression. Axon guidance includes any gene related to repulsion/attraction or extension of existing axons and consists of 11 DEG and 1 non-differentially expressed gene. Axon guidance activity is similar between time-points. Apoptosis included any genes associated with programmed cell death, neurite pruning and the intrinsic/extrinsic apoptotic pathways, and consists of 7 DEG, with activity peaking at 24 h for 5 of 7 genes showing differential expression. General learningrelated is a catch all group for important mammalian learning genes which have only a few gene members (e.g., Signaling pathways, Neurotropic factors) or do not fit well into the other groups (e.g., camk2a or ncs-1) and consists of 14 DEG. Given that these genes were selected based on their importance in mammalian learning, it is interesting that most of them are not differentially expressed after prey catching conditioning in B. orientalis.

DISCUSSION
Prey-catching conditioning in B. orientalis engaged distinct temporal patterns of brain gene expression relating to structural and synaptic plasticity, programmed cell death, and transcriptional regulation; patterns which generally agree with functional patterns of gene expression after a learning event in other vertebrates (Bolhuis et al., 2000;Cavallaro et al., 2002;D'Agata and Cavallaro, 2003;Velho  Each gene within a given pathway is shown with its expression fold-change vs. the control. Red is down-regulation and green is up-regulation. and Mello, 2008; Valles et al., 2011;Hertler et al., 2016). In B. orientalis, gene expression following conditioning starts with an increase in IEG expression at 2 h, which coincides with upregulation of gene expression associated with transcriptional regulation and structural plasticity (axon guidance, actin organization, synapse assembly). This is then followed at 4 h post-conditioning by an increase in axon guidance and cell adhesion activity, and a suppression of transcriptional regulation and axonogenesis. The last time point included in our experiment at 24 h is characterized by an increase in transcriptional regulation, apoptotic activity, and the formation of new neuronal projections. Functional group enrichment analysis highlighted structural plasticity, axon guidance, and programmed cell death as important changes in brain tissue induced by conditioning. While some similarities in the components of these pathways exist between mammals and the fire-bellied toad, and suggest phylogenetic conservation of learning pathways in tetrapods, many interesting differences are evident, specifically relating to the absence of well-known mammalian learning-related genes and the temporal patterns of axonogenesis.

Immediate Early Genes
For the most part, the expression patterns of differentially expressed IEGs in response to prey-catching conditioning are consistent with previous literature on activity-dependent molecular responses in vertebrates (Pinaud, 2004;Burmeister and Fernald, 2005;Pinaud et al., 2006;Wada et al., 2006;  Bifunctional apoptosis regulator (bfar) --533.74 The immediate-early gene (IEG) group represents any gene categorized as IEG. The General learning-related group represents genes chosen for their relevance to mammalian learning that do not fit into the other categories. The Structural group represents any gene involved in cytoskeletal modification. The Guidance group represents any gene involved in axon guidance. The Apoptotic group represents any gene involved in programmed cell death. Dashes (-) indicate no significant difference vs. the control. Green numbers are significantly upregulated while red numbers are downregulated. Where overlap in gene function exists, the gene was subjectively placed into the most relevant group. Valles et al., 2011;Minatohara et al., 2016). The expression of c-fos, egr-1, and acr/arg3.1 have been extensively studied and are consistently associated with neural stimulation and activity (Gallo et al., 2018). In rats peak c-fos expression generally ranges from 30 min to 2 h and returns to baseline around 4 h after a learning event (Cullinan et al., 1995;Barros et al., 2015). Similarly, c-fos peaks as soon as 30 min following song learning in zebra finches (Bolhuis et al., 2000). This resembles what we see in B. orientalis, with c-fos expression levels peaking early, and tapering to baseline by 24 h post-training. Similarly, arc/arg3.1 and nr4a1 expression levels, peak at 30-60 min, and lasts up to 4 h (arc/arg3.1: Plath et al., 2006;Bramham et al., 2010;nr4a1: Hawk et al., 2012). In B. orientalis these genes peak at 2 h, show reduced expression at 4 h and return to baseline by 24 h after training. Still, not all observed IEG expression patterns in B. orientalis are comparable to mammalian or bird patterns. The gene fosb is expected to show an early peak like c-fos (Nestler, 2015) but instead showed only high upregulation at 4 h in B. orientalis. Homer1, which is integral to the regulation of synaptic changes induced by LTP (Chowdhury and Hell, 2018), was not differentially expressed at any time-point after conditioning in B. orientalis, and neither was bdnf, a gene coding for an important neurotrophic factor associated with neurogenesis and longterm memory in mammals and song learning in birds (Lu et al., 2008;Wang et al., 2019). Therefore, despite observing important similarities in IEG expression between B. orientalis and other vertebrates, there were notable differences in that some IEGs showed different temporal patterns of expression while others were not engaged by the learning task used in the present study.

Structural Modification
Gene ontology (GO) enrichment analysis highlighted several pathways which change as time progresses after prey catching conditioning, such as cytoskeletal modification, cell guidance, and cell death/genesis. We see gene expression early on (2-4 h) related to structural modification (actin cytoskeletal organization, axon guidance, adhesion, synapse assembly), then later expression of genes which mediate neurite genesis and apoptosis (Supplementary Table S4C). It appears that, like mammals, the establishment of new neural connections is started early following learning and transitions into neurite remodeling mechanisms later. The early phase of memory formation in mammals involves actin cytoskeleton degradation and mobilization of glutamatergic receptors, followed by the consolidation of plasticity into longterm memory by the stabilization of these actin filaments (Abraham and Williams, 2003;Lynch et al., 2007). Interestingly, the early changes in our time course of gene expression revolve around structural modification and cytoskeletal organization, specifically in terms of actin organization. We see early expression of genes like the tropomodulins (tmod2; tmod3), which are involved in actin-capping and could help stabilize actin filaments, and the glutamate ionotropic receptor AMPA type subunit 1 (gria1), which is mobilized early in learning (Keifer and Zheng, 2010;Rao et al., 2014). Furthermore, several important genes involved in synaptic plasticity are upregulated early on: IQ motif and SEC7 domain-containing protein 3 (iqsec3), known to regulate inhibitory synapse formation (Um et al., 2016), neural cell adhesion molecule L1 (l1cam), involved in synaptic plasticity (Patzke et al., 2016), ephrin type-B receptor 2 (ephb2), involved in dendritic spine modification and excitatory synapse formation (Henkemeyer et al., 2003), and leucine-rich repeat-containing protein 24 (lrrc24), associated with synapse formation (Um et al., 2014). Since shortterm memory formation is independent of protein synthesis (actin de/stabilization, AMPA receptor mobilization, etc.) and relies on cell stores of existing proteins, the gene expression observed at the earliest time point may represent turnover of these components.
Understanding the late expression of actin cytoskeletal genes observed in B. orientalis, as highlighted by GO enrichment analysis, is not straightforward. The structural genes expressed at 2-4 h are associated with modification of the neuron or neurites, while at 24 h the structural genes may be involved in neurite genesis or destruction. At 24 h there are several neuron projection development (GO:0031175) genes involved in neurite genesis or destruction that are differentially expressed.

Shared Molecular Pathways Involving Axon Guidance
We see gene expression associated with axon guidance activity at all timepoints, a well-known mechanism underlying learning and memory processes (Kolodkin and Tessier-Lavigne, 2011;Batool et al., 2019). The four main protein families involved in mammalian axon guidance are netrins, slit homolog proteins, ephrins, and semaphorins. Each family has representatives differentially expressed in our experiment.
The semaphorins 7a, 6a, and 3d (sema7a/6a/3d), which act to repulse growth cones or promote axon outgrowth, show differential expression across the treatments with sema7a downregulated at 2-4 h (outgrowth promoting; Jeroen Pasterkamp et al., 2003), sema3d upregulated at 4 h (axon repulsion; Wolman et al., 2004), and sema6a upregulated at 24 h (axon repulsion; Haklai-Topper et al., 2010). The semaphorins act by binding the receptor families plexins, neuropilins and integrins, however, expression of these receptors in B. orientalis is quite different from the expected mammalian pattern (Neufeld et al., 2002;Drabkin et al., 2014;Alto and Terman, 2017;Boyer and Gupton, 2018). We either saw a different receptor type expressed in toads, for example integrin alpha-7 (itga7) is downregulated at 24 h but it is not associated with semaphorin activity in mammals, or we saw no differential expression, as was the case with plexins. Further, while semaphorin 3d is known to bind neuropilin 1 (nrp-1), nrp-1 was downregulated at 24 h in toads while the related nrp-2 was upregulated along with sema3d at 4 h. Other growth cone manipulating genes differentially expressed in our experiment include the ligand netrin-1 (net-1), which is down-regulated at 2 h, and while its receptors are not differentially expressed, actin-binding LIM protein (ablim) is upregulated at 24 and is a downstream constituent of the net-1 axon attraction pathway (Roof et al., 1997;Gitai et al., 2003). Lastly, the slit homolog 2 protein (slit2) ligand acts in axon repulsion and is highly upregulated at 24 h in toads, while its receptor roundabout homolog 2 (robo2) is downregulated at 2 h (Lopez-Bendito et al., 2007).
While axon repulsion seems to be a late stage event in conditioned B. orientalis, ephrin activity complicates this picture. The ephrin A5 ligand (efna5) is upregulated at 4 h and is known to bind to ephrin type-B receptor (ephb2; Himanen et al., 2004), which is highly upregulated at 2 h, while the ephrin type-A receptor 8 (epha8) is downregulated at 2 h. In mammals, Ephrin signaling was initially thought to be only repulsive but is now considered important for synaptic adhesion, and thus synaptogenesis (Lai and Ip, 2009). This is further evidenced by the fact that ephb2 is known to produce an increase in NMDA receptors, enhancing LTP (Nolt et al., 2011). Ephrin activity in mammals has also been shown to be temporally mediated in that axon attraction and outgrowth is seen early, but eventually transitions into axon repulsion or withdrawal later (Liu et al., 2018).
While axon guidance pathways in conditioned B. orientalis appear to share some similarities with mammals, in terms of the genes involved and late axon repulsion, the overall patterns of expression are inconsistent. Axon outgrowth and attraction are early events during learning in mammals, but they appear to be late events in B. orientalis. Ultimately, axon guidance pathways appear to play an important role in tetrapods learning; however, the exact function of these pathways is not clear.

Shared Molecular Pathways Involving Apoptosis
Functional group analysis showed enrichment of gene expression associated with apoptosis in the toad brains at 4 and 24 h post-conditioning. The genes involved are all part of the receptor-mediated "extrinsic" apoptosis pathway, which is well described in mammals (Elmore, 2007). Conversely, the "intrinsic" apoptosis pathway, which is based on the activation of mitochondrial pathways, does not seem to be part of learningassociated apoptotic processes. Interestingly, the apoptotic cysteine-aspartic proteases (caspases) and the B-cell lymphoma 2 (bcl-2) protein family, which are integral to apoptosis, are not differentially expressed across time points in conditioned toads. However, several genes in the apoptotic pathways show changes in expression.
The tumor necrosis factor receptor superfamily member 6 (fas) is a death receptor which mediates apoptosis through caspase 8 (Casp8) activation, or through the C-Jun N-terminal kinase (JNK) signaling pathway. The fas gene is strongly suppressed at 2-4 h post-conditioning in the toad brain, but FAS-associated factor 1 (faf1), which potentiates FAS-induced cell death (Chu et al., 1995), is highly upregulated at 4 h. Furthermore, the fas receptor-activated gene death domain-associated protein 6 (daxx), which mediates JNK pathway activation and thus apoptosis, peaks at 4 h and is then suppressed at 24 h. The downstream effect of fas-daxx-jnk pathway activation is transcription of the fas receptor and its ligand (fas-l; Zhang J. et al., 2000). Perhaps, this suggests that while Fas/Faf1-mediated apoptotic activity begins at 4 h and lasts at least until 24 h, new fas receptors do not need to be mobilized until very late. Interestingly, programmed cell death protein 4 (pdcd4) and protein kinase C delta type (prkcd), which peak at 24 h in the toad brain, are both associated with altered activity of the fas-daxx-jnk signaling pathway and the NF-kB signaling pathway, leading to pro-apoptotic activity (Emoto et al., 1995;Kilpatrick et al., 2002;Brodie and Blumberg, 2003;Lankat-Buttgereit and Göke, 2003).
Finally, the apoptotic gene death-inducer obliterator 1 (dido1) is involved in the regulation of pdcd4 and is also upregulated at 24 h post-training in the toad brain.
While we see evidence of activation of apoptotic pathways following conditioning in B. orientalis, the underlying purpose of this activation is unclear. One potential clue is that the apoptotic genes expressed in this study appear to modulate apoptotic pathways peripherally. Apoptosis is a well-known mechanism that shapes the developing nervous system, however, during adulthood, learning-related apoptosis activity may function as a pruning mechanism rather than a whole-cell death mechanism. This view is supported by findings showing that apoptosis mechanisms and neurite pruning share some molecular pathways (Riccomagno and Kolodkin, 2015;Mukherjee and Williams, 2017). For example, caspases 3 and 9 as well as the apoptosis regulator BAX (bax) are integral to axonal pruning (Simon et al., 2012;Cusack et al., 2013). Caspases are activated very quickly following induction of neurite pruning, but at a low level of expression, which prevents their full-blown apoptotic functions. While we did not see direct evidence of caspase or Bax/Bak activity in the toad brain, either the aforementioned fas apoptosis pathway or modulation of the Bax/Bak pathway could contribute to neurite pruning associated with learning. For example, the bifunctional apoptosis regulator (bfar) gene, which is suppressed at 24 h in toads, is known to interact with bcl-2 and works to suppress the intrinsic bax activity (Zhang H. et al., 2000;Roth et al., 2003), which may suggest propruning activity. Ultimately, while we see evidence that extrinsic apoptotic pathways similar to mammals are being activated, the results suggest that activation of these pathways in the late stages following prey catching conditioning underlie neuronal circuit refinement through neurite pruning rather than replacement of neurons through full blown cell death and birth of new neurons.
Interestingly, there is evidence that axon guidance gene families, such as the semaphorins and ephrins, also have neurite pruning functions. For example, nrp-2 is known to form a complex with a plexin activated by the sema3f ligand to initiate pruning activity (Bagri et al., 2003;Schuldiner and Yaron, 2015). It is interesting then, that we see nrp-2 upregulated with a semaphorin 3 family member at 4 h in this study. A reverse ephrin b3 (efnb3) signaling pathway (Eph − > Ephrin signaling, where the ephrin-B protein acts as a receptor; Xu and Henkemeyer, 2009), has also been implicated in neural pruning activity and while the efnb3 is not differentially expressed in this study, it has been known to bind to the receptor ephb2 (McClelland et al., 2010), which is upregulated at 2 h in the toad brain and may suggest pathway activation. Such early upregulation of ephb2 may be a preparatory step for late pruning activity.

Divergent Learning Pathways
Our study on a basal tetrapod provides further support for the hypothesis that basic learning pathways are conserved in animals. Several learning-related genes known from mammalian literature were expected to be differentially expressed in the late stages after conditioning in B. orientalis. However, while we saw differential expression of many of these genes, the expression patterns of some well-known learning associated genes were inconsistent with mammalian experiments. Genes like monocarboxylate transporter 1 (slc16a1; Takimoto and Hamada, 2014) and metabotropic glutamate receptor 5 (grm5; Jong et al., 2009) are differentially expressed at different times in toads in comparison to mammals, while calcium/calmodulin-dependent protein kinase 2A (camk2a), neuronal calcium sensor-1 (ncs-1), and nerve growth factor (ngf ), which are all important learningrelated genes in mammals (Weiss et al., 2010;Berry et al., 2012;Lisman et al., 2012), were not differentially expressed in this study. These inconsistencies could represent fundamental differences between the molecular mechanisms of learning in amphibians and other vertebrates. However, caution must be taken before concluding that such fundamental differences exist based on the present study. B. orientalis is a novel model in neurobiology and lacks a sequenced genome. Only 19% of expressed sequences and 30% of DEG were successfully annotated using our approach. This limited annotation likely had cascading effects on the pathway analyses presented here.
Poor genomic resources available for this species surely contributed to a low level of annotation, but Kwon (2015) also found thousands of taxonomically restricted genes in Xenopus tropicalis and X. laevis despite better genome resources. Interestingly, Kwon (2015) points out that Xenopus laevis and X. tropicalis genomes contains ∼18,000 taxonomicallyrestricted genes (Hellsten et al., 2010;Session et al., 2016), possibly suggesting an evolutionary divergence between amphibians and amniotes. Given the basal position of B. orientalis, there may be some divergence in the molecular mechanisms of learning compared to other more derived tetrapods; such divergence would be evidenced by the discovery of taxonomically-restricted anuran orthologs in B. orientalis learning-related gene expression. Indeed, a number of putative taxonomically-restricted gene orthologs were differentially expressed after conditioning in the B. orientalis transcriptome. Therefore, the presence of abundant taxonomically-restricted anuran genes could have contributed to the low annotation success of the B. orientalis brain tissue transcriptome. Additionally, differences in learning-related gene expression compared to other vertebrates may reflect divergent evolution of some learning and memory mechanisms in anuran amphibians.

Methodological Limitations
In addition to poor annotation and the possibility of genomic divergence in vertebrates, the behavioral paradigm used to elicit gene expression may limit what molecular pathways are engaged. The many different tasks used to study the molecular correlates of learning in animals might not all engage comparable gene expression as the prey catching conditioning task used here. Furthermore, time point selection may not have adequately covered important periods of gene expression. While we are confident that the 2 h time point represented the peak of IEG expression based on preliminary qPCR data (Lewis, unpublished), the later time points were based on mammalian literature and light/dark cycle constraints. The 12 h time point was omitted from our experiment due to potential confounding effects of nighttime darkness on gene expression in the diurnal fire-bellied toad. Additionally, because conditioning procedures are difficult to implement in amphibians, our control method allowed only the assessment of candidate gene expression that could be associated with learning. Some gene expression could have instead been associated with exposure to the experimental context or prey stimulus independent of the association between prey catching and food reward. Thus, our conclusions about learning-related gene expression had to rely extensively on comparison with studies done in mammals because relevant gene expression studies are scarce in other animals. Despite the above constraints, our results highlight interesting avenues for future investigation and provide useful genomic resources for comparative research in the molecular basis of learning.

CONCLUSION
This study represents the first in-depth transcriptional analysis of the temporal patterns of gene expression following conditioning in an amphibian. Bombina orientalis is a basal anuran which may provide clues about the ancestral state of learning and memory mechanisms in tetrapods and thus help evaluate the evolution of learning in vertebrates. Many similarities exist between the present results in B. orientalis and learningrelated gene expression in mammals, especially with respect to IEG activity and structural plasticity of the neuron, but many interesting differences were also highlighted. It will be important to determine which of the gene candidates outlined in this study are directly associated with learning processes. This could be done by using a different conditioning task that is amenable to unpaired controls, although it could prove difficult to find such a task for work on amphibians. Alternatively, because extended prey catching training produces an increasingly inflexible and automatic response (Ramsay et al., 2013), these different phases of learning could be used to track gene expression as conditioning events progressively lose their potency to influence behavior. Further study of differences in molecular correlates of learning could provide important clues to help explain differences in cognition among vertebrates. Future work on the brain regions and circuits involved in conditioning in amphibians should also help with the latter.

DATA AVAILABILITY STATEMENT
The datasets generated for this study can be found in the NCBI gene omnibus; accession: GSE135054.

ETHICS STATEMENT
The animal study was reviewed and approved by the University of Guelph animal care committee (animal utilization protocol #3590).