Gene expression profiles in Rana pirica tadpoles following exposure to a predation threat

Background Rana pirica tadpoles show morphological changes in response to a predation threat: larvae of the dragonfly Aeshna nigroflava induce heightened tail depth, whereas larval salamander Hynobius retardatus induce a bulgy morphology with heightened tail depth. Although both predators induce similar tail morphologies, it is possible that there are functional differences between these tail morphs. Results Here, we performed a discriminant microarray analysis using Xenopus laevis genome arrays to compare tail tissues of control and predator-exposed tadpoles. We identified 9 genes showing large-scale changes in their expression profile: ELAV-like1, methyltransferase like 7A, dolichyl-phosphate mannosyltransferase, laminin subunit beta-1, gremlin 1, BCL6 corepressor-like 1, and three genes of unknown identity. A further 80 genes showed greater than 5 fold differences in expression after exposure to dragonfly larvae and 81 genes showed altered expression after exposure to larval salamanders. Predation-threat responsive genes were identified by selecting genes that reverted to control levels of expression following removal of the predator. Thirteen genes were induced specifically by dragonfly larvae, nine others were salamander-specific, and sixteen were induced by both. Functional analyses indicated that some of the genes induced by dragonfly larvae caused an increase in laminins necessary for cell adhesion in the extracellular matrix. The higher expression of gremlin 1 and HIF1a genes after exposure to dragonfly larvae indicated an in vivo hypoxic reaction, while down-regulation of syndecan-2 may indicate impairment of angiogenesis. Exposure to larval salamanders caused down-regulation of XCIRP-1, which is known to inhibit expression of adhesion molecules; the tadpoles showed reduced expression of cα(E)-catenin, small muscle protein, dystrophin, and myosin light chain genes. Conclusion The connective tissue of tadpoles exposed to larval salamanders may be looser. The differences in gene expression profiles induced by the two predators suggest that there are functional differences between the altered tail tissues of the two groups of tadpoles. Electronic supplementary material The online version of this article (doi:10.1186/s12864-015-1389-4) contains supplementary material, which is available to authorized users.


