Transcriptional response of Mexican axolotls to Ambystoma tigrinum virus (ATV) infection

Background Very little is known about the immunological responses of amphibians to pathogens that are causing global population declines. We used a custom microarray gene chip to characterize gene expression responses of axolotls (Ambystoma mexicanum) to an emerging viral pathogen, Ambystoma tigrinum virus (ATV). Result At 0, 24, 72, and 144 hours post-infection, spleen and lung samples were removed for estimation of host mRNA abundance and viral load. A total of 158 up-regulated and 105 down-regulated genes were identified across all time points using statistical and fold level criteria. The presumptive functions of these genes suggest a robust innate immune and antiviral gene expression response is initiated by A. mexicanum as early as 24 hours after ATV infection. At 24 hours, we observed transcript abundance changes for genes that are associated with phagocytosis and cytokine signaling, complement, and other general immune and defense responses. By 144 hours, we observed gene expression changes indicating host-mediated cell death, inflammation, and cytotoxicity. Conclusion Although A. mexicanum appears to mount a robust innate immune response, we did not observe gene expression changes indicative of lymphocyte proliferation in the spleen, which is associated with clearance of Frog 3 iridovirus in adult Xenopus. We speculate that ATV may be especially lethal to A. mexicanum and related tiger salamanders because they lack proliferative lymphocyte responses that are needed to clear highly virulent iridoviruses. Genes identified from this study provide important new resources to investigate ATV disease pathology and host-pathogen dynamics in natural populations.


Background
Emerging infectious diseases (EIDs) pose a serious threat to the health, stability, and persistence of human and wildlife populations [1][2][3][4]. Genetic and genomic tools have been incredibly useful for discovery of genes associated with host response and variation in resistance or susceptibility to a variety of pathogens [5][6][7]. The advent of genomic tools such as microarray analysis has offered new insights into host-pathogen systems. Additionally, their application to genomic response to host disease response allows rapid characterization of candidate genes for further research into control and eradication methods.
EIDs are a leading hypothesis for the global decline of amphibians and two pathogens in particular, Batrachochytrium dendrobatidis and Ranaviruses have been implicated in worldwide epizootics. Although studies are beginning to investigate possible mechanisms of resistance to these pathogens [8], in general, very little is known about the immune response of amphibians to EIDs. This is because most natural amphibian species are not used as laboratory models and we lack fundamental molecular tools to investigate disease pathology and host-pathogen interactions at the molecular level for all but a few species (e.g., Ambystoma tigrinum spp., Xenopus spp.).
Over the last 15 years, Ranavirus infections have been associated with marked increases in morbidity and mortality in fish, reptiles, and amphibians [9]. Ranaviruses are globally-distributed double-stranded, methylated DNA viruses of fish, amphibians and reptiles and are implicated in amphibian epizootics worldwide [9][10][11]. Both encapsulated and non-encapsulated forms can be infectious. The virus enters the cell via receptor mediated endocytosis or via fusion with the plasma membrane; and DNA and RNA synthesis occur in the nucleus, while protein synthesis occurs at morphologically specific assembly sites in the cytoplasm [9]. In North America, ranaviruses have been isolated from the majority of recent documented amphibian epizootics [12], including from tiger salamander (Ambystoma tigrinum) epizootics in Saskatchewan, Canada [13], Arizona [14], North Dakota, Utah, and Colorado, USA [15,16]. The viral variant that infects tiger salamanders, ATV, is transmitted either via direct contact with an infected animal or immersion in water that contains virus and infected individuals exhibit systemic hemorrhaging, edema, ulceration, and necrosis of the integument and internal organs [13,17,18]. In cases where ATV infection leads to mortality, it usually occurs within 2-3 weeks of exposure, with animals displaying symptoms often between 8-10 days post-exposure. Thus, ATV can rapidly overwhelm the tiger salamander immune response. However, mortality is not always a pathological endpoint because virulence and resistance are known to vary among ATV strains and tiger salamander populations, respectively, as indicated by both laboratory experiments and field observations [19]. Research characterizing the tiger salamander genomic response to ATV is needed to better understand the pathology, virulence, and possible mechanisms of resistance to this emerging disease.
The tiger salamander species complex includes A. mexicanum (Mexican axolotl), a model organism with a growing genomic and informatics resource base [20]. The immune system of the Mexican axolotl has been extensively studied using several classical approaches. Relative to other vertebrate models, the axolotl immune response has been described as immunodeficient [21,22]. There are several reasons for this characterization, including: production of only two immunoglobulin (Ig) classes, only one of which regulates the humoral response and neither of which is anamnestic [23,24]; no response to soluble antigens [25]; poor mixed lymphocyte reactions [26,27]; and lack of cellular cooperation during the humoral immune response as indicated by enhanced humoral immunity following thymectomy or X-ray irradiation [28,29]. Weak immune responses are known for salamanders in general, and the Mexican axolotl and related tiger salamanders are especially susceptible to ATV infections with high observed mortality rates both in the laboratory and in the field. Indeed, an outbreak of ATV in 2003 at the Indiana Axolotl Colony significantly reduced adult stocks before the virus was contained. By way of comparison, adult Xenopus effectively clear close-related FV3 Ranavirus with an immune response that includes an early T-cell proliferative phase in the spleen [30].
To further investigate the axolotl immune response to ATV, we used an Affymetrix custom microarray to identify genes that were significantly, differentially expressed in the spleen. We then compared these genes to a list of genes associated with regeneration that were previously identified from A. mexicanum using the same microarray platform. We reasoned that such a comparison would allow us to filter gene expression responses of humoral cells induced generally in response to injury and stress from those expressed specifically in response to ATV infection. Also, this comparison would potentially identify gene expression signatures associated with cell proliferation in response to ATV, as we have previously identified many cell proliferation probe sets on the Ambystoma genechip that are differentially regulated during spinal cord regeneration [31]. The genes that we describe provide mechanistic insights and new tools to investigate salamander antiviral responses in the laboratory and in natural populations.

