Development of a Pacific oyster (Crassostrea gigas) 31,918-feature microarray: identification of reference genes and tissue-enriched expression patterns

Background Research using the Pacific oyster Crassostrea gigas as a model organism has experienced rapid growth in recent years due to the development of high-throughput molecular technologies. As many as 56,268 EST sequences have been sequenced to date, representing a genome-wide resource that can be used for transcriptomic investigations. Results In this paper, we developed a Pacific oyster microarray containing oligonucleotides representing 31,918 transcribed sequences selected from the publicly accessible GigasDatabase. This newly designed microarray was used to study the transcriptome of male and female gonads, mantle, gills, posterior adductor muscle, visceral ganglia, hemocytes, labial palps and digestive gland. Statistical analyses identified genes differentially expressed among tissues and clusters of tissue-enriched genes. These genes reflect major tissue-specific functions at the molecular level, such as tissue formation in the mantle, filtering in the gills and labial palps, and reproduction in the gonads. Hierarchical clustering predicted the involvement of unannotated genes in specific functional pathways such as the insulin/NPY pathway, an important pathway under study in our model species. Microarray data also accurately identified reference genes whose mRNA level appeared stable across all the analyzed tissues. Adp-ribosylation factor 1 (arf1) appeared to be the most robust reference for normalizing gene expression data across different tissues and is therefore proposed as a relevant reference gene for further gene expression analysis in the Pacific oyster. Conclusions This study provides a new transcriptomic tool for studies of oyster biology, which will help in the annotation of its genome and which identifies candidate reference genes for gene expression analysis.


