Global analysis of gene expression in mineralizing fish vertebra-derived cell lines: new insights into anti-mineralogenic effect of vanadate

Background Fish has been deemed suitable to study the complex mechanisms of vertebrate skeletogenesis and gilthead seabream (Sparus aurata), a marine teleost with acellular bone, has been successfully used in recent years to study the function and regulation of bone and cartilage related genes during development and in adult animals. Tools recently developed for gilthead seabream, e.g. mineralogenic cell lines and a 4 × 44K Agilent oligo-array, were used to identify molecular determinants of in vitro mineralization and genes involved in anti-mineralogenic action of vanadate. Results Global analysis of gene expression identified 4,223 and 4,147 genes differentially expressed (fold change - FC > 1.5) during in vitro mineralization of VSa13 (pre-chondrocyte) and VSa16 (pre-osteoblast) cells, respectively. Comparative analysis indicated that nearly 45% of these genes are common to both cell lines and gene ontology (GO) classification is also similar for both cell types. Up-regulated genes (FC > 10) were mainly associated with transport, matrix/membrane, metabolism and signaling, while down-regulated genes were mainly associated with metabolism, calcium binding, transport and signaling. Analysis of gene expression in proliferative and mineralizing cells exposed to vanadate revealed 1,779 and 1,136 differentially expressed genes, respectively. Of these genes, 67 exhibited reverse patterns of expression upon vanadate treatment during proliferation or mineralization. Conclusions Comparative analysis of expression data from fish and data available in the literature for mammalian cell systems (bone-derived cells undergoing differentiation) indicate that the same type of genes, and in some cases the same orthologs, are involved in mechanisms of in vitro mineralization, suggesting their conservation throughout vertebrate evolution and across cell types. Array technology also allowed identification of genes differentially expressed upon exposure of fish cell lines to vanadate and likely involved in its anti-mineralogenic activity. Many were found to be unknown or they were never associated to bone homeostasis previously, thus providing a set of potential candidates whose study will likely bring insights into the complex mechanisms of tissue mineralization and bone formation.

Conclusions: Comparative analysis of expression data from fish and data available in the literature for mammalian cell systems (bone-derived cells undergoing differentiation) indicate that the same type of genes, and in some cases the same orthologs, are involved in mechanisms of in vitro mineralization, suggesting their conservation throughout vertebrate evolution and across cell types. Array technology also allowed identification of genes differentially expressed upon exposure of fish cell lines to vanadate and likely involved in its anti-mineralogenic activity. Many were found to be unknown or they were never associated to bone homeostasis previously, thus providing a set of potential candidates whose study will likely bring insights into the complex mechanisms of tissue mineralization and bone formation.