Animal care and surgery protocols
Inbred A. mexicanum eggs from a single full-sib mating were obtained from the Ambystoma Genetic Stock Center at the University of Kentucky. Each A. mexicanum egg and larva was reared in an individual container in aquifer water treated with ReptiSafe and changed weekly. Individuals were fed brine shrimp ad libitum for the first four weeks post-hatching and blackworms (Tubifex) ad libitum thereafter. Animals were reared in an environmental chamber on a 12:12 h light:dark cycle at 20°C. At 4.5 months of age, 12 individuals were injected with 100 μl of 10 6 p.f.u./ml of ATV isolated from the axolotl colony and suspended in cell culture medium. This amount of virus was determined to be the minimum lethal dose via injection in previous unpublished experiments (Storfer, unpublished data) and the strain utilized in the experiment was extracted from axolotls that had previously been infected and killed by the virus. Simultaneously, four uninfected (control) individuals were sacrificed in MS222 for spleen and lung removal. Spleens from all animals were flash frozen in liquid nitrogen. The same surgical procedure was performed on four infected individuals following 24, 72 and 144 hours of infection. Spleen tissue was utilized due to its previously noted importance in CD8+ T cell immune responses to Ranaviruses, particularly FV3, in frogs [30] Additionally, spleen is an important immune organ as antigens from the blood are processed in the spleen. Lung tissue was removed for viral quantification as it is an internal organ that can be utilized in early stage virus quantification (Stewart, unpublished data).
During the infection period behavioral observations were taken opportunistically. Total RNA was extracted from spleen with TRIzol (Invitrogen) according to the manufacturer's protocol. RNA isolations were further purified using RNeasy mini columns (Qiagen). The amount of RNA present in each isolate was determined via UV spectrophotometry, and RNA quality was inspected via a 2100 Agilent Bioanalyzer. Sixteen high quality isolates (four replicates at each of four sampling times: 0 (controls), 24, 72, and 144 hours post-infection) were used to make individual-specific pools of biotin labeled cRNA probes. Each of the 16 pools was then independently hybridized to an Amby_001 custom Affymetrix GeneChip (for a more detailed description of the microarray platform see [31] and [32]). The University of Kentucky Microarray Core Facility generated cRNA probes and performed hybridizations according to standard Affymetrix protocols.