Background
The Pacific oyster Crassostrea gigas, is one of the world's most economically important bivalves (4.2 million metrics tons annual production, worth $ 3.5 billion) [1]. It is also one of the best characterized models for biochemical, molecular and genetics studies among the Lophotrochozoa. Crassostrea gigas, and other oyster species, play an important role in estuarine and marine coastal habitats, where increasing human activity is causing environmental degradation [2]. In these ecosystems, Pacific oysters suffer summer mortalities due to increased stresses and disease outbreaks [3]. At the other extreme, C. gigas can become invasive in new habitats (e.g. in northern Europe) [4]. The scientific rationale of obtaining a wide variety of sequences and developing genomic tools for C. gigas is to take advantage of its membership of the Lophotrochozoa, an understudied clade of bilateria. Its study will contribute to knowledge in the fields of functional, comparative and evolutionary genomics by throwing new light on genome function and diversity [5][6][7][8][9], the evolution of sexuality [10][11][12][13] or immunity [14][15][16]. In the context of environmental genomics, oysters might also prove to be an attractive model organism for understanding population responses to environmental stresses and adaptation through genetic change, as well as for deciphering the physiological bases of complex traits (growth, reproduction, survival) [17][18][19][20].
In July 2010, 56,268 EST sequences were compiled in the publicly accessible GigasDatabase version 6 housed at Sigenae (http://public-contigbrowser.sigenae.org:9090/ Crassostrea_gigas/index.html) [21]. This genomic resource allows the development of genome wide microarrays for screening physiological traits and responses to the environment. DNA microarray technology is a high-throughput method for measuring the expression levels of thousands of genes simultaneously. This approach has been applied extensively to establish gene expression patterns in many organisms, including yeast, worm, human, fruit fly and rice [22][23][24][25][26]. Investigators have already analyzed gene expression in Crassostrea virginica and Crassostrea gigas using cDNA microarrays constructed from 6,780 [27] and 9,058 [21] genes. These microarrays were used to describe gene expression patterns in response to heat stress [28] and hypoxia [29], and to compare oysters that were resistant or susceptible to summer mortalities [21]. Since these studies were done, a large database of ESTs has become available, providing the information necessary to design a more comprehensive microarray suitable for characterizing genome-wide transcriptome profiles in Crassostrea gigas.
Here, for the first time, we describe the design of an oligo-microarray covering a total of 31,918 contig sequences, and its application to monitoring gene expression profiles in multiple tissues of adult oysters. In addition to the first obvious advantage offered by this new microarray, which is the significant increase in the number of contigs per array, this Agilent oligonucleotide microarray uses shorter probes (60 mer instead of an average mean of 413 bp in C. gigas cDNA arrays), resulting in higher density arrays and cheaper manufacture. In addition, with respect to cDNA arrays, oligonucleotide microarrays are reported to provide higher sensitivity and reproducibility [30]. This work aims to survey variation in gene expression across multiple oyster tissue types, to provide experimental evidence for gene function assignments, and to identify gene clusters related to tissue-specific processes. Knowledge of tissue-specific expression is a prerequisite to any studies of organism development, normal functioning and response to injury and disease. Furthermore, the comparison of significantly different transcriptomes allowed the identification of genes with stable expression across tissues. Quantification of gene expression has become a crucial step to investigating any molecular mechanism or physiological process in oysters. However, the quality of normalized quantitative RT-PCR data is dependent on the accuracy of the normalizer itself. To avoid misinterpretation of the differential expression profile of a target gene, it is therefore necessary to identify reliable reference genes and test their stable expression across the different experimental conditions tested. Because microarray screening simultaneously unravels the expression of a huge number of genes, it offers the opportunity to discover genes with transcriptional stability that can serve as pertinent internal reference standards for future quantitative RT-PCR studies on oysters.
The present genome-wide analysis of gene expression in C. gigas tissues does not only establish lists of genes expressed in the surveyed tissues, but also provides evidence for the level of expression of these genes in different tissues. Importantly, this approach may provide clues for elucidating the functions of genes underlying specific processes and identify candidate genes predicted to regulate traits of interest. This argument is particularly relevant for the unannotated ESTs, for which such a survey can provide clues to resolving the orphan status.

Robustness of C. gigas microarray data
We designed a microarray using 31,918 oligonucleotide 60-mer probes representing the presently known genes in the Pacific oyster Crassostrea gigas (GigasDatabase, version 6). Of the 31,918 probes printed onto the array, 20,170 (63.2%) corresponded to annotated contigs in the GigasDatabase. This microarray is intended for use in studies of oyster transcriptomes at the whole-genome level. In our study, the microarray was used to analyze the transcriptome of 9 tissues using one-channel hybridization (Cy3). The custom microarray design and the raw and normalized data is available from the Gene Expression Omnibus (GSE26265; http://www.ncbi.nlm.nih.gov/geo/ query/acc.cgi?acc=GSE26265).
The homogeneous hybridization of the arrays was first estimated by comparing the expression profiles of duplicated probes. Thirty-four probe pairs were printed onto the arrays (Additional file 1). Expression values were highly reproducible with a mean R 2 of 95.7% (Additional file 1). An R 2 of 98.4% was obtained by removing data points with expression values below background levels.
Next, to estimate the robustness of the microarray data, we compared the expression data obtained from different probes with identical annotation. Indeed, the oyster transcriptome is currently not completely sequenced and assembled and, therefore, some non full-length mRNAs are represented by numerous oligonucleotide probes printed on the array. This redundancy was used to confirm inter-and intra-array reproducibility. The consistency of the expression profiles obtained from different probes expected to come from the same mRNA confirmed the reproducibility of the data. For instance, among the 16 contigs originally annotated as similar to the myosin heavy chain of striated muscle from Aequipecten irradians (also known as Argopecten irradians) [GenBank: P24733] (GigasDatabase, version 6), 14 showed extremely similar expression profiles (Mean R 2 = 97.94%) ( Figure 1A). The two contigs with different expression profiles showed very low gene expression values, close to or below background level, suggesting a low probe efficiency or a low expression in the different tissues. The most recent assembly of the C. gigas sequences (version 8, March 2011) constituted three contigs and one singleton from these 16 contigs ( Figure 1B). Two contigs (BQ426757 and CU686207) encode two different myosin heavy chain isoforms (83% nucleotide identity). Contig CU686461 and singleton AM853364 encode two different isoforms of catchin transcript (86% nucleotide identity). Catchin is a protein present in molluscan catch muscles, which is produced by alternative splicing of myosin heavy chain [31,32]. In our dataset, myosin heavy chain and catchin present very similar expression profiles across tissues and a higher expression in catch muscle ( Figure 1A).
Such strong profile consistency was frequently observed in our dataset such as the 3 sequences that matched the uncharacterized protein ZK643.6 (Caenorhabditis elegans) [GenBank: P30652]. Thus, our microarray data may prove useful in providing additional clues suggesting which contigs could constitute a single gene. The genes represented by multiple contigs can also be used to assess the reproducibility of the data and the homogeneity of the hybridization in future analysis. In Additional file 2 we provide a list of 200 putative genes significantly differentially expressed among tissues (ANOVA p < 0.05) and represented by 2 to 11 probes on the array. Over all contig pairs, we obtained a strong correlation of expression profiles (Mean R 2 = 83.2%).
Quantitative RT-PCR was performed for four selected candidate reference genes and six tissue-enriched genes (primers listed in Additional file 3). The patterns of transcript abundance detected for these genes in the array and in quantitative RT-PCR showed extremely similar profiles (mean R 2 = 85%) (Additional file 4). The correlation increased to 93% when we removed the hemocyte-enriched gene (BQ426482), which presented the lowest measured correlation between microarray and quantitative RT-PCR data (63%). Quantitative RT-PCR is commonly used as a validation tool for confirming gene expression results obtained from microarray analysis; however, microarray and quantitative RT-PCR data often disagree [33]. Using multiple probes corresponding to the same contig and comparing their resulting expression profiles may thus be a more accurate way of individually assessing microarray results and oligonucleotide specificity. A strong correlation between multiple probes corresponding to the same gene, as obtained in the present work, reflects high microarray quality.

Prevalent gene expression patterns
We applied principal component analysis (PCA) to the whole dataset. The PCA of all genes in the surveyed tissues showed that datasets acquired from the same tissue type grouped together, indicating that there is a strong homogeneity in the gene expression pattern in each tissue type ( Figure 2). Furthermore, tissues with similar functions had a lower scattering, as observed for the transcriptomes of gills and labial palps. Indeed, gills and labial palps are two pallial organs, both secreting mucus and possessing similar structures involved in particle selection, which is an essential part of bivalve filter-feeding.
It was not possible to extract a sufficient amount of total RNA for microarray hybridization from a single visceral ganglion or from hemolymph of a single individual. Therefore, we randomly pooled samples from 6 individuals. Pooling results in biological averaging for most genes. Consequently, the variation among samples in these two tissue categories is much lower than the variation among individuals for other tissues ( Figure 2). However, Kendziorski et al. [34] demonstated that the construction of a large number of pools enables the biological variation to be properly assessed and comparative studies to be made.
To further explore relationships among tissue samples, we performed a hierarchical clustering analysis on all the genes printed on the array. The results were extremely similar to those of the PCA analysis ( Figure 3A). Gonads grouped together on a single branch, separated from all other tissues with 100% bootstrap support, demonstrating that the transcriptome of gonads (regardless of gender) is significantly different from the transcriptome of somatic tissues. Hierarchical clustering then further separated individual gonad samples according to their gender. The data also showed that the digestive gland and adductor muscle clearly had different transcriptomes, with 100% bootstrap support. Mantle tissue clustered apart from visceral ganglion and hemocyte samples with 60% bootstrap support. However, the hierarchical cluster analysis could not differentiate between labial palps and gills. Interestingly, labial palps and gills from the same individual always clustered together (1Pm with 1GIm, 5Pf with 5GIf and 6Pf with 6GIf) on the hierarchical clustering tree, once again demonstrating the strong similarities between these two tissues ( Figure 3). Further analysis revealed that these organs co-expressed genes in relation with their common biological functions (see below).

Tissue-enriched gene expression
In the absence of a functional assay, tissue-specific gene expression features have often been viewed as indicators of tissue-specific function. In order to identify genes with differential tissue-enriched expression patterns, we initially performed a one-way ANOVA with a p-value  threshold < 0.01, adjusted with Bonferroni's correction. A total of 7,586 contigs were differentially expressed between tissues (Table 1). These differentially expressed contigs were further filtered to determine whether they were overexpressed in a particular tissue. We required that the mean of the log 10 normalized expression level of the replicates in a given tissue was at least 1.1 fold higher than the means in each of the other tissues (see Methods). A 1.1 fold increase of log 10 values is equivalent to a 1.26 fold increase of the raw data. By this process, tissue-enriched contigs were identified for all tissues analyzed ( Figure 3B). The number of tissue-enriched contigs varied substantially, ranging from 23 contigs in labial palps to 471 in the digestive gland. The number of tissue-enriched contigs and the proportion of contigs annotated in the GigasDatabase are shown in Table 1. Figure 3B shows the heat map of the tissue-enriched contigs, and the complete list with their descriptions are given in Additional file 5. The highest percentage of contigs with annotation was obtained for genes specifically expressed in the digestive gland (51%). This is almost certainly because digestive processes are ancient and highly conserved. In digestive glands, many enzymes are widely distributed and conserved over a large range of species. For example, alpha-amylase, a key carbohydrase enzyme for starch digestion has been characterized from bacteria, fungi, plants and animals [35] and high sequence identity was observed for the domains responsible for its enzymatic activity [36]. It is also possible that there were qualitative differences in the original sequencing of mRNAs from the different tissues, which would result in poor BLAST matches. The constant updating of GigasDatabase will allow improvement of the sequence descriptions.
The annotated tissue-enriched contigs showed a strong relationship with the physiological functions of the corresponding tissues (Additional file 5). For example, numerous genes involved in digestion were highly expressed in the digestive gland (clh2, chia, plb1, pnlip, muc2, prss7,    , some of the non-annotated ganglia-specific contigs were found to encode the precursors of putative novel neuropeptides [CU995163, CU992964, AM868259] that the BLAST tool did not identify, probably due to limited sequence conservation between peptides or their precursors [37]. These putative novel neuropeptides displayed the conventional dibasic (KR, RR) cleavage sites for prohormone convertases and potentially downstream GKR sequencing serving as combined amidation and proteolytic signals. The annotation of a high proportion of the genes printed onto this array facilitates the exploration of specific biological processes and the identification of genes potentially involved in the same pathways. For example, of the genes specifically expressed in the visceral ganglia, we found that neuropeptide Y [CU983945] and insulin [CU992662, AM855819 and CU987888] clustered closely ( Figure 3D). NPY and insulin are involved in the processes of the food intake regulation, allocation of energy to growth, reproduction and basal metabolism in both vertebrates and invertebrates [38][39][40]. Such regulatory molecules deserve particular attention, especially in the context of the summer mortality syndrome in which reproductive effort plays a crucial role through its effects on the capacity of oysters to survive [13,15,41]. Overall, genes without significant annotation were specifically attributed to a tissue based on their expression pattern, which can provide valuable help towards the assignment of gene function. Indeed, even if gene function is primarily assessed on the basis of altered phenotypes associated with gene disruption, knowledge of spatio-temporal expression patterns of a given gene represents valuable information for assigning a putative function.
Some functional clusters also showed similarities between tissues, such as male and female gonads, both of which expressed significantly higher levels of 341 genes than all other tissues, or labial palps and gills, which showed significantly higher expression of a common set of 101 genes. Hemocytes and visceral ganglia shared 54 overexpressed genes ( Figure 3C). The complete list of dual tissue-enriched contigs and their descriptions are provided in additional file 6. As expected, male and female gonads showed the highest expression levels of all cell cycle, mitosis and meiosis regulation genes, including cell cycle checkpoint (rad1, fbxo43, mad2a, clspn and cdkn3), control of replication (mcm3, -2 and -7, rfc3, -4 and -5, Gins1, -2 and -3) and DNA repair genes (mre11, xrcc1 and -3, rad54b, -51b, -50 and -21, msh6, lig1 and top2b). Gills and labial palps, the main feeding organs of oysters, expressed 37 genes potentially involved in epithelia morphogenesis (mark2, plscr2 and pfn4), cilia movement (vil1, cam and srI) and in detoxification and defense mechanisms (gst, ces3, gpx1, adamts18 and pox1). The visceral ganglia and hemocytes mainly shared a high expression in topoisomerase I (6 genes), two G protein-coupled receptors (nmur1 and agtr1) whose vertebrate orthologues bind neuropeptides mediating blood pressure regulation and stress response [42,43], and a set of defense-related genes (big defensin, pxn-1, hsp27 and piap) whose expression in the central nervous system suggests a concerted activity with hemocytes in defense and stress response [14,44].

Reference genes
The term "housekeeping gene" is generally applied to genes ubiquitously expressed at the same level in all tissues. These constitute the basal transcriptome for the maintenance of basic cellular functions. Their uniform expression means they can be used as reference genes for quantitative RT-PCR [45]. Thus, we investigated the stability of gene expression levels to identify appropriate reference genes. To allow for comparisons of means and standard deviations, logarithmic normalized data were centered on their mean (fixed at 4). A candidate reference gene was defined as a gene that was not differentially expressed when tested by ANOVA (p-value < 0.05) and which had a coefficient of variation (CVarray, the ratio of the standard deviation to the mean) of < 5%, a maximum fold change (MFC, the ratio of the maximum and minimum values observed within the dataset) of < 1.1, and a mean expression level higher than the average Remarkably, none of the 3 most commonly used reference genes in oyster (actin, g3apdh and ef1a) ranked within the top 125 identified candidate reference genes.
Actin and g3apdh had p-values < 0.05 when tested by ANOVA and were not ranked (Table 3). Ef1a did not show any significant variation when tested by ANOVA but, similarly to actin and g3apdh, it had an MFC > 1.1. Even if the MFC had not been taken into account, ef1a would still have ranked only 256 st among the candidate genes (Table 3).
To demonstrate the importance of identifying appropriate reference genes, we designed primers for four of the 15 top-ranked novel reference genes (hkg4, hkg8, hkg12 and hkg14) and compared their mRNA levels estimated by quantitative RT-PCR with those of the three commonly used reference genes (actin, ef1a and g3apdh) in all tissue samples. To normalize data across samples, the mean of Ct values across all genes was calculated for each sample. The log values (mean Ct /Ct sample) were then compared for each gene. Box plots showing the minimum, maximum and mean value of Ct and log value (mean Ct/Ct sample) for each gene are shown in Figure 4A and 4B. In both datasets, hkg4, annotated as an adp-ribosylation factor 1 (arf1) in GigasDatabase, and g3apdh were the most stable genes, with coefficients of variation (CV; ratio of the standard deviation on the mean) of 0.04 for hkg4 and 0.05 for g3apdh. They also had the lowest CVs among all genes tested ( Figure 4C). Applied to the raw data, Normfinder identified g3apdh as the best reference gene over the ensemble of the tissues, whereas it identified hkg4 as the best reference gene when considering normalized data. Using geNorm analysis on non-normalized expression levels, we determined the gene-stability expression measures (M) for the seven candidate reference genes. Again, hkg4 and g3apdh were the most stable reference genes, with an average expression stability value of 0.041 ( Figure  4D). Arf1 has also been identified as a housekeeping gene in humans [46]. ADP-ribosylation factors are highly conserved 21 kDa GTPases involved in vesicular trafficking in all eukaryotes [47]. Glyceraldehyde-3-phosphate dehydrogenase (GAPDH) is one of the most commonly used housekeeping genes for comparisons of gene expression data, though variability of GAPDH expression between human tissue types was recently shown [48]. In C. gigas, g3apdh was found to be significantly more expressed in the adductor muscle (ANOVA p < 0.05) and should, thus, be avoided when studying gene expression across different tissues. However, its expression appeared stable within tissues suggesting that it may be suitable for other experimental designs. We therefore identified ADP ribosylation factor 1 (arf1) as being a more relevant reference gene for use in normalization of real time PCR for tissueexpression studies than those typically used. In addition, the list of the 125 potential reference genes will clearly help in the choice of references genes in further gene expression studies, but we still emphasize the need to experimentally validate reference genes for specific experimental designs [45].

Conclusions
Research on gene function and oyster physiology using real time PCR and microarray is becoming increasingly important and widespread. Because of their economic importance, evolutionary position and use as a model organism in environmental genomics, more research is now being conducted on oyster species. A large EST project lead by French and US teams generated a genomewide EST database, GigasDatabase, that is now freely available and housed on the Sigenae website. In the present study, the GigasDatabase was used to generate oligomicroarrays, providing the genomic tools necessary to undertake high-throughput studies of significantly important phenomena, such as summer mortalities in bivalves, host-pathogen interactions, reproduction and the effect of polyploidy on physiology of Pacific oysters. However, in order to make sense of gene expression data, we developed a 32 K microarray containing all known available contigs of C. gigas, thus aiming to represent the whole transcriptome. This study fulfilled 3 main goals. (1) It designed and presented a genome-wide oyster microarray that can now be employed in different gene expression studies in Pacific oysters (GSE 26265). The microarray is publicly available, following a request to the Agilent company and an authorization given by our team to share and use the oyster design. (2) It described the basic transcriptome of nine different tissues, providing important information for gene annotation and function prediction for genome sequencing. (3) It identified a large array of genes that can be used as reference genes to normalize gene expression data.

Animal sampling
Twenty-four adult oysters (3 years old) were sampled in July 2010 in Baie des Veys (Normandy, France). First, hemolymph was extracted from the adductor muscle using a 23-gauge needle attached to a 1-ml syringe, then gonad, gills, digestive gland, the posterior adductor muscle, labial palps, mantle and visceral ganglions were dissected. At this stage, pools of 6 individuals were randomly constituted for hemolymph and visceral ganglion samples. All samples were snap-frozen in liquid nitrogen, and stored at -80°C. Total RNA was extracted from individual tissues (3 female gonads, 3 male gonads, 4 gills, 4 digestive glands, Table 2 The 49 annotated candidate reference genes ranked in order of their coefficient of variation (spot: location on microarray, ID genbank: accession number, stdev: standard deviation, description: best BLAST hit description as described in GigasDatabase, e value: best BLAST hit e-value, CV: coefficient of variation, max: maximum value, min: minimum value, MFC: maximum fold change, p value: one-way ANOVA P > 0.05) (Macherey Nagel) isolation columns. After the homogenisation of the N 2 grinding powder in Tri-Reagent and the first centrifugation, the aqueous phase was extracted to directly isolate the total RNA on the column, as described Table 2 The 49 annotated candidate reference genes ranked in order of their coefficient of variation (spot: location on microarray, ID genbank: accession number, stdev: standard deviation, description: best BLAST hit description as described in GigasDatabase, e value: best BLAST hit e-value, CV: coefficient of variation, max: maximum value, min: minimum value, MFC: maximum fold change, p value: one-way ANOVA P > 0.05) (Continued) by the manufacturer. Absence of trace genomic DNA was confirmed by gel electrophoresis. RNA concentrations were determined using a ND-1000 spectrophotometer (Nanodrop Technologies) at 260 nm, using the conversion factor 1 OD = 40 μg/ml RNA. RNA integrity was verified on an Agilent bioanalyzer using RNA 6000 Nano kits (Agilent Technologies), according to manufacturer's instructions, without consideration for the RNA Integrity Number (RIN) [49]. Indeed, in molluscan RNA, the comigration of the 28S rRNA fragments with the 18S rRNA prevented us from using the RIN, thus only the absence of RNA degradation can be considered [50][51][52]. A detailed description of the validation of RNA integrity in oysters is provided in Additional file 8. Samples were stored at -80°C until use.

Oligonucleotide microarray design and construction
The sequences used to design oligonucleotide probes were obtained from the publicly accessible GigasDatabase version 6, available on the Sigenae website (http://public-contigbrowser.sigenae.org:9090/Crassostrea_gigas/index.html) [21]. In July 2010, the final dataset contained 31,952 contigs assembled from 56,268 sequences selected from cDNA libraries derived from a wide variety of oyster tissues and developmental stages (Additional file 9). Contigs had previously been annotated with BLAST × and added to GigasDatabase [21]. From the resulting 31,952 contigs, 60-mer oligonucleotidic probes were designed using eArray software with the default parameters (Agilent Technologies). After repeat masking, low complexity filtering and similarity searches, 31,918 sequences contained exploitable probes. One probe per sequence was then printed onto a 4 × 44 K custom GE Microarray (Agilent Technologies). This automatic probe design generated 34 duplicated probes from highly similar sequences (Additional file 1). These probes were used to evaluate the reproducibility of hybridization (see results). Probes for negative controls were also printed onto the array in situ and used to test for low background before analysis (see Correction and Normalization). A new assembly (Gigas Database version 8) has been generated since the end of this study. The retrieval of sequences is possible using a Genbank name query, followed by .p.cg.8. All sequences used in this study are also available in Genbank and EMBL and are given in Additional file 9.
Oligonucleotide microarray hybridization RNA Amplification, labeling and hybridization For microarray hybridizations, 200 ng of total RNA was indirectly labeled with Cy3 using the Low Input Quick Amp labeling kit (Agilent Technologies), according to the manufacturer's instructions. Qiagen RNeasy mini spin columns were used for purifying amplified RNA (aRNA) samples. After purification, aRNA amplification and dye incorporation rates were verified using a ND-1000 spectrophotometer (Nanodrop Technologies) and shown to lie between 200 and 500 ng/μL (aRNA concentration) and between 20 and 50 pmol/μg aRNA (dye incorporation). Hybridization was performed using the Agilent Gene expression hybridization kit (5188-5242), as described by the manufacturer, with 1.65 μg of aRNA labeled with Cy3. Tissue samples were randomly hybridized onto 10 different slides, which were subsequently treated with Gene expression wash buffer solution (5188-5327; Agilent Technologies), and Stabilization and Drying solutions (5185-5979; Agilent Technologies). Finally, slides were scanned on an Agilent Technologies G2565AA Microarray Scanner system at 5 μm resolution, using default parameters.

Correction and Normalization
Feature extraction and data normalization were conducted with Agilent Feature Extraction software 6.1. (Agilent Technologies), using the default/recommended normalization method. Initial microarray testing revealed cross hybridization of C. gigas aRNA on the spiked-in probes. Thus Agilent spike-ins were not included in this study and are not recommended for future microarray analysis of C. gigas. The array contained 153 negative control spots with a mean gene expression of 0.84 (variance 19.6). Mean expression level of negative control probes was calculated for each tissue. The maximum mean obtained for each  tissue × negative control probe was then used as a threshold level to determine C. gigas gene expression. According to this, 3,630 (11.4%) C. gigas probes had a gene expression below the background level in all tissue samples. A matrix of gene expression levels was generated, where each row corresponded to a different C. gigas oligonucleotide probe and each column to one tissue sample. The expression level of each contig was then logarithmically transformed, centered around 4 (the relative mean of all contigs within all tissues), and reduced so that relative variations rather than absolute values could be used for interpretation. After normalization, probes that showed gene expression at the background level had log 10 values ranging from 2.15 to 3.15, with a mean of 2.65.

Microarray data analysis
We initially applied a principal component analysis (PCA) using geneANOVA software [53] to assess the internal consistency of different transcriptional data sets and to obtain the proportion of variance for each principal component. Although component 1 had the highest proportion of variance (87.4%), it did not discriminate tissue samples from each other. All samples had exactly the same proportion of variance within this component. PC1 would represent the great majority of genes that do not show tissue-related differential expression. Components 2, 3 and 4 discriminate the tissues and were used to draw a 3D scatter plot (XLstat; Addinsoft) and to organize the 33 tissue transcriptomes along the principal components. Hierarchical clustering using Pearson's correlation was then performed using TMeV [54,55] and the results were compared to the PCA. The values of normalized signal intensities were further analyzed using one-way analysis of variance (ANOVA), with a p-value threshold of less than 0.01 and adjusted with Bonferroni correction (TMeV). Data were further filtered to determine whether their expression was specific to a tissue type. For each gene, the means of normalized expression from tissue replicates were calculated and compared to the means of normalized expression obtained for each of the other tissues.

Quantitative RT-PCR
For quantitative RT-PCR analysis, RNA samples were treated with DNAse I (Promega; 1 U/μg total RNA) for 30 minutes at 37°C according to the manufacturer's instructions. After DNAse I treatment, absence of genomic DNA was confirmed by qPCR on a total RNA sample. Reverse transcription was then carried out using 1 μg of total RNA from each sample, 1 μg of random hexanucleotidic primers (Promega), 0.5 mM dNTPs and 200 U M-MLV reverse transcriptase (Promega) at 37°C for 1 hour in the appropriate buffer. Amplification reactions were performed in 1 × Absolute blue qPCR SYBRGreen fluorescein mix (Thermo scientific) with 5 ng of cDNA template and 300 nM of each primer in a final volume of 15 μL using an iCycler iQ™ thermocycler (Bio-Rad). The comparative threshold cycle (C T ) method was used to quantify copy number of the target gene in the tissues. Primer sequences used in this study are listed in additional file 3. The PCR amplification efficiency (E; E = 10 (-1/slope) ) for each primer pair was determined by linear regression analysis of a dilution series [56]. The specificity of the primer pairs was confirmed by melting curve analysis at the end of each qPCR run and each amplicon was verified by gel electrophoresis. For selected tissue-specific genes, the ΔCt was calculated using arf1 (adp-ribosylation factor 1) as a reference gene and compared to the normalized expression values from the microarray analysis.

Reference genes
Genes expressed in all tissues with a mean normalized intensity that exceeded 4 (the mean of all normalized expression values) were considered as potential reference genes for quantitative RT-PCR. One-way ANOVA with a p-value threshold < 0.05 was then performed and candidate reference genes were defined as genes that were non-differentially expressed. For each gene, the coefficient of variation (CV, the ratio of the standard deviation to the mean) was calculated. CV was used to compare the degree of variation between the genes with different means of expression [57], and the Maximum Fold Change (MFC, ratio of the maximum and minimum values) was used to reflect the minor variations of those candidates within all biological samples [58]. Genes were thus ranked in order of increasing CV to identify the top most stable candidates.
Quantitative RT-PCR was performed to identify the best candidate control genes. The Ct values were determined at a specific threshold (fixed at 500 RLU). Ct values were normalized across samples by calculating the ratio of the Ct value to the mean of the Ct value obtained for all genes for a given sample and logarithmically transformed (base 10). The CV was calculated for each gene using normalized Ct values. Normfinder [59] was used on raw and normalized Ct values to measure stability values and identify the best housekeeping gene candidates. Finally, we used geNorm on raw Ct values to determine the expression stability (M value) [60]. Results obtained from the different methods were compared to validate the selection of ideal reference genes for comparative studies among tissues.

Additional material
Additional file 1: List of 34 duplicated contigs and the gene expression values obtained within each tissue. This table provides for each probe the sequence: Probe sequence, Row: Row on array, Col: Column on array, ProbeName, ContigName: Sigenae accession number,