Background
Vertebrate skeleton is a multifunctional organ providing protection for soft tissues, structural support for muscle and connective tissues, storage for calcium and phosphorus, and a site for haematopoietic cell production and B lymphocyte maturation in adults [1]. Skeletogenesis, i.e. skeleton formation during development, is a process involving complex cellular and molecular mechanisms associated with ossification and bone remodeling. Both processes are known to be responsible for maintaining bone mass and skeletal integrity throughout life [2]. Skeletogenesis requires concerted interplay between various cellular activities (e.g. osteoblast and chondrocyte differentiation [1]) and numerous molecular determinants (e.g. skeletal proteins, growth factors, transcriptional regulators and signaling pathways [2]). Although human and mouse genetics have greatly contributed to unveil the mechanisms involved in skeletogenesis, information remains often insufficient for the successful development of therapies targeting skeletal diseases. Recent studies, reviewed by McGonnell and Fowkes [3], have demonstrated the suitability of fish models to investigate vertebrate development and in particular skeletogenesis. The resemblance of biochemical and physiological processes from fish to mammals, the presence in fish of orthologs for most mammalian genes, the similarities in organ morphology and systems composition are among the traits that contributed to the recent and rapid interest in fish models. Those combined with technical advantages, e.g. large progeny, external reproduction fast growth and translucent larvae, and the power of fish genetics, have definitively transformed fish systems as promising alternatives to mammalian systems. In fact, various fish mutants have already been used to model human skeletal diseases: osteogenesis imperfecta (an autosomal dominant disorder characterized by extreme bone fragility) can be modeled by zebrafish chihuahua mutant [4]; craniofacial syndromes (holoprosencephaly, campomelic dysplasia and Ehlers-Danlos syndrome) can be modeled by zebrafish sonic you, jellyfish, and b4galt7 mutants [5,6], idiopathic scoliosis can be modeled by guppy curve back mutant [7] and important efforts have consequently been made towards the development of fish biochemical, molecular, and cellular tools [8]. These include: (i) the sequence of genomes of various fish models, e.g. zebrafish, green-spotted puffer fish, Japanese medaka and stickleback; (ii) the development of large collections of expressed sequence tags (EST) which were produced for several fish species (e.g. Atlantic salmon [9], rainbow trout [10], Atlantic halibut [11] or channel and blue catfish [12]); (iii) the development of several microarray platforms to explore these EST collections (e.g. Nimblegen Technology high-density oligo-array for the catfish [13], parallel synthesis technology high-density DNA microarray for the Atlantic halibut [14], Agilent Sure-Print™ Technology oligo-array for the largemouth bass [15], the rainbow trout [16], and the gilthead seabream Sparus aurata [17]); and (iv) the development of several fish-derived cell lines (the Fish Cell Line Database, http://www.fcma.ualg.pt/edge/FICELdb.mht). To investigate mechanisms of tissue mineralization, various bonederived cell lines of fish origin are available [18] and gilthead seabream VSa13 and VSa16 cell lines (derived from vertebra) are of particular interest due to their pre-chondrocyte and pre-osteoblast phenotypes. Although they are both capable of mineralizing their extracellular matrix, VSa13 and VSa16 cell lines behave differently regarding their degree of mineral deposition [19], levels of alkaline phosphatase activity [19], expression of mineralogenic genes, e.g. matrix gla protein (MGP), osteocalcin (OC), osteopontin (SPP1), bone morphogenetic protein-2 (BMP-2) [19][20][21], and susceptibility to mineralogenic or anti-mineralogenic molecules, e.g. insulin, IGF-1, vanadate [22,23] and retinoic acid [unpublished results]. They are therefore considered as different bone cell types and common genes differentially expressed during in vitro mineralization should represent key mineralogenic genes, while other differentially expressed genes would represent genes involved in cell type-specific processes not related to mineralization. In this work, we have used an Agilent Sureprint 4 × 44K oligo-array, containing two non-overlapping probes for each of 19,734 unique gene transcripts [17], to analyze global gene expression during in vitro mineralization of these two cell lines. We have also identified in VSa13 cell line the presence of mineralogenic genes whose expression was altered upon cell exposure to vanadate, a molecule with anti-mineralogenic activity [23].

Results
Genes differentially expressed during in vitro mineralization of gilthead seabream vertebra-derived cell lines Confluent cultures of VSa13 and VSa16 cells were cultivated during 4 weeks under control (regular medium) or mineralizing (regular medium supplemented with calcium, phosphate and L-ascorbic acid) conditions. Deposition of mineral nodules within extracellular matrix was confirmed by von Kossa staining in cells exposed to mineralization cocktail ( Figure 1) and total RNA was extracted from three biological replicates per condition. After proper amplification and labeling, each RNA sample was hybridized against the oligo-array and raw expression data were extracted and filtered using Figure 1 Pictures of pre-chondrocyte (VSa13) and preosteoblast (VSa16) cells undergoing mineralization. Cultures were treated for 4 weeks under control or mineralizing conditions then von Kossa stained to reveal mineral deposition. Agilent Feature Extraction 9.5.1 software, then normalized. Quantile normalization showed the highest agreement among replicates (i.e. lowest variation of normalized fluorescence distribution among replicates; data not shown) and was consequently used to normalize raw data sets. Normalized data sets were analyzed through significance analysis of microarray (SAM) and genes differentially expressed in control versus mineralizing conditions were identified. False discovery rate (FDR) threshold was set at 5% and only probes with fold change (FC) over 1.5 were considered. A total of 4,777 and 4,554 probes -corresponding to 3,011 and 3,049 unique genes -indicative of an up-regulated expression were identified from VSa13 and VSa16 RNAs, respectively. Among these genes, 1,489 were shown to be common to both cell lines ( Figure 2). Similarly, 2,359 and 1,642 probes -corresponding to 1,212 and 1,098 unique genes -indicative of a down-regulated expression were identified from VSa13 and VSa16 RNAs, respectively. Among these genes, 469 were shown to be common to both cell lines ( Figure 2). In VSa13 and VSa16 cells, 69% and 49% of differentially expressed genes were simultaneously detected by two non-overlapping probes, respectively, which could indicate the high occurrence of alternative splicing. It could also be the consequence of slight differences of hybridization between the two probes. Indeed, most genes detected by only one probe seem to follow the same pattern of expression when compared to the second non-significant probe. Raw and normalized fluorescence data have been deposited in the GEO database under accession numbers GSE18915 (VSa13) and GSE18941 (VSa16).

Ontology of mineralization-related genes
Genes differentially expressed during in vitro mineralization were classified according to their putative Gene Ontology (GO) -biological process (BP), molecular function (MF) and cellular component (CP) -using OBO-Edit software, AmiGO database [24] and SAPD database [25], based on significant similarity with known genes in public data bases. In general, a GO term could be associated with less than 25% of mineralizationrelated genes (% of BP/MF/CP was 19.2/24.3/18.5 in VSa13 cells and 16.6/20.7/10.3 in VSa16 cells). Comparative analysis of GO classifications indicated that the same type of genes, but not necessary the same genes, were involved in extracellular matrix (ECM) mineralization of both cell lines ( Figure 3). Moreover, similar occurrence of BP, MF and CP classes was observed in genes common to both cell lines (Additional files 1, 2 and 3, Tables S1, S2 and S3). The most represented BP classes ( Figure 3A and Additional file 1, Table S1) were: (i) metabolism (52.9% and 55.9% in VSa13 and VSa16 cells, respectively), (ii) establishment of localization (16% and 16.1%; mostly related to transport), (iii) cellular processes (14.4% and 12.4%, mostly related to signaling) and (iv) regulation (7.8% and 6.9%, mostly associated to cell cycle). Most MF classes were related to binding activity (47.2% and 45.4%; mostly to nucleotides, ions and nucleic acid) and catalytic activity (36.5% and 38.4%) ( Figure 3B and Additional file 2, Table S2). A smaller part was involved in transport, enzyme regulation and molecular transduction activities (total of 10.1% and 10.7% in VSa13 and VSa16 cells, respectively). Finally, most CP classes were related to cytosol or membrane compartments (61.1% and 59.0%) and to specific organelles (14.6% and 17.6%; nucleus, endoplasmic reticulum or Golgi complex). A smaller fraction was associated to macromolecular complexes (9.9% and 8.0%), extracellular region (6.6% and 8.0%) and organelle parts (6.7% and 6.4%). In a general manner, the identification of numerous and diverse genes in both cell lines suggested that ECM mineralization is a complex process that requires tight regulatory mechanisms. In order to understand which type of differentially expressed genes were enriched in these cell lines, we performed a functional annotation analysis using the Database for Annotation, Visualization and Integrated Discovery (DAVID v6.7, available at http://david.abcc.ncifcrf.gov [26,27]) and A two class SAM test was performed; FDR and FC parameters were lower than 5% and higher than 1.5, respectively. Size of diagrams is proportional to the size of gene pools.
focused on processes with fold enrichment > 1.1 and significance p-value < 0.05. In both cell lines, results indicated an enrichment of genes associated with biosynthetic processes, cell cycle and growth, signaling, protein synthesis, stress response, biological regulation and metabolism (Table 1). Other processes such as protein interaction, DNA/RNA metabolism, development, adhesion and bone, were cell type-specific (only identified in VSa16 cells), confirming that although both cell lines are mineralogenic, they also represent different bone cell types.
In an attempt to pinpoint the most relevant mineralogenic genes, we looked at genes with a FC > 10 (FDR was maintained to < 5%). A total of 46/48 up-regulated and 49/49 down-regulated genes were identified in VSa13 and VSa16 cells, respectively. Most of them (61% in VSa13 cells and 59% in VSa16 cells) did not match any known gene (Additional files 4, 5, 6 and 7, Tables S4, S5, S6 and S7). Of these genes, 33 were common to both cell lines, GO analysis indicated that up-regulated genes were mainly associated with (i) transport, (ii) matrix/membrane, (iii) metabolism and (iv) signaling, while down-regulated genes were mainly associated with (i) metabolism, (ii) calcium binding, (iii) transport and (iv) signaling.

Proliferative and anti-mineralogenic effects of vanadate: identification of genes with significant patterns of expression
In order to further investigate genes involved in ECM mineralization, we exposed dividing and mineralizing VSa13 cells to vanadate, an ultra-trace metal with antimineralogenic effect in fish vertebra-derived cell lines, and analyzed global gene expression. Any gene differentially expressed upon vanadate treatments would be potentially important for differentiation/mineralization mechanisms. As expected from previous studies [23,28], cell proliferation was stimulated by vanadate at concentrations up to 7.5 μM, while ECM mineralization was inhibited by vanadate at concentrations up to 5 μM ( Figure 4). Total RNA was collected from three biological replicates of proliferating and mineralizing VSa13 cells exposed to vanadate or left untreated. After proper amplification and labeling, each RNA sample was hybridized against gilthead seabream oligo-array and data were extracted, normalized, and analyzed as previously described. Differentially expressed genes were identified through a two-class SAM analysis with FDR < 5% and FC > 1.5 in the following data sets: i) control versus mineralization: 4,223 differentially expressed genes (3,011 up-and 1,212 down-regulated genes), ii) mineralization versus mineralization + vanadate: 1,136 differentially expressed genes (406 up-and 730 down-regulated genes), and iii) proliferation versus proliferation + vanadate: 1,779 differentially expressed genes (496 up-and 1283 down-regulated genes). Among those genes, 342 were common to conditions i) and ii) ( Figure 5A). In order to identify key mineralogenic genes in VSa13 cells, we looked at genes which expression was oppositely regulated during in vitro mineralization and upon vanadate treatment. These genes and their respective GO categories are listed in Table 2 and could be classified according to a score calculated as the ratio between FC M (control versus mineralization) and FC MV (mineralization versus mineralization + vanadate). Most genes Figure 3 Pie chart representations of GO entries occurrence among genes detected in VSa13 and VSa16 cells. Pie charts represent biological processes, molecular function and cellular component GO entries occurrence among differentially expressed genes in control versus mineralized VSa13 and VSa16 cells. Raw data were normalized using the quantile method and a two class SAM test was performed; FDR and FC were lower than 5% and higher than 1.5, respectively. Biological processes GO entries were associated with 19.2% (VSa13) and 16   were shown to be associated with metabolism, cell matrix/adhesion, signaling and calcium binding. Fewer genes were associated with apoptosis, transport, proteolysis, structural activity, translation and growth. Finally, proliferative and anti-mineralogenic effects of vanadate suggest that its mechanism of action may be associated with cell differentiation. Therefore, genes that followed i) opposite expression between mineralization and mineralization upon vanadate treatment and ii) concordant expression between vanadate treatments during proliferation and vanadate treatments during mineralization should be of particular interest. A total of 136 genes were differentially expressed when comparing data sets of dividing and mineralizing cell cultures exposed to vanadate. Among those genes, 87 genes were also differentially expressed during in vitro mineralization without vanadate ( Figure 5B) and genes that fit into the pattern of expression described above were associated to their respective GO categories and listed in Table 3 according to highest FC values. Most genes were associated with metabolism, signaling and cell matrix adhesion. Fewer genes were related to proteolysis, calcium binding, cell cycle and DNA replication.