Quality Control and Data Processing
All quality control and processing analyses were done in R [33]http://www.r-project.org. We used the Bioconductor package "affy" http://www.bioconductor.org to perform several quality control analyses at the individual probe level [34,35]. These analyses included: (1) viewing images of the log(intensity) values of the probes on each Gene-Chip to check for spatial artifacts, (2) investigating measures of central tendency and dispersion by viewing boxplots and histograms of all the GeneChips, (3) viewing pair-wise M versus A plot matrices for replicate Gene-Chips, and (4) viewing an RNA degradation plot [35] that enables the visualization of the 3' labeling bias associated with all GeneChips simultaneously. Upon conducting these probe level analyses, we background corrected, normalized, and summarized all sixteen GeneChips using the Robust Multi-array Average (RMA) algorithm [36]. Following this, we calculated correlation matrices for replicate GeneChips (four correlation matrices with four GeneChips per matrix; all r from replicate GeneChips > 0.980) on the summarized probe-set level data. The strong correlations observed between replicate GeneChips suggests that we were able to obtain a high degree of repeatability within treatments.

Data Filtering
Microarrays may not accurately quantify the abundance of lowly expressed genes [37]. Calculating statistical tests for such genes adds to the multiple testing burden that is inherent to microarray studies. To address this issue, we filtered genes whose mean intensity across all 16 Gene-Chips was greater than the mean of the lowest quartiles (25 th percentiles) across all GeneChips (n = 16, mean = 5.83, SD = 0.06; data presented on a log 2 scale). Upon imposing this filtering criterion, 3619 probe-sets were available for significance testing.

Identifying Differentially Expressed Genes
We used the Bioconductor package LIMMA [38,39] to generate moderated t-statistics for all six of the possible pair-wise contrasts of the four sampling times investigated in our study. LIMMA employs an empirical Bayes methodology that effectively shrinks the sample variances towards a pooled estimate. This approach reduces the likelihood of obtaining large test statistics due to underestimation of the sample variances. The moderated t-statistics generated by LIMMA test the null hypothesis that the difference between the two groups being compared is zero (i.e., group 1 -group 2 = 0). LIMMA also generates moderated F-statistics that test the null hypothesis that none of the contrasts within a family of contrasts are statistically significant. We corrected for multiple testing by applying the step-up algorithm [40] to the P-values of the moderated F-statistics associated with our six contrasts. Upon correcting for multiple testing, we identified 2322 genes (probe-sets) that were statistically significant. To prioritize amongst differentially expressed genes, we focused on probe-sets that exhibited two-fold or greater changes at any time-point relative to controls. Any gene that was non-significantly down-regulated but significantly up-regulated at one or more time points was considered up-regulated, and vice versa for classification of up-versus down-regulation. We also required that these probe sets have moderated F-statistics greater than or equal to the 50 th percentile of the 2322 F-statistics from the statistically significant probe-sets (F ≥ 12.68). We further limited our analysis to only those probe sets that exhibited significant sequence identity with a human reference sequence. We note that 263 probe-sets with no functional annotation were statistically significant, differentially expressed by ≥ two-fold, and had F-values ≥ 12.68.

Clustering
Hybridization intensities were averaged within treatment groups (0, 24, 72, and 144 hrs post-infection) and log 2 ratios were calculated for each non-zero sampling time relative to 0 hours post-infection. Genesis v. 1.6.0 [41,42] was used cluster these log 2 ratio data and to generate heat maps. Clustering was conducted using a Self Organizing Map (SOM) algorithm. Default conditions were used with the exception that the SOM was allowed to run for 263,000 iterations. The dimensions of the final SOM are 2 x *1 y . These dimensions were determined by comparing output from several different combinations.