Background
Phenotypic plasticity is the ability to produce different phenotypes under different environmental conditions and to respond to changes in environmental conditions [1][2][3]. This phenomenon has long been studied by evolutionary biologists interested in its adaptive significance [4,5]. One particular type of phenotypic plasticity is that displayed as an inducible defense, which is stimulated directly by cues associated with a predation threat [6][7][8][9][10][11][12]. Predator-induced phenotypic plasticity in anuran tadpoles has been extensively studied. Anuran tadpoles exhibit a range of inducible morphological changes, such as heightened tail depth, in the presence of a threat by various types of pond dwelling predator such as dragonfly larvae [13][14][15]. Tadpoles with a heightened tail depth phenotype show higher survival rates when dragonfly larvae are present [16,17]. Rana pirica tadpoles display a unique bulgy morph when exposed to their main predator, larval salamander Hynobius retardatus [10]. The inducible bulgy morphology is believed to be an evolutionary defense against the gape-limited H. retardatus larvae under an intimate predator-prey relationship [18][19][20][21]; the bulgy morph is only induced by a predation threat from larval salamanders and functions to prevent the tadpoles from being swallowed [10,17].
R. pirica presumably evolved phenotypic plasticity in tadpoles as a defense against specific predators: the inducible bulgy body with heightened tail against larval salamanders and the heightened tail morph against the larvae of the dragonfly Aeshna nigroflava. However, dragonfly larvae are top predators in natural ponds [22], and the survival rate of Rana pirica tadpoles with the salamanderinduced bulgy morph is lower than that of tadpoles with the dragonfly-induced heightened tail when exposed to predation by dragonfly larvae. Furthermore, bulgy morph tadpoles have the same survival rate as non-induced tadpoles when placed with dragonfly larvae. Tadpoles with the dragonfly-induced higher tail morphology are less vulnerable to predation by larval salamanders than non-induced tadpoles, indicating that the higher tail phenotype has adaptive advantages compared to other phenotypes under conditions of salamander and dragonfly predation [17]. Moreover, in the presence of dragonfly risk cues, tadpoles with a bulgy morph and heightened tail induced by an earlier exposure to a salamander reduce only the bulgy body but retain the heightened tail [23].
Thus, the tadpoles can regulate their morph according to differences in predator cues; this behavior offers a valuable system for investigating switching mechanisms in adaptive phenotypic plasticity against predators. The evolutionary aspects of these complex phenotypic changes have been addressed in a number of studies over the last decade [20,[24][25][26][27].
Although some studies have identified genetic variation and geographic differentiation in anuran tadpoles with respect to inducible anti-predation defenses, and have shown that these traits are heritable [20,27], we have very little understanding of the genetic mechanisms involved in the morphogenetic alterations associated with these defense traits.
We previously conducted cDNA subtraction and speciesspecific microarray analyses of epithelial tissues from R. pirica tadpoles, both bulgy morph and non-bulgy morph [28], and identified key genes relating to morphogenetic changes [29]. A larger functional species-specific microarray (3 k array) was prepared in the previous study [29] and used in induction-reversion experiments to analyze mRNAs extracted from the facial tissues of tadpoles. These analyses identified a novel uromodulin-like gene, pirica, that was increasingly up-regulated as the period of exposure to larval salamanders lengthened. The pirica protein was found to contain a zona pellucida domain similar to proteins that function to control water permeability; the pirica gene was shown to be expressed in the superficial epidermis of the tadpole skin [29]. We also demonstrated that water retention in the connective tissue and maintenance of a constant osmotic pressure were important factors for bulgy morph formation, supporting the interpretation that predator-induced expression of pirica in the skin causes retention of absorbed water [30]. The immuno-related proteins hyaluronic acid, histone H3 and 14-3-3 zeta were the most abundant constituents of the liquid aspirated from the connective tissue. These findings suggested that formation of the bulgy morph might also require activation of the innate immune system [30].
We concluded from our previous studies that evolution of the inducible bulgy morphology against the gapelimited H. retardatus larvae involved changes to the control of body water dynamics and that some key genes were involved in production of the bulgy body. As mentioned above, the heightened tail morph induced by dragonfly larvae appears to provide greater protection against predators, although Rana pirica tadpoles are able to adopt another morphology in response to an intermediate predator, i.e., larval salamanders. This ability to switch response suggests expression of both a common set of genes against predators and also a predator-specific set of genes. In the present study, we sought to answer two particular questions: first, do the bulgy morph and the heightened tail morph induced by salamander and dragonfly larvae, respectively, differ with respect to gene expression patterns; second, which genes are induced by all predators and which are predator-specific?

Ethics statement
All animal experiments were conducted by trained personnel in accordance with the guidelines of the Animal Care Committee, Nihon University.

Experimental animals and design of the microarray
Eggs of R. pirica and larvae of the dragonfly Aeshna nigroflava and the salamander H. retardatus were collected from a pond in Hokkaido, Japan, and placed in separate 12 liter aquaria. After hatching, R. pirica tadpoles were fed rabbit chow ad libitum. Larval H. retardatus and A. nigroflava were fed small R. pirica tadpoles ad libitum. Water in all aquaria was changed every second day. The experiment was conducted in a laboratory at 18°C, using a natural day/night (about 14/10 hours) regime. Experiments were performed using 4 liter aquaria (29 × 16.5 cm surface area, height 9 cm) each containing 2 liters of aged tap water. Fifty 10-day-old tadpoles of similar sizes (mean body length = 7.90 ± 0.38 mm (standard deviation), n = 48) were randomly chosen from the holding tank and placed in each aquarium. The water in all aquaria was changed every second day throughout the experimental period.
The protocol for experiments involving continuous or short-term exposure of tadpoles to a predator is summarized in Figure 1. Four predator exposure conditions were used, along with a non-exposed control: 20 aquaria were set up for predation-exposed and control tadpoles (we also set up a further 4 aquaria as a backup only for the salamander treatments). Thus, each treatment group had 4 replicate aquaria. The dragonfly and salamander predation experiments were initiated by the introduction of either a large dragonfly larva in a cage or three larval salamanders (about 18 mm length). In two groups, the predators were removed after 4 days and the tadpoles were allowed to revert to a normal phenotype (−Drago and −Salam). At 8 days, we Figure 1 Experimental design showing control, continuous exposure treatment to dragonfly larvae or larval salamanders, and removal of the predation threat to allow the tadpoles to recover. These five treatment groups were used to produce the RNAs used in the microarray analysis with the Xenopus genome array. In the group with continuous exposure to predators (Drago or Salam), the tadpoles were under predation threat for the full 8 days. In the groups exposed for a limited period and then allowed to recover (−Drago or −Salam), tadpoles were initially kept with dragonfly larvae or larval salamanders for 4 days to produce the predator-induced phenotype; after 4 days, the predation threat was removed and the tadpoles were allowed to revert to their normal phenotype for 4 days. The control group of tadpoles was not exposed to a predation threat. Four replicate groups were used for each treatment; tadpoles were sampled on day 8. Tail tissues from the tadpoles were used for RNA extraction for the microarray analysis. The microarray analysis was performed in triplicate. Abbreviations: Cont, control; Drago, 8 day dragonfly exposure; Salam, 8 day salamander exposure; −Drago, 4 day exposure to dragonfly threat and 4 day recovery; −Salam, 4 day exposure to salamander threat and 4 day recovery.
collected tail tissue samples from tadpoles in the aquaria of each experimental group, and performed the microarray analysis using these triplicate samples.
To minimize loss of tadpoles in the salamander treatment groups (Salam), the larval salamanders were replaced daily by others that had been kept in a holding tank with sufficient R. pirica tadpoles to allow easy feeding. The replacement predators were chosen randomly from each holding tank. Every second day, we counted the number of tadpoles in each aquarium.
For the microarray, RNAs were extracted and prepared from 10 tadpoles from each 4 liter aquaria. Therefore, a total of 40 tadpoles were used for each treatment group. Tail tissue was dissected from each tadpole and cut into small pieces; the tissues from all 40 tadpoles of each treatment group were combined and used for RNA extraction. The extracted RNA from each group of 40 tadpoles served as replications for each comparative experiment.

RNA extraction
The tail tissue samples from tadpoles in each treatment group were placed in a single tube containing RNAlater™ (QIAGEN RNA stabilization reagent). The samples were stored at −80°C until RNA extraction. Total RNA was concentrated using an RNeasy Midi Kit® according to the manufacturer's protocol. Total RNA purity was determined by measuring the A 260 /A 280 absorbance ratio and by visualizing the 18S and 28S ribosomal RNA using gel electrophoresis. RNA samples that had an A 260 /A 280 ratio between 1.7 and 2.1 and passed a visual examination of 18S and 26S ribosomal RNA were further processed for the microarray analysis.

Cross-species microarray
Due to the lack of a large-scale microarray from R. pirica, we chose to make use of a cross-species protocol. Microarrays from Xenopus laevis and X. tropicalis are commercially available (Affymetrix). A comparison of e-values from microarrays of the two species after probing with R. pirica mRNAs indicated there were no statistically significant differences (see Additional file 1 and Additional file 2: Figure S1). We made use of the X. laevis microarray for the analyses of predator-induced gene expression in R. pirica reported here.
Samples were prepared for microarray analysis according to Affymetrix GeneChip® protocols. An aliquot containing 5 μg total RNA was converted to double-stranded cDNA with an Affymetrix One-Cycle cDNA Synthesis kit (Affymetrix). Biotin-labeled cRNA was generated from the cDNA samples by in vitro transcription with T7 RNA polymerase using an Affymetrix IVT Labeling kit (Affymetrix). The biotin-labeled cRNAs were fragmented to an average size of 35-200 bases by incubation at 94°C for 35 min in the fragmentation buffer provided in the Affymetrix sample cleanup kit. Each fragmented sample (15 μg) was hybridized for 16 h at 45°C onto individual Xenopus laevis genome arrays, which enable analysis of 15611 transcripts. The gene chips were washed and stained in an Affymetrix Fluidics Station 450 according to the Gen-eChip®Expression Analysis Technical Manual. The chips were scanned with the GeneChip®Scanner 3000 (Affymetrix), and the data was imported into the GeneChip Operating software (GCOS v1.2). The microarray data analysis was carried out using Agilent GeneSpring v11.5, and the data was normalized using the MAS5 summarization algorithm and baseline to median of all samples by baseline option. The data discussed in this publication have been deposited in the NCBI Gene Expression Omnibus (GEO) and are accessible through GEO Series accession number GSE33250.

Analysis of the microarray data
In this study, it was not appropriate to use normal data processing procedures, such as false discovery rate and other strict statistical methods for screening, due to the low signal intensities obtained in the cross-species hybridization [31]. Thus, to identify genes that were influenced by the predation treatments among the 15611 genes on the microarray, we assumed that the response of each gene was independent of that of other genes. It was necessary to repeat the 15611 independent statistical tests, e.g., analysis of variance, to identify sets of genes that gave different population mean responses depending upon the level of the factor. However, if more than 10000 independent statistical tests with a 5% significance level are performed, then the probability that at least one of the tests will be significant is almost 1. Therefore, we performed multiple comparison tests. As 15611 independent statistical tests, with a family-wise error rate of 5% (i.e., the probability that at least one of the statistical tests will be significant is 5% when all the population means are the same), were performed, then the significance level of each test is 1-(1-0.05) 1/15611 = 0.00032857%, i.e. 0.00032857%. The family-wise error rate has also been termed the experiment-wise error rate (see Hsu, 1996 [32]). Usually, the number of tests in a multiple comparison test is up to about 20, or possibly 100 in an extreme case. The number of statistical tests corresponding to the number of genes is over 10000, which is an extremely large number. If the family-wise error rate is 5%, then the significance level for each test of 0.00032857% is too low to permit detection of differences among population means. However, if we use a family-wise error rate of 90% for each of the 15611 tests, then a significance level of 1-(1-0.90) 1/15611 = 0.0001474867 (or 0.01474867%) is obtained. As a compromise, we performed ANOVA with a significance level of 0.01474867% to identify a set of genes from the total of 15611 on the microarray that might be influenced by the predation treatments.
The second step in the data analysis was validation of the selection of the gene set in the first stage of analysis. To achieve this validation, we applied a discriminant analysis to obtain a function or functions that differed among treatment groups. The canonical discriminant functions are calculated by maximizing a fraction, the projection of a linear function of the variance covariance matrix of between groups to that of the variance covariance matrix within the groups, i.e. maximizing the separation of the groups. As another validation step, we applied a hierarchical cluster analysis; we chose single linkage with Euclidean distance, group average with Euclidean distance and Ward's method with squared Euclidean distance. Single linkage is a suitable method to identify outliers and the other methods are suitable to identify like clods. The question of which method should be used depends on the definition of preferred clusters (see [33]). However, if clusters are clearly separated, i.e. within clusters are concentrated on their center, and between clusters are clearly distinct, then whichever method of cluster analysis is used, we obtain the same conclusion; this shows the robustness of the selection of the gene set at the first analytical step.
We compared predator-specific and non-specific gene expression patterns in tadpoles allowed to recover from exposure dragonfly or larval salamanders. In this analysis, we performed data mining to identify changes in expression greater than five-fold compared to the controls. The fold change calculation was performed for each probe set by comparing the mean value for the predator treatments samples to the mean value for the control samples.

Pathway analysis of genes selected by discriminant analysis
Functional information on the nine genes selected by the discriminant analysis (Table 1) was used in an Ingenuity Pathway Analysis (IPA, Qiagen).
cDNA cloning from R. pirica tadpoles and quantitative real-time PCR cDNA fragments from R. pirica tadpoles were selected on the basis of the results of the cross-species microarray. Mixed total RNA (1 μg) purified from the tail parts of tadpoles exposed to different predation threats were used for cDNA synthesis using 10 pmol of anchored T17ADP primer (5′-GAG TCG ACT CGA GAA TTC T 17 -3′) and Superscript™III (Invitrogen). The selected cDNA fragments were produced by PCR using gene specific primers and KOD FX Neo (TOYOBO). The amplification procedure was performed according to the manufacturer's protocol. The gene specific primers were designed using information obtained from homology searches of cDNA databases of X. tropicalis, X. laevis, and two or three other species (NCBI). Each cDNA fragment was cloned into pCRII-TOPO (Invitrogen), and sequencing was carried out using BigDye terminator version 3.1 and an ABI 3100 DNA sequencer (Applied Biosystems). Quantitative real time PCR was performed using Rotor Gene Q MDx 5plex HRM (QIAGEN). In the qPCR, 18S ribosomal RNA obtained from R. pirica was used as the internal control. qPCR was performed using Rotor-Gene SYBR® Green PCR Kit (QIAGEN) as described in the manufacturer's protocol. For each sample, duplicate qPCR amplifications were performed. A 40-cycle program with a hot start was used; the PCR cycles consisted of 5 seconds at 95°C for denaturation, and 10 seconds at 60°C for annealing and elongation. The qPCR primers for ELAV-like1(LC011860) and 18S (LC011910) are shown in Additional file 3: Table S1.

Discriminant analysis in a cross-species microarray
Microarray hybridization was performed as described in Figure 1. The symbols Cont to −Salam in Figure 1 indicate the comparative design of the microarray analysis, which was performed in triplicate. Hybridization conditions were checked using control genes; the chip from control 1 was eliminated from the analysis owing to hybridization failure ( Figure 2a). A1_at., were selected by ANOVA of the microarray data from the 5 treatment groups (the control and four different predator conditions) in the present experiment ( Figure 1); multiple comparison tests were then performed. For the validation of the selection, discriminant analysis of the cross-species microarray was carried out. In order to avoid the singularity of the variance covariance matrix of the within groups, 9 genes were selected. The first 2 discriminant functions give 97% of the total discrimination using all the discriminant functions, which is sufficient to discriminate among the 5 groups.
The canonical discriminant functions (1) and (2) are: x ¼ −6:973−8:135X :5177:1:S1 at þ 7:910X j j :318:1:S1 at   Figure 2b shows the territorial map based on the two canonical discriminant functions. The scattered points in Figure 2b are the projection points from the nine genes by the two canonical discriminant functions and the numbers 1 to 5 indicate the five treatment groups. The seven lines, i.e., y (i,j) , in Figure 2b indicate the boundary lines between the groups i and j. For instance, the boundary line between groups 2 and 5, i.e., continuous exposure to larval salamanders and removal of the dragonfly larvae, is given by y (2,5) The nine genes selected by this discriminant analysis are summarized in Table 1 and also shown in the hierarchical combined tree in Figure 3a. Clusters were created using single linkage method with Euclidean distance, and the clusters indicated clear separation of experimental condition such as control, −Drago and so on. This validates the selection method for the nine genes. Further, the same separation was also obtained using 2 other clustering methods: group average method with Euclidean distances and Ward's method with squared Euclidean distances.
The results indicated that expression profiles of the 9 genes are robust.

Pathway analysis of the nine genes selected by discriminant analysis
The nine genes selected by discriminant analysis were subjected to an IPA pathway analysis using the corresponding human gene. The pathway analysis showed that five of the genes had a protein-protein interaction through ubiquitin C. There was also an mRNA-protein interaction among ELAV-like1, methyltransferase 7A and a hypothetical protein (RAD9-HUS1-RAD1) ( Figure 4). The interaction of these proteins belongs to one of the pathways involved in developmental disorders, hereditary disorders and neurological diseases in humans. The postulated interactions regarding these diseases are summarized in Additional file 4: Table S2.

Screening of specific genes induced by exposure to a predator
Standard methods were used to determine fold changes in expression of genes following exposure to a predator. We set a threshold change of over 5 fold compared to the control and identified 316 and 301 genes that were induced by exposure to dragonfly larvae or larval salamanders, respectively (Figure 5a and b). These data are also described in Additional file 5: Table S3 and Additional file 6: Table  S4. Of these genes, the identities of only 80 and 81 genes, respectively, were known and these were selected as predator-responsive genes. The expression profiles of these genes were compared by a hierarchical clustering analysis using single linkage with Euclidean distance. This analysis indicated that gene #72 (targeting protein for Xklp2) showed an 11.5 fold increase in expression and gene #23 (UPF0534 protein) showed an 11.3 fold decrease in expression compared to control after exposure to a dragonfly larva ( Figure 6). Gene #72 (Frzb-1 protein) showed a 12.7 fold increase in expression and gene #21 (14-3-3 protein) a 10.9 fold decrease compared to control after exposure to a larval salamander (Figure 7).
We also identified predator-responsive genes in the treatment groups in which the predators were removed after 4 days (Figures 6 and 7). The expression profiles of genes showing more than 5 fold changes in expression are shown in Figures 6b and 7b. Selection of genes that are predatorspecific can be achieved using the procedure illustrated in Figures 6c and 7c. Thus, for example, the expression level of a predation-threat responsive gene might return to the control level upon removal of the predation threat. Selection of genes that show this pattern of expression should identify those that are predator-specific. Using this approach, we obtained 13 and 9 genes that were specifically responsive to a predation threat by dragonfly and larval salamanders, respectively. These genes are indicated by the pink and blue boxes in Figures 6a  and 7a, respectively. A further 16 genes were common to both dragonfly larvae and larval salamanders (Figures 6 and 7); these genes are indicated by the green boxes. The results from hierarchical cluster analyses using single linkage with Euclidean distance of the common and predation-threat responsive genes are summarized in Figure 8a.
In these analyses, the continuous predator exposure groups and the removal of predator groups each formed clusters, and the control formed an outer cluster.
Salamander and dragonfly-specific genes are also summarized (Figure 8b and c). Among those that were real salamander-specific were gene #1 (cold inducible RNAbinding protein), gene #3 (protein kinase), gene #22 (clo-neIM28/E3-111), gene #23 (enhancer of split related 10), gene #40 (transcription factor 12), gene #67 (tetratricopeptide repeat protein 30), and gene #69 (ribonucleoprotein: Elav-like family member 3) (Figure 8b). Although these 9 genes showed similar patterns, cold inducible RNA-binding protein, enhancer of split related 10, and tetratricopeptide repeat protein 30 showed opposite expression patterns between Salam and Drago treatment groups. With regard to the dragonfly-specific genes in cDNA cloning and quantitative real-time PCR Eight genes from R. pirica tadpoles showed E-values for their sequence data of less than 10 −72 compared with to X. laevis. These data suggesting that these R. pirica cDNA sequences were similar to those of X. laevis. Furthermore, qPCR analysis of ELAV-like1 likewise showed a similar tendency to that observed in the cross-species microarray using a uniformly most powerful invariant test. The E-values of the sequence data are given in Additional file 7: Table S5 and the statistical analyses of the qPCR are presented in Additional file 8: Figure S2 and Additional file 9.

Nine genes identified in the discriminant analysis
Predator-induced phenotypic plasticity results in morphological adaptations that confer protection against predators. The alterations to the tail in tadpoles in the presence of either the top predator (dragonfly larvae) or the intermediate predator (larval salamanders) may provide protection through enabling an increased swimming speed to escape from the former or by prevention against swallowing by the latter. It is a reasonable presumption that the  Screening of predator-induced genes showing a greater than 5 fold difference compared to control. The selected genes are depicted by volcano plotting, and the threshold change for gene screening was set as 'more than 5 fold change compared to control'. In total, 316 and 301 genes respectively were identified as induced by dragonfly larvae (a) and larval salamanders (b). Fold change is expressed as log 2 X in the X axis.  phenotypic changes in the predation-related effects are mediated by stress responsive genes.
To determine the identities of the stress responsive genes and to characterize their function, we made use of a Xenopus laevis gene chip to screen for functional genes. The use of this cross-species microarray approach was due to the unavailability of a large scale microarray from Rana pirica. Although cross-species microarrays suffer from signal reduction, they have nevertheless been used successfully in comparative evolutionary and ecological studies, and for gene-expression profiling of many species that lack representative microarray platforms [31].
In total, nine genes were identified here as predationresponsive, and their functions have been determined through database searches. The ELAV (embryonic lethal, abnormal vision)-like 1 (HuR) protein has various roles in different cellular mechanisms including apoptosis [34,35], oncogenesis [36], cellular responses to stress [37], and the stabilization of several cellular mRNAs [38]. ELAV affects cell fate by regulating the stability and translation of mRNAs that encode cell stress response proteins [37] and stress-induced cell death proteins [35] and also determine responses to DNA damaging agents such as UVC and actinomycin D [39]. Although the available information indicates that HuR is highly related to stress conditions, a dramatic change in ELAV-like 1 (HuR) gene expression (about a 3.3 fold decrease compared with control) was observed following the removal of the larval salamanders ( Figure 2 and Table 1). Since the gene showed an increased expression under the other experimental conditions, including removal of the dragonfly larvae (Table 1), then it is likely that the gene is not involved in a response to predation stress here.
Methyltransferase-like 7a was down-regulated in all treatments compared with controls ( Figure 2 and Table 1). In particular, a 6.2 fold reduction was observed following removal of the larval salamanders. The gene is a methyltransferase and may function in histone methylation or other metabolic processes. To date, however, there are no detailed descriptions of the precise function of this gene. Interestingly, expression of the gene is altered in posttraumatic stress disorder and shows a significant decrease in human patients [40]. As described earlier, dragonfly larvae are top predators and very efficient killers of tadpoles when they are kept in the same aquarium. Although larval salamanders also predate tadpoles, they are gape-limited and restricted in the size of the prey they can swallow [10,17]. Thus, the difference between the predation threat posed by the dragonfly and salamander is simply one of relative effectiveness. In the present study, some of the tadpoles were bitten by larval salamanders and a few were killed; however, as the dragonfly larvae were caged, then no tadpoles were harmed in these treatment groups. Therefore, in the salamander treatment groups, it is possible that the tadpoles were showing a stress response to the threat of predation.
Laminins are large, heterotrimeric glycoproteins that are major components of both the intestinal extracellular matrix and the basal lamina [41]. These glycoproteins are involved in cell adhesion [42], neurite outgrowth [43], the formation of epithelial structures [44], and are especially important for the formation of an initial polymeric scaffold of cell-attached matrix via interactions with integrins, lutheran, agrin, dystroglycan, and other components. Laminins, type IV collagen and nidogens are thought to constitute the central basement membrane scaffolding [45]. Thus, laminins clearly interact with other proteins to contribute to cell adhesion and, therefore, presumably promote basement membrane anchorage. Interestingly, we found here that laminin subunit beta-1-like was increased 8.4 fold in tadpoles continuously exposed to dragonfly larvae and 4.3 fold in those in which the larvae were removed (Table 1). Moreover, laminin B1 and laminin LII were also increased 5.6 and 7.5 fold compared to the control in the treatment group with continuous exposure to dragonfly larvae ( Figure 6). However, in the treatment group with continuous exposure to larval salamanders, a 2.46 fold increase compared to control was observed; in the group in which the larval salamanders were removed a 1.3 fold decrease was found. It will be of interest to determine whether the altered patterns of laminin expression are related to changes in cell-cell adhesion in tadpoles with heightened tails.
Gremlin 1 is a bone morphogenetic protein (BMP) antagonist, and has an important role in regulating BMPs during lung, limb, and kidney development as well as during neural crest cell differentiation [46,47]. Gremlin 1 is overexpressed in various human tumors [48,49], and induces cell migration, proliferation, and invasion through induction of a fibroblast-like morphology and decreased E-cadherin expression [50]. Gremlin also down-regulates chondrogenesis and programmed cell death [51]. In our experiments, dragonfly treatment induced greatest expression of the gremlin gene (Figure 3), down-regulation of caspase-9, and increased expression of apoptosis inhibitor 5 ( Figure 6). In the salamander treatment, a slight increase in (See figure on previous page.) Figure 6 Hierarchical clustering of 80 known genes showing greater than 5 fold difference compared to control that were induced by dragonfly larvae, and selection of dragonfly specific genes. (a) Expression profiles of the genes using hierarchical clustering by single linkage with Euclidean distance. (b) Expression profiles of the 80 genes and (c) procedure used for selection of predation-threat responsive genes. Genes enclosed by pink boxes are dragonfly-specific, and those enclosed in green boxes are those commonly observed after the salamander treatment. gremlin gene expression was observed, and apoptosis inhibitor 5 was also clearly increased ( Figure 7).
Under hypoxic conditions, HIF-lα isoforms are stabilized and induce HIF-dependent transcription of Vegf, Greml, and Noggin [52]. Gremlin 1 expression is increased in the walls of the small vessels of the pulmonary circulation in vivo during the development of hypoxic pulmonary hypertension [53]. In our microarray experiment, two types of HIF1a gene were used as a probe, and the largest increases in expression (3 fold and 1.3 fold compared to control) were observed in the dragonfly treatment groups; similarly, the highest increase in expression of gremlin 1 (4.5 fold compared to control) occurred in a dragonfly treatment (Table 1). Exposure to a predation threat from larval salamanders also induced an approximately 1.5 fold increase in the expression of gremlin 1. These results may indicate that predation stress involves hypoxic reactions and programmed cell death in vivo.
BCL6 corepressor-like 1 (BCoR-L1) is a transcription corepressor and functional studies have shown that it can bind to class II histone deacetylases, and interact with the CTBP1 (C terminal binding protein 1) corepressor to repress expression of E-cadherin [54]. BCoR-L1 displays homology to several proteins involved in DNA damage repair pathways and in transcription regulation (BCoR) [54]. It is also involved in the regulation of proliferation and apoptosis [55]. Interestingly, we found higher expression of the gene following the removal of the dragonfly larvae or larval salamanders (1.7 and 1.4 fold higher expression than control, respectively; Table 1); only continuous exposure to larval salamanders caused a decrease in gene expression.
The remaining three genes were of unknown identity, although gene #9 showed the highest increase in expression in the salamander treatment group among the nine genes. Therefore, the function of gene #9 might be important in the formation of the bulgy morph in response to the threat of salamander predation.
The IPA analysis showed that 7 of the selected 9 genes were connected directly or indirectly to ubiquitin C ( Figure 4). Therefore, these genes involve the ubiquitinproteasome system. As described above, ELAV-like1 is a multifunctional protein, has been reported to be in the ubiquitin-proteasome pathway in the regulation and function of HuR (an ELAV-like1 protein) [56], and is also related to ubiquitination in predator-induced morphological changes. Further, since ELAV-like 1 protein binds methyltransferase-like 7A and RAD9-HUS1_RAD1 mRNA, then expression of ELAV-like1 may be regulated by these two genes.
A combination of discriminant and IPA analyses identified the involvement of the ubiquitin-proteasome system, suggesting that ubiquitination may have a role to play in predator-induced morphological changes.
Genes showing 5 fold differences in expression compared to control.
Overall, 80 and 81 genes were selected as responsive to the dragonfly larvae and larval salamander predation treatments, respectively (Figures 6 and 7). In these experiments, we were seeking to identify predation-threat responsive genes and 9 were selected as salamanderspecific ( Figure 8b). Two of these genes, zinc finger protein 750 and apoptosis inhibitor 5, were also included in the set of common genes. Three genes, cold inducible RNA-binding protein (5.7 fold reduction compared to control), enhancer of split related 10 (5.9 fold reduction), and tetratricopeptide repeat protein 30 (6.6 fold increase) showed the reverse expression profile in tadpoles exposed to larval salamanders compared to dragonfly larvae. Therefore, these three genes seemed to be further specific for the salamander treatment.
Cold inducible RNA-binding protein (CIRP) is induced by cellular stress such as cold shock [57], UV irradiation [58], oxidative stress or hypoxia [59], and inhibition of neural apoptosis [60]. Down-regulation of CIRP results in decreased cell proliferation [61]. The CIRP homologue in Xenopus, called XCIRP-1 [62], was used as the probe on this Xenopus microarray. Down-regulation of XCIRP-1 inhibits expression of adhesion molecules such as αEand β-catenin, Cand E-cadherin in Xenopus embryos [63]. As we found that exposure to larval salamanders enhances radical production in tadpoles compared to controls (unpublished data), the down-regulation of XCIRP-1 in the salamander treatment group was unexpected. Expression of C-cadherin was decreased in all treatment groups except control, and was 2.5 and 2.0 fold lower than control in the salamander and dragonfly treatment groups, respectively. α(E)-catenin (Xl.6961.1.S1_at) also showed the greatest decrease (2.4 fold compared to control) in the salamander treatment group.
Interestingly, muscle related genes (gene #2, small muscle protein; gene #7, dystrophin; gene #12, myosin light chain) were down-regulated in tadpoles exposed to larval salamanders (Figure 7). In the dragonfly treatment groups, gene #25 (tropomyosin) was down-regulated; however, gene #51 (collagen alpha 2), gene #62 (adhesion regulating molecule), gene #38 (similar to myosin IC) showed increased expression. These results point to a complex response, possibly related to the structural changes in tail morphology, in tadpoles subjected to a predation threat.
Gene #23 (enhancer of split related 10) was downregulated only by the salamander treatment (Figure 8b). The Xenopus homologue of this gene is an isoform of the human gene Hairy and enhancer of split homologue 5 (Hes5) [64], and this negatively regulates cell differentiation during embryogenesis [65]. It is not presently clear whether the predator-induced morphological changes to tadpole tail morphology involve a cell differentiation process.
Tetratricopeptide repeat protein 30 has a role in tubulin glycosylation and glutamylation for maintaining cilia structure and motility [66]. Cilia and basal bodies play important roles in many physiological processes, including cell and fluid movements, sensory perception and development [67], and ion channels that function as chemo-, osmo-, or mechanosensors [68,69]. This gene showed opposite responses in tadpoles of the dragonfly and salamander treatment groups (Figure 8b): a 6.6 fold increase was seen in the salamander treatment group, while a 6.5 fold increase occurred following the removal of the dragonfly larvae compared to control. In the salamander treatment, tadpoles maintain their osmotic pressure by uptake of Na and Cl ions to increase the amount of bodily fluids for bulgy morph formation [30]. Possibly, non-motile or sensory cilia may be involved in this process.
With regard to dragonfly-specific genes, centrin 3, proteasome subunit, and syndecan-2 showed different expression patterns between salamander and dragonfly treatment groups. Interestingly, centrin 3 also showed a somewhat opposite reaction in the salamander and dragonfly treatment groups (Figure 8c). Centrin is conserved in a variety of eukaryotes [70], and is involved with inhibition of cell proliferation [71]. However, it is unclear whether the up-regulation of centrin3 in tadpoles exposed to dragonfly larvae is associated with a change in the rate of cell proliferation.
Syndecans (types 1-4) are a family of heparan sulfatebearing transmembrane proteoglycans involved in cell adhesion, mobility and growth factor interactions [72][73][74]. The syndecan-2 ectodomain promotes focal adhesion and stress fiber formation in fibroblasts in a distinctly different pattern to fibronectin and independently of heparin sulfate requirement [75]. Syndecan-2 also controls assembly of laminin and fibronectin into a fibrillar matrix [76]. However, syndecan-2 down-regulation impairs angiogenesis in human microvascular endothelial cells, and reduces spreading and adhesion of endothelial cells, thereby enhancing their migration but also impairing the formation of capillary-like structures [77]. The data in the microarray analysis showed an approximately 6.2 and 3.2 reduction in syndecan-2 expression compared with the control in the dragonfly and removal of salamander treatment groups, respectively. With regard to bulgy morph formation, it is reasonable to expect the tadpoles to develop a vascular system for maintaining body fluids. However, down-regulation of angiogenesis might be required when the salamander predation threat is removed. Exposure to dragonfly larvae predation may also induce down-regulation of angiogenesis to cause edema; the dragonfly is a top predator and, therefore, tadpoles may need to increase their swimming ability to avoid capture. In addition, as described above, the highest expression of HIF genes was observed in the dragonfly treatment group. This result might also be connected to the impairment of angiogenesis.

Conclusions
In this study, we have identified some key genes involved in the adaptation of tadpole bodies in response to specific predators. These predation-threat responsive genes seem to function in producing morphological changes that depend on the nature of the predation threat; the selected genes might also include causal or associated genes for adaptation to predators. Recently, numerous reports have shown that signal transduction in gene expression pathways are connected directly or indirectly. Therefore, predatorspecific genes might be responsible not only for responses to predation threat but also to many types of stress. Therefore, the common genes listed in Figure 8a might have the potential to be related or indirect genes induced by both predators. However, the single gene depicted as salamander and dragonfly specific genes (Figure 8b and c) might have the potential to be a causal gene for the specific responses to each predator.
Although the threat of predation by larval salamanders or by dragonfly larvae induces the formation of a heightened tail, gene expression profiles in tadpoles exposed to these predators were different suggesting functional differences in the modified tail tissue. It will be of interest to determine what factors stimulate these different patterns of expression of key genes.