Validation of microarray data by quantitative real-time PCR analysis of gene expression
Twelve genes differentially expressed during in vitro mineralization as for microarray data (6 up-and 6 down-regulated genes) were selected for validation by quantitative real-time PCR (qPCR) according to the following criteria: two genes with FC > 10, two genes with 10 > FC > 2, and two genes with 2 > FC > 1.5 for each cell line. Comparative analysis of microarray and qPCR expression data is presented in Figure 6. Analysis of Pearson correlation revealed coefficients higher than 0.9 for both cell lines (strong correlation) when comparing data from microarray probes 1 and 2, and coefficients between 0.4 and 0.5 for VSa13 cells (moderate correlation) and between 0.5 and 0.8 for VSa16 cells (moderately high correlation) when comparing data from microarray and qPCR. The individual analysis of each gene revealed that FC measured by qPCR were always higher (p < 0.05) than that measured by microarray (with the exception of RAR-b in VSa13 cells), suggesting a higher sensitivity of qPCR analysis.
Similarly, 6 genes differentially expressed in mineralizing cells exposed to vanadate (3 up-and 3 down-regulated genes) were selected for validation of microarray Functional annotation tool from the database for annotation, visualization and integrated discovery (DAVID; v6.7) has been used to identify significant biological processes among differentially expressed genes. p < 0.05 (P) and fold enrichment (FE) was higher than 1.1.  data by qPCR (Figure 7). Analysis of Pearson correlation revealed coefficients higher than 0.99 (strong correlation) when comparing data from microarray probes 1 and 2, and higher than 0.9 (strong correlation) when comparing data from microarray and qPCR. In general, correlation coefficients were shown to be higher among genes analyzed in vanadate-treated VSa13 cells than among those analyzed in mineralizing VSa13 and VSa16 cells. This difference could in part be explained by the lower FC values observed among genes identified in vanadate-treated cells and consequent lower tendency for FC compression of oligoarray data.