Enrichment Analyses
Functional annotation of genes by gene ontology was performed using the Database for Annotation, Visualization and Integrated Discovery (DAVID, [43]). Functional annotation clustering was performed using the default settings with the exception of using the highest classification stringency.

Quantitative real-time PCR
We used quantitative real-time PCR (qPCR) to confirm the results of the microarrays. We estimated a fold change for 24 and 72 hr time points using the ΔΔct method of relative quantification [44], utilizing ribosomal protein L 19 as an endogenous control gene. The same total RNA that was used for microarray analysis was used to create cDNA for qPCR using the BioRad iScript cDNA synthesis kit, following manufacturer instruction. Primers for the qPCR were designed using Primer Express 2.0 (Applied Biosystems). Primers were designed to encompass the sequence of GeneChip probe sets (Additional file 1). qPCR was accomplished using SYBR Green chemistry.
To verify that exposed animals were infected and to quantify viral load and replication over time, we performed qPCR on lung tissue with TaqMan chemistry following the protocol detailed in [45]. ANOVA with a Tukey's HSD correction for all pairwise comparisons was performed to determine if viral loads were significantly different across time points.

Viral load and disease pathogenesis
Viral load for each animal was estimated using qPCR and then averaged for each time point (Fig 1). The significant increase in viral load across time points indicates that animals were infected and that viral replication was occurring. ANOVA with a Tukey's correction for multiple comparisons confirmed that viral load increased linearly between 24, 72, and 144 hours post-infection, and all time points were significantly different from all other time points (F 3,44 = 242.56; p ≤ 0.01).
No animals displayed any gross symptoms of ATV infection in terms of hemorrhaging, lesions or edema, either externally or on any internal organs upon euthanasia and subsequent surgery. Similarly, there were no notable changes in behavior observed during the period of infection. This is likely due to the relatively short infection period utilized in this experiment. As noted in the introduction, infected animals often take 8-10 days, or more, to become symptomatic upon infection.

Gene clustering and functional annotation
We identified 263 probe sets with statistically significant differences in mRNA abundances between Day 0 and any other subsequent time point (Tables 1, 2). We assume that statistically significant probe sets correspond to genes that were differentially regulated after ATV infection. Cluster analysis of the statistically significant genes identified two groups that exhibited similar changes in mRNA abundance. After ATV infection at Day 0, 158 putative genes showed a significant increase in mRNA abundance at subsequent time points (Figure 2), while 105 transcripts showed a significant decrease ( Figure 3). Thus, more genes were up-regulated than down-regulated in response to ATV infection. Overall, DAVID categorized statistically significant genes among 44 different groups that correspond to different biological processes. Eight of these groups contained more genes than would be expected by chance sampling of genes from the microarray (geometric mean p-value < 0.05); these groups were considered significantly enriched with candidate genes relative to other groups ( Table 3). Four of these significant groups contain gene ontologies related to immune response and pathogen response, including innate immunity, complement activation, lysosome function, and antigen processing and presentation. The most enriched functional group contains genes primarily related to immune function and defense responses. The remaining four functional groups contain gene ontologies related to ion binding, ion transport, vitamin metabolism, and response to an unfolded protein. Many genes that were classified in broader biological process categories that are not directly immunityrelated are nonetheless associated with immunity in vertebrates [e.g. [46][47][48]].