Gilthead seabream oligo-array is a suitable tool to analyze global gene expression of mineralogenic gilthead seabream cell lines
In the present study, a comprehensive set of data has been produced following hybridization of an Agilent 4x44K oligo-array with RNA samples prepared from control, mineralizing and vanadate-treated gilthead seabream vertebra-derived cell samples. A two-class SAM analysis identified mineralogenic genes, some of which had already been previously investigated during in vitro mineralization of VSa13 and VSa16 cells through a candidate gene approach, e.g. tissue non-specific alkaline phosphatase (TNAP; unpublished data), BMP-2 [21], SPP1 [20], MGP and OC [19]. At first glance, a good correlation (i.e. similar type and extent of regulation) between microarray and pre-existing data was observed. This good correlation was later confirmed through qPCR analysis, although a tendency for FC value compression of oligo-array data was also noted. A similar fold-change compression was observed by Wang et al. while validating two commercial long-oligonucleotide microarray platforms, Applied Biosystems and Agilent, through large scale qPCR analysis of gene expression [29], and more recently by Ferraresso and colleagues while using Agilent gilthead seabream oligo-array [17]. Although our microarray data exhibited a dynamic range inferior to that of qPCR data, a situation which was somehow expected, they were generally accurate; Agilent seabream oligo-array therefore represent a valuable tool for simultaneous analysis of the expression of thousands of genes.