Genes Up-regulated in Response to ATV
Across all time points, the majority of up-regulated genes were related to immune response or other related functions, such as inflammation and apoptosis. Other up-regulated genes pertained to gene functions such as ion binding and transport, membrane related functions, and protein binding and modification. Twenty-three genes (represented by 26 probe sets) demonstrated 2-fold or greater changes at 24 hours post-infection, all of which were up-regulated. Ten of these 23 have functions pertaining to immune response. Of the remaining highly expressed genes, one was associated with inflammation, two to regulation of apoptosis, three to ion binding, three to protein binding and modification, one to transport, one to the extracellular constituent, and one to membrane and glycolipids. Many of these genes showed increasing transcript abundances over time. At 72 hours post infection, 43 genes had a greater than 5-fold change, and 40 genes had a greater than 5-fold change at 144 hours. The highest expression level, 91-fold increase at 144 hours, was observed for interferon-induced protein with tetracopeptide repeats 5 (IFIT5).

Genes Down-regulated in Response to ATV
In contrast to the very high fold changes observed among up-regulated genes, the largest fold change observed among down-regulated genes was approximately 4.9-fold, in chondroitin sulfate proteoglycan (NCAN). Five down-regulated genes each code for regulation of transcription and translation. An additional 15 down-regulated genes correspond to 20 probe sets that have functions associated with cell division and mitosis, which was not observed in the up-regulated genes. Other notable down-regulated gene ontologies include one gene corresponding to pinocytosis and endocytosis, and one gene related to natural killer cell mediated cytotoxicity.

Validation of Microarray Results Using Quantitative Realtime PCR
We used qPCR to estimate fold changes for nine genes to verify our microarray data (Table 4). For five of the nine genes investigated (56%; Myxovirus resistance 1, Macrophage receptor with collagenous structure, Complement component 3, Cyclin dependant kinase inhibitor 1B, Vaccinia related kinase 1) there is good agreement between the microarray and qPCR data. In genes where the microarray estimates of fold change were modest (Serine dehydratase like, Hemoglobin gamma alpha, Glycogen synthase kinase, Programmed cell death 8) there is poorer agreement between fold change estimates from these two technologies. However, for this latter group of genes with modest fold change values, the microarray and qPCR data were always within four fold of each other. These results demonstrate that we were able to verify robust differences that were suggested by the microarray data.
Log values of viral particles quantified with quantitative real-time PCR across all time points          Expression profiles for cluster 1

Analyses to identify proliferation gene expression signatures
Comparison of gene expression after ATV and tail amputation identified 25 genes that are significantly up-regulated in both experimental frameworks (Table 5). No significantly down-regulated genes were identified in common. Several of the commonly up-regulated genes appear to be related to humoral immunity, and membrane and extracellular matrix related functions. Additionally, general stress response genes such as heat shock 70 kDa protein 5 were similarly regulated. None of the cell cycle genes that are significantly up-regulated during tail