In vitro mineralization recruits similar genes in fish and in mammalian bone-derived systems
Putative mineralogenic genes, i.e. those differentially expressed during in vitro mineralization of VSa13 and VSa16 cells, presented a similar distribution into GO categories. Approximately half of those genes were shown to be common to VSa13 and VSa16 cells, indicating that both cell lines, although representing different cell types, recruit similar genes and processes during mineralization. Some biological processes were however differentially enriched in both cell lines. In particular, genes related to GO categories bone mineralization, biomineral formation, remodeling, ossification and skeletal development were enriched in mineralizing     bone-derived cell lines. Global analysis of gene expression of ATDC5 cells -mouse pre-chondrocytes similar to VSa13 cells -and MC3T3-E1 cells -mouse pre-osteoblasts similar to VSa16 cells -identified mineralogenic genes associated with catalysis, signal transduction, transport, transcription, structure and motor activity [30] and metabolism, cell cycle, signaling, extracellular matrix, immune response and transcription [31][32][33], respectively. Similarity in patterns of gene expression of mammalian and fish pre-chondrocyte and pre-osteoblast cell lines, suggested that mechanisms of tissue mineralization might be conserved among vertebrates but also among mineralogenic cell types.

Anti-mineralogenic activity of vanadate as a way to identify key/novel genes involved in mineralization
Although most of the genes identified in the first step of this analysis certainly play a role during in vitro mineralization, some of them must be more important than others. These key genes were identified from the initial bulk of genes by using the anti-mineralogenic activity of vanadate [22,23,34]. Indeed, those genes which expression levels were oppositely regulated during in vitro mineralization and upon treatment with vanadate were considered as good candidates. Vanadate stimulates proliferation of VSa13 cells and strongly inhibits its ECM mineralization, and these processes seem to involve MAPK and putative PI-3K\Ras\ERK pathways [22,23]. Genes differentially expressed under these conditions could represent new candidate genes crucial for bone formation. Moreover, vanadium compounds have long been known for their insulin-like properties [35] and role in bone formation [34,36] but to the best of our knowledge, their effects on gene expression have never been investigated, in particular in relation to bone. DAVID functional annotation tool for KEGG pathways identified genes in vanadate-treated VSa13 cells associated with insulin signaling pathway: 3-phosphoinositide-dependent protein kinase-1 in proliferating cells and Ras homolog gene family (member Q) in differentiating cells. The involvement of signaling pathways related to insulin activity is consistent with insulinmimetic properties of vanadium compounds [35] and recent studies showing that mechanisms of action of vanadate and insulin are similar in fish VSa13 cells and that both molecules exhibit an anti-mineralogenic activity [23]. Among the genes oppositely regulated during in vitro mineralization and upon treatment with vanadate (see Table 3), two have been associated to extracellular matrix and previously shown to play an important role in ECM structure: tenascin (normally expressed in mesenchymal stem cells and osteoblasts) and thrombospondin (normally expressed in mesenchymal stem cells and chondrocytes). Both have been associated to fracture healing, spinal curvature and craniofacial defects in knock-out mice [37]. Rap1b, intermediate of MAPK (among other pathways), ADP-ribosylation factor 5, GTP-binding and effector of phospholipase D signaling, and cyclin-dependent kinases regulatory subunit 1, a Ras effector protein, were also among the genes listed in  in bone formation [1], and were recently shown to modulate/antagonize bone morphogenetic protein activity in transgenic mouse and zebrafish [40,41]. Our data showed an opposite regulation of SCUBE-like and BMP-2 (associated with bone formation [21,42]) gene expression in mineralizing and vanadate-treated cells, suggesting that SCUBE-like protein may play a key role in antimineralogenic activity of vanadate through its action on BMP-2 gene and/or protein. Further studies should however be carried out in order to confirm this hypothesis. Notably, numerous genes detected in this study were classified as unknown. Absence of orthologs in other vertebrate species, high divergence of fish genes and/or low level of annotation in fish sequence databases are likely to contribute to explain this situation. In addition, the fact that numerous genes identified throughout this work have not been previously linked to bone formation, suggests that genetic mechanisms involved in ECM mineralization and bone formation, whether in mammalian or fish species, are still poorly understood.