Discussion and conclusion
Emerging infectious diseases are implicated in the global decline of amphibians and other animals [3,[49][50][51]. There is urgent need to develop understanding of amphibian immunological responses to pathogens and to identify host genes that may be important in disease resistance.
Our study shows that functional genomics provides a means to rapidly meet these needs. We infected Mexican axolotls from the Ambystoma Genetic Stock Center with a viral pathogen that is clearly affecting tiger salamander populations in nature [10,[13][14][15]19]. Our results show that ATV infection induces transcriptional changes of genes that are known to function in vertebrate immunity. Below we discuss the transcriptional response in more detail and suggest hypotheses to explain why ATV is often lethal to axolotls and other tiger salamanders.
We detected significant gene expression changes 24 hours post infection. Many of these gene expression changes likely reflect transcription within lymphocytes, as they are the predominant cell type in the spleen of juvenile and adult axolotls [52]. Indeed, the functions of many of these genes are associated with neutrophil, dendritic, and macrophage cell functions, including cytokine signaling (chemokine (C-X-C motif) receptor 4), phagocytosis and destruction of phagocytised particles (disabled homolog 2, mitogen-responsive phosphoprotein, neutrophil cytosolic factor 2, lysosomal-associated membrane protein 1, RAS homolog gene family, member B), complement (complement factor B, complement component 3), and inflammation (pentraxin related gene, rapidly induced by IL-1 beta, cytochrome B-245 beta polypeptide, n-myc and STAT interactor). Up-regulation of complement components that are known to function in the removal of viral particles, and up-regulation of the stress-associated transcription factor jun-b, clearly shows that ATV induced a humoral gene expression response in the axolotl. Further support for this idea was obtained by comparing ATV-induced gene expression changes to changes identified from a previous microarray experiment using A. mexicanum and the same microarray platform. Twenty-five genes that were up-regulated in response to ATV infection were also identified as significantly up-regulated during regeneration [31]. In both microarray studies, blood was not perfused from tissues prior to tissue collection and it is known that leukocyctes express genes during the early wound-healing phase of spinal cord and limb regeneration. Thus, it seems likely that many of the early gene expression changes that we observed in response to ATV-infection reflect a general, humoral transcriptional response to stress.
In addition to this general humoral response, the gene expression patterns that we observed suggest that the Mexican axolotl manifests an antiviral transcriptional response that is not unlike that observed in other vertebrates. For example, ATV infection clearly induces an interferon-mediated, antiviral response. Although probe sets for interferon genes are not represented on the Gene-Chip, we estimate based upon literature surveys that at least 20% of the significant genes that we identified are known in other systems (in vitro and in vivo) to be involved in interferon-mediated transcription [53][54][55]. These genes exhibited some of the largest fold-changes and include two primary transcription factors that compete to activate (interferon regulatory factor 1, up-regulated) and repress (interferon regulatory factor 2, down-regulated) transcription of interferon-alpha and beta (Type 1 interferon), and inferon-inducible genes that recognize and degrade intra-cellular viral nucleic acid (interferon induced with helicase C domain 1). Considering further that four of the most highly enriched functional groups also contained genes relating to the immune response and pathogen response, the results show that axolotls mount a robust anti-viral response from 24-144 hours post-infection.
Given the robust immunological transcription response that we observed, it is curious why ATV is so virulent to tiger salamanders. In the closely related Ranavirus frog virus 3 (FV3), larval Xenopus laevis succomb to FV3 but adults effectively clear virons and develop lasting resistance to future infection [56]. Adult resistance in X. laevis is correlated with a significant proliferation of cytotoxic CD8 + T cells in the spleen upon infection (within 6 days), as well as increased mortality upon CD8 + T cell depletion [30,57]. Mortality events due to ATV are more significant among larvae in natural tiger salamander populations, however metamorphosed adult tiger salamanders are more susceptible than larvae to ATV infection in the lab [18]. It is well established that Mexican axolotls have a less complicated immune system and never develop the type of mature immune response typical of amniote vertebrates [21][22][23][24][25][26][27][28][29]52]. We did not observe any gene expression changes that would indicate proliferative leukocyte responses in axolotl spleen. Perhaps this is because we used juvenile axolotls that are incapable of such a response. However, it is also possible that ATV maybe more resistant to the immune response mounted by A. mexicanum than FV3 is to the Xenopus immune response.
Phylogenetic analyses indicate ATV is more closely related to iridoviruses found in fish than to FV3, which suggests a relatively recent host switch occurring with the introduction of sportfish to areas of the southwestern United States [15]. Iridoviruses found in sportfish have a larger genome and contain more ORFs related to immune evasion than FV3, which could also be related to improved performance of this virus on the salamander host [15]. Further studies are needed to better understand the ontogeny of immunological responses in axolotls, the virulence of different ranaviruses, and the role of innate versus adaptive immunity in ATV infection.
Our study has identified hundreds of new candidate genes for laboratory and field studies of stress and disease in tiger salamanders. Significantly more gene candidates will undoubtedly be discovered using a higher content, 2 nd generation microarray that is currently under development. Genomic and bioinformatics tools make Ambystoma a powerful system for wildlife disease research.
In particular, molecular information can be quickly crossreferenced from a genetically homogeneous strain that is available for laboratory studies (Mexican axolotl), to other closely related tiger salamander species in North America [20]. Such power is needed to quickly understand how ATV and other pathogens are overwhelming amphibian immune responses and causing population declines in nature.