Conclusions
Global gene expression has been analyzed during ECM mineralization of gilthead seabream vertebra-derived cell lines using a recently developed oligo-array. A considerably high number of differentially expressed genes was detected, and occurrence of GO categories was found to be similar in both cell lines, with approximately half of the genes common to both cell lines. When comparing occurrence of GO categories in VSa13 and VSa16 with those found in mammalian systems, the similarities found suggested conservation in mineralization-associated processes across vertebrates. Interestingly, enrichment for bone-related genes was observed in VSa16, but not in VSa13, thus reinforcing the previously described association of this cell line to osteoblast lineage. Furthermore, analysis of genes differentially expressed upon exposure to vanadate, a known anti-mineralogenic molecule, has permitted the identification of key/novel mineralogenic genes, which could be classified as: i) annotated genes with known roles in bone formation (e. g. tenascin and thrombospondin), ii) annotated genes with unknown roles in bone formation (e.g. SCUBE-2) and iii) unknown transcripts. Although further analyses are required for genes included in the last two categories, the large number of transcripts detected in this study should bring new insights into the process of mineralization.

Cell culture and ECM mineralization
VSa13 and VSa16 cells were cultured in Dulbecco's modified Eagle medium (DMEM) supplemented with 10% fetal bovine serum, 2 mM L-glutamine, antibiotics and antimycotics (all from Invitrogen), as described previously [19]. ECM mineralization was induced in confluent cultures by supplementing medium with 50 μg/ml of L-ascorbic acid, 10 mM of β-glycerophosphate and 4 mM of calcium chloride (all from Sigma-Aldrich). Culture medium was renewed every 3.5 days. At appropriate times, mineral deposition was evaluated through von Kossa's staining and densitometry analysis [43].

Preparation of vanadate solution
Vanadate stock solution (5 mM, pH 6.7) was prepared from ammonium metavanadate (Sigma-Aldrich) and stored at 4°C.

RNA extraction and purification
Total RNA was extracted from cell cultures as described by Chomczynski and Sacchi [44]. RNA samples were purified using QIAGEN RNeasy Mini kit then treated with QIAGEN RNase-free DNase according to manufacturer's instructions. RNA concentration was determined by spectrophotometry (NanoDrop ND-1000, Thermo Scientific) and RNA integrity evaluated by electrophoresis (2100 Bioanalyzer, Agilent Technologies). RNA integrity number (RIN) index was calculated for each sample using Agilent 2100 Expert software. Only RNA samples with a RIN >8 were further processed.

RNA amplification, labeling and array hybridization
Total RNA (500 ng) was supplemented with a mixture of 10 different viral polyadenylated RNAs (Agilent Spike-In mix) then linearly amplified and labeled with Cy3-dCTP using Agilent One-Color Microarray-Based Gene Expression Analysis protocol. Labeled cRNA was purified using QIAGEN RNeasy Mini kit, and sample concentration and specific activity (pmol Cy3/μg cRNA) were measured by spectrophotometry. Labeled cRNA (1650 ng) was fragmented by adding 11 μl of 10X blocking agent and 2.2 μl of 25X fragmentation buffer, heating at 60°C for 30 min, and finally adding 55 μl of 2X GE hybridization buffer. Hybridization solution (100 μl) was placed in the gasket slide and assembled to the microarray slide (each slide containing four arrays). Slides were incubated for 17 h at 65°C in an Agilent hybridization oven, then dissociated from the hybridization chamber and quickly submerged in GE wash buffer #1 for 1 min. An additional wash was performed in prewarmed (37°C) GE wash buffer #2 for another 1 min. Slides were scanned using Agilent G2565BA DNA microarray scanner. Scan resolution was set to 5 μm and two different sensitivity levels (XDR Hi 100% and XDR Lo 10%) were used. Both images were analyzed simultaneously using the standard procedures described in Agilent Feature Extraction software 9.5.1.