Transfer of T. nigroviridis from water containing 2.9 mM Ca²⁺ to water containing 0.01 (LowCa) and 10 mM Ca²⁺ (HighCa) resulted, respectively, in a small decrease and increase in blood plasma total [Ca²⁺], both at 2 and 12 h, which was significantly different between the two experimental groups (Figure 1). At the same time the expression level of gill epithelial calcium channel mRNA (ecac, official name trpv6) changed in the opposite direction, i.e., up-regulation in LowCa and down-regulation in HighCa which were statistically significant (P < 0.05) between these experimental groups 12 h after transfer. The expression of CaSR transcripts was also detected at very low levels in the gills of all experimental groups by quantitative real-time (q) PCR (and by SuperSAGE in the LowCa2 h library, one tag with count 1), but no changes in gene expression in response to water [Ca²⁺] were detected (not shown). These results are in agreement with previous studies 10 11 12 13 25 and confirmed that the experimental approach induced a challenge to calcium balance mechanisms in T. nigroviridis gills.

Five gill SuperSAGE libraries from control fish at 2 h (C2 h) and 12 h (C12 h), from LowCa at 2 h (LowCa2 h) and 12 h (LowCa12 h) and from HighCa at 12 h only (HighC12 h) were sequenced in a single half run of 454 pyrosequencing and yielded 302,033 reads (ditag sequences). Stringent quality control of ditag sequences yielded 344,375 26 bp-SuperSAGE mono-tags which were retained for further analysis (Table 1). Each library contained on average of ~69,000 tags of which ~26,500 were unique (unitags, corresponding to unique transcripts; Table 1).

Approximately 70% of the unitags were singletons (tags with 1 count), while the number of unique transcripts decreased for classes with higher abundance (Table 1). Singleton tags were included in the statistical analysis 26 to avoid losing information from low abundance transcripts, but those appearing only in one library (probably derived from sequencing errors) failed the stringent multiple tests to assess differential expression.

In order to eliminate genes potentially responding to handling stress (given the relatively brief duration of experimental exposure), the gill expression profiles of Ca²⁺ challenged libraries were always compared to the gill profile of fish transferred to control water at the corresponding time-point (control libraries), while two types of statistical tests and multiple library comparisons (see methods) were also carried out to increase stringency and to account for time-effects. A total of 1,426 unique tags (transcripts) showing differential expression between control and [Ca²⁺] treatments were identified of which 1,339 were up- or down-regulated at least 2-fold and selected for further analyses.

The annotation of the differentially expressed tags was carried out by BlastN mapping to available DNA sequences for this species, followed by BlastX assignment of these longer (tag-matching) cDNAs or 1000-bp DNA fragments (extracted upstream from the tag-matching sites in the genome) to the highly curated non-redundant Swiss-Prot protein database.

In the BlastN step, carried out under high stringency (26/26 nucleotide matches), 516 tags (38.5%) matched to one or more cDNAs available in GenBank (NCBI cDNAs), 243 tags (18.1%) matched cDNAs predicted from the T.nigroviridis genome (Ensembl cDNAs) and 692 tags (51.7%) matched the genome (Additional file 1), making a total of 869 tags with a DNA match to available T.nigroviridis DNA datasets (Table 2). Three hundred and fifty three tags (26.4%) mapped only to the genome or Ensembl genome predicted cDNAs (Additional file 1) and may represent novel transcripts. It was not possible to calculate the rate of unambiguous tag mapping to unique transcripts, since the available cDNA dataset contained multiple redundant expressed sequence tags (ESTs) and no unigene 27 dataset is available for T.nigroviridis. However, only a reduced number of tags matched more than one location in the genome (28 tags, 4%), suggesting unambiguous annotation for most 26-bp tags.

Four hundred and seventy differentially expressed tags (35.1%) could not be mapped to any available DNA sequence and may represent genes or splice variant transcripts not previously covered by the sequencing of genomic and expressed sequence tag (EST) data. Depending on the DNA dataset, lowering stringency to 24/24 nucleotide matches improved annotation up to 19% (Additional file 2), suggesting that some tags contain nucleotide mismatches to available DNAs, potentially derived from sequencing errors or single-nucleotide polymorphisms (SNPs) in transcript sequences. However, even in these circumstances 60% (279) of the non-annotated tags could not be assigned, supporting the notion that they are not represented in available sequence databases. Therefore the 26/26 (perfect match) stringency rule was maintained in this analysis to ensure a less ambiguous annotation of differentially expressed tags.

In the second step of annotation (BlastX of tag-matching DNAs), the following preference order was used when assigning a protein ID for the tags matching different DNA datasets: the NCBI cDNA BlastX hit was used preferably compared to Ensembl cDNA hits (which are predicted cDNAs lacking verification of structure and expression), and only when no tag-matching cDNAs were available or no significant BlastX hit occurred (significance at expect value, E < 10^-5) was genomic match used. This choice was supported by the lower BlastX annotation rate for genomic 1000-bp fragments (28.2%) compared to NCBI and Ensembl cDNAs (>70%; Additional file 1), perhaps as a result of introns, regions with low sequencing quality or errors in assembly in genomic sequences. Nevertheless, significant BlastX hits for tags matching more than one DNA dataset (223 tags) were found to have the same annotation in 186 cases (83.4%), to belong to the same protein family in 32 cases (14.4%) and to differ only in 5 cases (2.2%). These tags and those for genes selected for qPCR were subjected to careful manual annotation. Detailed annotation results can be found in Additional file 3, which summarizes all data concerning the 1339 differentially expressed tags.

In summary, we could map with high stringency 869 tags (65%), 372 to anonymous DNAs and 497 to curated Swiss-Prot protein entries (Table 2). The number of different genes identified (413) was lower than the number of tags annotated to protein, as 54 genes (13.1%) had multiple tags assigned to the same gene (average 1.2 differentially expressed tags per gene, maximum 7), suggesting they were products of sequence polymorphisms and/or alternative mRNA processing (i.e., produced by alternative splicing or use of alternative polyadenylation sites, which may affect the position of the most 3'-end NlaIII site and thus generate different tags). In addition, 43 (8.6%) of the annotated differentially expressed tags and 37 of the genes (9%) were identified as putatively derived from natural antisense RNA transcripts, as they had an inverted match to a sense cDNA or a direct match to an antisense cDNA.

We first looked at the 10 transcripts that underwent the largest up- or down-regulatory change for each challenge that could be annotated to well characterized genes (Table 3, see Additional file 3 for all differentially expressed transcripts). These included transcripts for a range of calcium-binding proteins or participants in CaBP protein complexes involved in muscle contraction or cytoskeleton organization, such as actn3, mle3, mlrs, tnni2, tpm1, prvb and telt, which were rapidly and strongly up-regulated in gills of fish exposed to water with low [Ca²⁺]. Two genes, entk (involved in proteolysis) and rl11 (translation), were among the top 10 up-regulated transcripts at 12 h for both LowCa and HighCa, while co8a1 (an extracellular matrix constituent) and frim (iron ion oxidation and transport) were among the most down-regulated for both LowCa and HighCa, suggesting that some common mechanisms were activated.

K-means clustering was used to group transcripts according to their patterns of expression between control and low Ca²⁺ water over time (LowCa analysis) and across the three Ca²⁺ availability conditions at 12 h (12 h analysis, to compare the effects of high vs. low Ca²⁺ water challenges).

For each analysis, six clusters (optimal number of K = 6 obtained by Gap statistics) were obtained (Figure 2A and 2B), containing 192 ± 44 (mean ± standard deviation) tags each, which showed clear distinguishable patterns between them. Since most transcripts were common to both analyses, the number of tags common between clusters of the two analyses is also represented in Figure 2C. To understand the biological significance of the identified differentially expressed transcripts, the gene ontology 28 terms (GO) significantly enriched in the whole lists of differential expressed tags or in each of the clusters compared to the normal gill transcriptome (i.e., the annotated SuperSAGE profiles obtained in control water at 2 and/or 12 h) were identified and are listed in Additional files 4 and 5.

The biological process categories most significantly over-represented in the global lists of differentially expressed tags were mainly related to muscle contraction and cytoskeleton: categories "regulation of muscle contraction", "muscle contraction" and "cytokinesis" for both LowCa and 12 h analysis; the "positive/negative regulation of ATPAse activity" and "response to calcium ion" in the LowCa list and the "actin filament-based movement" in the 12 h list. The "phosphocreatine biosynthetic process", involved in cellular energy homeostasis, was enriched in LowCa at 2 but not in the global 12 h list. Among the enriched molecular functions in both analyses were the "structural constituent of muscle" and "motor activity categories", while the "calcium channel regulator activity" was enriched in the LowCa list; finally, enriched cellular components included "sarcomer", "myofibril" and the troponin/myosin complexes.

When detailing the expression patterns of transcript tags whose expression was altered by water with low [Ca²⁺], cluster 6, grouping 223 tags with a strong up-regulation at 2 h but not at 12 h (Figure 2A), had the highest number of enriched categories (273) with the highest significance probabilities (adjusted p value down to 10^-52) and was the main contributor to enrichment in the global list of differentially expressed transcripts (containing 89% of the enriched categories in the global LowCa list, Additional file 4). This suggests a strong coordinated response of functionally related genes in response to the LowCa stimulus at 2 h. Enriched biological processes were mainly related to Ca²⁺ signaling or homeostasis, muscle contraction and cytoskeleton, including all the categories indicated above plus specific categories such as "calcium-mediated signaling", "calcium ion homeostasis" and "actin cytoskeleton reorganization". Transcripts grouped in this cluster and responsible for the enrichment in these categories include actin α skeletal muscle (acts), several myosin (myss, mle3, myl4, mlrs) and troponin forms (tnnc1, tnnc2, tnni1, tnni2, tnni3, tnnt2, tnnt3), α-tropomyosin (tpm1), the sarcoplasmic/endoplasmic reticulum (SR/ER) calcium ATPase 1 (at2a1 or serca1), the sarcomeric mitochondrial creatine kinase (kcrs) and creatine kinase muscle type 1 (kcrm1). Glucose metabolism also appears to have been affected in LowCa, with "gluconeogenesis" enriched in this cluster (e.g. genes f16p2 and g3p), as well as up-regulation of genes involved in glycolysis, such as aldoa, enob and g3p. This suggests an increase in the conversion of non-carbohydrate substrates (e.g. amino acids or lipids) into glucose and glucose oxidation, probably to meet the energetic needs of the acclimation response.

In addition to several molecular functions related to muscle contraction/cytoskeleton (Additional file 4), the "calcium ion binding" category was enriched, with half of the CaBPs identified grouping in this cluster. The "creatine kinase activity" (involved in cellular ATP homeostasis) was also enriched with several creatine kinase genes grouping in the cluster (kcrb, kcrf, kcrm1, kcrs), and enriched cellular components included sarcomer, troponin and myosin complexes, myofibril and cytoskeleton.

Changes in expression of 30% of the differentially expressed transcript tags in cluster 6 appear to be specific to the early-response to LowCa, as they did not meet the criteria for the 12 h differential expression among water types (Figure 2C).

Cluster LowCa-3 contains 106 transcript tags also up-regulated by LowCa but only after 12 h (Figure 2A). Cytoskeleton remodeling appears to be induced also at 12 h, with enrichment in "intermediate filament-based process" and up-regulation of keratin 13 (k1c13) and the intermediate filament protein ON3 (ion3), while "glutathione transferase activity" was also enriched (genes gsta and mgst3), suggesting an increase in detoxification activity. Most transcripts in this cluster (66%, e.g. k1c13 and ion3) fall into cluster 4 of the 12 h analysis (Figure 2C), indicating they were regulated by LowCa but not HighCa at the 12 h exposure, while 29% (e.g. gsta and mgst3) were also regulated by HighCa (cluster 12 h-6).

Cluster LowCa-2 includes transcripts up-regulated by LowCa at both 2 and 12 h (Figure 2A), mainly enriched in categories related to energy generation such as "mitochondrial electron transport, NADH to ubiquinone", "ATP synthesis coupled electron transport", "NADH dehydrogenase (ubiquinone) activity", "mitochondrion", and contains transcripts for several NADH dehydrogenases and NADH-ubiquinone oxidoreductases. Thirty eight percent of the cluster transcripts (including nu2 m and nua4l) showed a general up-regulation by both low and high calcium challenges (cluster 12 h-6, Figure 2).

Interestingly, clusters LowCa-1 and -5 grouped tags with changing expression levels in the control water over time (Figure 2A), but down-regulated by LowCa at 2 or 12 h, respectively. Whether or not these patterns of temporal variation in expression in control water reflect a response to stress, a natural control by circadian rhythms (our sampling times of 2 or 12 h coincided with the beginning and end of the light period, at 8 a.m. and 8 p.m., respectively) or other type of control remains to be determined. Cluster 1 was enriched in the biological processes "mitotic G2 checkpoint" and included genes related to cell division such as ccnd2, and "defense response" (e.g. genes ha2b, hg2a), while cluster 5 was enriched in the processes "response to unfolded protein", "sleep", "response to UV" and "positive regulation of Wnt receptor signaling pathway" (genes hsp71, grp78). Ninety percent of the transcripts in cluster LowCa-5 were down-regulated by both Low and HighCa water at 12 h and fell into clusters 1 or 5 from the 12 h analysis (Figure 2B and 2C). Finally, cluster LowCa-4 contained transcripts down-regulated by LowCa at 2 h and up-regulated at 12 h and was mainly enriched in "proteolysis" (including genes for several proteasome subunits psa1, psa6, psa7, psb8 and prs7; the calpain subunit can2 or the enteropeptidase entk) and "tissue regeneration" (e.g. fibulin 5, fbln5), probably representing genes involved in tissue/cell remodeling during the responses towards acclimation to external Ca²⁺ changes.

The most interesting clusters for the comparison between the effects of low and high water [Ca²⁺] exposure at 12 h (Figure 2B) are clusters 5 (transcripts down-regulated by both water types) and 6 (up-regulated by both), and clusters 2 and 3 (up/down-regulated by HighCa but not LowCa) (Additional file 4). Cluster 5 is enriched in "regulation of striated muscle contraction" (genes tnni1, tnnt3) while cluster 6 is enriched in "collagen catabolic process", containing genes for proteases such as pepd but also one collagen form, co8a1. Cluster 2 includes transcripts only up-regulated by HighCa at 12 h but not by LowCa, 43% of which are also not regulated by LowCa in the 2 h exposure time (Figure 2A, B and 2C). However, no GO category was significantly enriched in this cluster, indicating that it contains a heterogeneous group of non-functionally related genes. Examples of genes only up-regulated by HighCa include ubiquitin (ubiq, proteolysis), the GTP-binding nuclear protein Ran (ran, nucleocytoplasmatic transport) and serine/threonine-protein kinase sgk1 (stress response and negative regulation of apoptosis).

Cluster 3, containing genes down-regulated by HighCa at 12 h but not (or only slightly) by LowCa is significantly enriched in 124 categories (Additional file 4), 66% of which were also enriched in cluster LowCa-6, such as "muscle contraction", "cytokinesis", "positive regulation of ATPAse activity", "actin filament-based movement" and "phosphocreatine biosynthetic process"/"creatine kinase activity". This result indicates a difference in the responses to low or high Ca²⁺, with the same processes that were up-regulated in LowCa at 2 h being down-regulated in HighCa at 12 h (although we do not know if they were regulated in HighCa at 2 h, since this group was not included in the study). Specific categories enriched in this cluster compared to clusters of the LowCa analysis include "protein secretion", "cell motility" or "ephrin receptor signaling pathway" (with down-regulation of the ephrin type-A receptor 2, epha2, and the proto oncogene tyrosine-protein kinase yrk).

To focus on the effects of water Ca²⁺ availability on the gill intracellular signaling machinery and on the gill effector mechanisms, we highlight in Table 4 five selected functional categories and related genes that were identified as differentially expressed through transcript tag annotation. Thirty genes (41 transcripts) encoded different CaBPs and mainly grouped to clusters LowCa-6 (51%) and -4 (24%), indicating that the expression of most identified CaBPs was rapidly up-regulated in response to LowCa water (e.g. calsequestrin, parvalbumin β and several myosin forms), while 19% was down-regulated by HighCa at 12 h (cluster 12 h-3, e.g. calsequestrin). One tag mapping to the important intracellular Ca²⁺-sensing and signal transduction protein calmodulin (CaM) was slightly (1.5-1.7-fold) but significantly (p < 0.05) down-regulated by LowCa at both 2 and 12 h, but is not included in our dataset of differentially expressed transcripts as its expression ratios are below our cut-off value of 2-fold.

Fourteen genes for transcription factors were regulated by Ca²⁺ with a variety of induction profiles (up or down-regulation by Low or HighCa), while 14 additional regulators of transcription were identified, most of which with a rapid response to LowCa (43% in cluster LowCa-6). A large number of transcripts for components of intracellular signaling pathways were also identified (54 differentially expressed transcripts, 47 different genes) with varying expression patterns (Table 4, Additional file 3).

Eighteen genes (23 transcripts) coding for ion transporters/ion homeostasis regulators and other proteins described to influence ionic changes in gills were identified among the differentially expressed transcripts (Table 4), including five tags matching members from the claudin (CLDN) superfamily of transmembrane tight junction proteins. Annotation of these and 63 additional (non-differentially expressed) cldn transcript tags via Swiss-Prot BlastX mainly matched mammalian claudins 3 and 4, but it was not possible to distinguish between fish claudin paralogues, as this gene family has been largely expanded in teleost fishes 29 . Cldn tags were thus annotated by comparison with F. rubripes proteins, allowing the assignment of differential expression to fish cldn3a, cldn28a, cldn32a and cldn8d genes [following the nomenclature established in 29]. Similarly, a transcript for an aquaporin form has been annotated to fish aqp3, using BlastX against fish proteins only.

The same mRNAs used for SuperSAGE were individually analyzed by qPCR for eight genes detected to be differentially expressed by SuperSAGE. Gene selection took into consideration the coverage of a wide range of transcript (tag) abundance and patterns of differential expression (Additional file 6). It also included genes containing a variable number of alternative tags although qPCR primers were not designed to distinguish alternative transcripts. A statistically significant positive correlation was obtained between the SuperSAGE expression (sum of counts of all tags matching each gene) and qPCR expression levels (relative to rps18) for 5 of the genes tested (see Additional file 6). Furthermore, there was a highly significant positive correlation (r = 0.932, p = 3.7 × 10^-11) for the relative expression of the 8 genes relative to control (Figure 3) between qPCR and SuperSAGE, indicating an overall concordance between the two techniques in detecting differential expression.

T. nigroviridis exposed to lower [Ca²⁺] had lower total calcium plasma levels than those exposed to high [Ca²⁺] at both 2 h and 12 h post-transfer (Figure 1), indicating a short-term dependence of total circulating levels on external [Ca²⁺], as observed in other fish species 7 35 36 . These changes are likely the result of the coordinated action of the gill, intestine and kidney, with decreased gill calcium uptake in LowCa water. This is supported by the up-regulation of the (apical side) epithelial Ca²⁺ channel (ecac) mRNA detected by qPCR at 12 h (Figure 1) which suggests that this compensatory effector mechanism was activated in order to increase transepithelial Ca²⁺ transport and restore Ca²⁺ homeostasis. The counteracting mechanism was also observed in HighCa water with ecac mRNA expression down-regulated. The plasma calcium and ecac response to changing water calcium are similar to what has been previously described in other fish species 10 11 12 13 and acted as a positive control in our experiment.

With approximately 70,000 transcripts analyzed per library and assuming around 300,000 transcripts are expressed per cell 26 37 , our analysis provides ~1-fold coverage for transcripts expressed at >0.9 copies per cell and thus the produced profiles are comprehensive of the gill transcriptome. The tag extraction yield was 57% and approximately 70% of the unitags were singletons. These statistics are in agreement with other reports using the same methodology 26 37 .

The rate of annotation was high with 65% of the 1,339 differentially expressed SuperSAGE tags which could be annotated (57% of which matched known proteins). The specificity of tag mapping (only 2-4% tags mapped to multiple genomic locations) was also high and within the level expected for 26-bp tags 38 . The assignment of multiple tags to 13.1% of the genes annotated highlights the capacity of SAGE to detect and discriminate a diversity of mRNA transcripts 22 26 39 40 . However, for 470 tags (35%), some of which among the top 10 most-differentially expressed tags for each challenge, no DNA or protein could be assigned. These results also reflect the capacity of SAGE to identify novel transcripts with potential relevance to the physiological process in study, given the magnitude of their differential expression. These can be isolated and identified [Additional file 7, 41 ] and will be an objective for future studies.

Finally, there was general concordance between SuperSAGE and qPCR in detecting differential expression of genes within a wide range of abundance levels, patterns of differential expression and variable number of alternative tags (Figure 3 and Additional file 6), as found in other SAGE studies using similar analysis 26 27 42 43 44 . The lack of correlation for a minority of genes tested by qPCR may have been caused by their low expression levels and presence of multiple tags in all of them, a fact that has also been previously reported in other studies [e.g., 45 46 47 ].

The main biological processes affected by short-term changes of water [Ca²⁺] were those related to muscle contraction, cytoskeleton organization and cytokinesis, with a rapid up-regulation of several actin, myosin and troponin forms, α-tropomyosin and α-actinin3, among others, 2 h after exposure to low [Ca²⁺]. Some of the same genes were down-regulated 12 h after exposure to high [Ca²⁺].

The concerted altered expression of those genes suggests that cytoskeletal modulation and/or cell proliferation are part of the process of acclimation of T. nigroviridis gills to changed [Ca²⁺]. Furthermore, it is possible that MR cells are key elements of these changes, since several studies have shown modifications to the MR machinery in the gill and skin epithelia in several fish species in response to manipulations of environmental [Ca²⁺] 48 49 50 51 52 . An important role for cytoskeletal actin in the MR cell short-term response to salinity change was previously suggested 53 , since disruption of actin polymerization was able to block the rapid apical crypt closing in response to hypotonic shock, while a thick annular actin ring was detected at the opening of these crypts in MR cells and was suggested to control its opening/closure 54 55 . Although other explanations cannot be excluded, this could indicate that the transcriptomic changes we observed may be, at least to some extent, part of a common response to situations challenging ion homeostasis (such as salinity changes and specific changes in [Ca²⁺], [Na⁺] or both) in which cell number, area and morphology of gill epithelial cells (particularly the specialized ion transporting MR cells) are modified with consequent modulation of ion fluxes [e.g. 46 54 56 ]. While most studies dealt with long-term acclimation periods of days or weeks, there are also reports showing MR morphological/proliferative modifications shortly after the challenge [e.g., 45 46 47 ]. Additionally, other components of the cytoskeleton were also regulated by low or high Ca²⁺ levels, including proteins forming the intermediate filaments (several cytokeratins and ION3) and the microtubule component β5-tubulin, consistent with a proposed role for the MR cell microtubule network in Ca²⁺ uptake in tilapia larvae 57 . Thus, the role of the cytoskeleton and the underlying processes during the early cellular adjustments of MR cells to ion (in particular calcium) and osmotic regulation should be further investigated.

Transcriptional changes of the actin cytoskeleton components could be also linked to changes in the structure and organization of the gills. Muscle cells are present in the gill vasculature and in adductor/abductor muscles that control the angle between gill filaments 8 58 and up-regulation of cytoskeletal actin, troponin and myosin complexes, and myofibril in low calcium may reflect the key role of Ca²⁺ in muscle contraction. Indeed, the concomitant up-regulation of cytoplasmatic CaBP PVB, creatine kinase, Ca²⁺ pump AT2A1 and the enzyme CKM1, all described to play a role in muscle relaxation 59 60 , suggest a coordinated action in muscle function and/or cytoskeleton reorganization in gill tissue.

One of the most enriched GO categories induced by low water Ca²⁺ availability was, as expected, calcium-ion homeostasis. A rapid up-regulation after 2 h exposure to LowCa water was detected for AT2A1 (SERCA1), the Ca²⁺ pump responsible for re-sequestration of Ca²⁺ into the SR/ER that, together with SR/ER Ca²⁺-release channels and cytoplasmatic/ER luminal CABPs, allow the tight regulation of cytosolic [Ca²⁺] in numerous cell types 61 . In addition, several genes for Ca²⁺-buffering CABPs were also up-regulated by LowCa water after 2 h (Table 4), including the ER lumen CaBPs calsequestrin-1 and sarcalumenin, or the cytosolic CaBPs parvalbumin β and ictacalcin 59 61 , revealing for the first time a tight and rapid up-regulation of several elements of the Ca²⁺ homeostasis control machinery at the transcriptional level in fish gills. Whether the detected changes in expression are a general cellular response of fish gill to low calcium and/or are more specifically related to the Ca²⁺ uptake pathway to maintain body calcium balance remains to be determined. However, the current accepted model for transepithelial Ca²⁺ uptake across fish gill MR cells includes a transcellular transport step of Ca²⁺ bound to CaBPs, such as Ca²⁺-calmodulin, parvalbumin and ictacalcin, through the cytoplasm or its sequestration into organelles (mainly the ER) 62 63 . The up-regulation after 2 h exposure to LowCa water of a number of CaBPs previously identified as Ca²⁺-buffers supports their possible role in Ca²⁺ uptake across fish gills, a hypothesis that may be validated by confirming these transcriptional changes specifically in the gill MR epithelial cells responsible for Ca²⁺ uptake.

Very low and unchanging gene expression levels of calcium sensing receptor (CaSR), the membrane receptor which senses extracellular calcium levels, were detected by qPCR and SuperSAGE (not shown). CaSR has been previously localized to the MR cells in fish gill epithelia and shown to activate phospholipase C and mitogen-activated protein kinase (MAPK) signaling in response to external [Ca²⁺] changes 16 . Our results are consistent with reports in which the mRNA and protein expression of CaSR were unaffected by calcium constraint (in water or food) or by salinity changes 4 64 .

SuperSAGE detected the differential expression of a large number of CaBPs (Table 4), including four proteins from the S-100 family (HORN, S10A5, S10AD and S10I), 75% of which were up-regulated after 2 h exposure to low water [Ca²⁺], while CaM was down-regulated at both 2 and 12 h. The immediate early gene response (induction at <2 h) response of these CaBP genes support the notion that they may be part of an intracellular Ca²⁺ sensor network 65 66 in fish gills involved in the initial perception of Ca²⁺ availability and in the consequent activation of signaling pathways and target proteins.

Only a small number of ion channel and transporter genes were found (Table 4) and no differential expression of transporters such as pmca or ncx was detected. This is in agreement with the theory that ECaC is the only (or main) plasma membrane Ca²⁺ transporter regulated by environmental [Ca²⁺] in gills, while the steady state expressions of PMCA and NCX may accommodate the needs for transepithelial Ca²⁺ transport under different [Ca²⁺] 10 . However, it is also possible that the exchange activities of these and other transporters (and not their mRNA expression levels) are modulated by [Ca²⁺] through Ca²⁺ sensing/signaling mechanisms, such as the above mentioned activation by CaBPs or cytoskeleton reorganization. For instance, the main factors regulating PMCA activity in eukaryotes appear to be Ca²⁺-bound CaM and phosphorylation by protein kinases A/C, but recent evidence in human erythrocyte membranes also point to effects of the interaction with actin cytoskeleton elements on PMCA activity 67 68 .

Additional ion homeostasis related genes were found to be significantly regulated by water Ca²⁺ levels (Table 4), including four members of the large multi-gene family of tight-junction proteins claudins (CLDNs), which have been implicated in the maintenance of hydromineral balance across osmoregulatory epithelia of euryhaline fishes 33 69 70 71 72 73 . It has been suggested that a decline in their expression may be related to the reshaping of the gill epithelia of euryhaline fish upon SW-acclimation and account for its "leakier" (higher ion permeability) properties. For example, both cldn8d and cldn3a expression were reduced in a salinity-dependent manner during T. nigroviridis acclimation for 2 weeks to FW and SW 69 70 . In tilapia and flounder gills CLDN3 proteins were reduced 4-8 days after FW-SW transfer and increased after SW-FW transfer 71 73 . In our study a decrease in cldn8d gill mRNA expression was detected after 12 h exposure to low and high [Ca²⁺] water compared to control water, while the expression of cldn3a was rapidly up-regulated by low [Ca²⁺] at 2 h and unchanged in both low and high [Ca²⁺] water at 12 h post-transfer. As for cldn28a and 28b no overall effects of FW-SW transfers were observed in salmon or tilapia 71 72 , but an up-regulation of cldn28a was detected in SW-FW tilapia 1-7 days post-transfer 71 . In our study cldn28b was also unaffected by water [Ca²⁺] and the slight down-regulation of cldn28a detected by SuperSAGE in LowCa at 2 h was not confirmed by qPCR. Finally, differential expression was detected for two alternative tags annotated to the fish cldn32a isoform, one of the fish cldn genes whose orthologs appear to have been lost in the mammalian lineage 29 , and these results provide the first evidence for regulation of its gill mRNA levels in response to water ion availability. Altogether, these observations suggest short-term alterations in paracellular epithelial permeability, the details of which are not very well known.

Two tags annotated to α carbonic anhydrases (Table 4) were up-regulated in LowCa but not HighCa, consistent with the reported up-regulation of CAH1 and CAH2 mRNA in zebra fish gills after 5-6 days exposure to soft-water 12 . Cytoplasmic carbonic anhydrases in gills may assist in apical ion uptake by providing HCO₃ ^- for apical Cl^-/HCO₃ ^- exchanger or protons for apical Na⁺/H⁺ exchanger 12 . Their up-regulation suggests an impact of Ca²⁺ availability in the gill uptake of other ions.

Finally, one tag for the water pore aqp3, a protein that localizes in the basolateral membrane of MR cells and expresses at lower level in SW than FW 74 75 , was down-regulated by both HighCa and LowCa challenges. This could indicate lower water permeability perhaps to minimize the potential increase in paracellular permeability evidenced in the changes in CLDNs gene expression.

From the statistical analysis of differential expression, clustering and GO enrichment analysis, it became evident that the most significant changes occurred in the transfer to low [Ca²⁺] at 2 h. Whether there was an equivalent pattern of response to high [Ca²⁺] water it was not possible to confirm because it was only analyzed at 12 h. However, the analysis of this SAGE library allowed the comparison of the impacts of low vs. high water [Ca²⁺] on the gill transcriptome, summarized in Figure 4. It turned out that although specific biological processes may be regulated, the majority of the affected processes and genes respond to both challenges deviating from normal conditions of water Ca²⁺ availability, although possibly with different time-frames, magnitudes and signal. For example, a large proportion of genes identified as up-regulated by LowCa water at 2 h were down-regulated by HighCa water at 12 h, including those more highly represented in GO terms, related to actin cytoskeleton, muscle contraction and citokinesis.

All animal maintenance and handling procedures were carried out in compliance with the recommendations of the Association of Animal Behavior 76 and national legislation. Wild green spotted puffer fish, originally imported from Thailand, were purchased from a local pet shop and maintained in 200 L closed circuit aquaria in brackish water (prepared from mixing fresh water and natural sea water to 10 ppt salinity) at 26 ± 1°C, pH 7.5 ± 0.5 and natural photoperiod. They were fed once a day with frozen mollusks. After >2 months of acclimatization, in April 2007, fish were divided in three groups of six (body weight 4.4 ± 1.6 g) and transferred to 25 L aquaria containing artificial brackish water (ABW) with 10 ppt salinity and different Ca²⁺ concentrations: 2.9 mM (control group), 0.01 mM (as in ion-poor fresh water, LowCa group) or 10 mM Ca²⁺ (as in seawater, HighCa group). ABW (10 ppt) was prepared from deionized water and salts (Sigma-Aldrich, Madrid, Spain) with a composition based on that for artificial seawater: 131.4 mM NaCl; 2.9 mM KCl; 7.1 mM MgCl₂; 7.1 mM MgSO₄ and 0.01, 2.9 or 10 mM CaCl₂, for the Control, LowCa or HighCa ABW, respectively. LowCa ABW was supplemented with 5.7 mM choline chloride (C₅H₁₄NOCl) to maintain osmolarity relatively to other ABW types, and pH adjusted to 7.5 with 1 mM NaOH in all water types. Osmolality, measured with a vapor pressure osmometer (Fiske One-Ten Osmometer, Fiske, VT, USA), was 275-300 mOsm/kg; total ammonia (NH₃/NH₄ ⁺) and nitrites (NO₂ ^-) in water, measured using TetraTest kits(Tetra Werke, Melle, Germany), were at recommended levels and temperature was kept at 26°C. No significant changes in water pH, temperature and [Ca²⁺] occurred during the time course of the experiment.

Three fish were removed from each tank after 2 h and 12 h, anesthetized (0.02%-phenoxyethanol, Sigma-Aldrich), blood sampled and fish decapitated. Blood was sampled from the cardiac dorsal vessel using heparinised glass capillary tubes and centrifuged for 5 min at 10,000 g for recovery of plasma. Complete gills (arches and filaments) were removed and snap frozen in liquid nitrogen until RNA extraction. The experiment was repeated twice in the same week (total n = 5-6 fish/group). Plasma total calcium concentration (mM) was measured using a colorimetric assay (o-Cresolphtalein, calcium kit, Spinreact, Girona, Spain) in duplicates, and differences between groups were analyzed by two-way ANOVA followed by Tukey' test at a statistical significance level of p < 0.05. The statistical software used was SigmaStat v.3.00 software (SPSS Inc, Chicago, USA).

Total RNA was extracted from individual frozen gills with TRI Reagent (Sigma-Aldrich) and analyzed for integrity and purity by agarose gel electrophoresis and spectrophotometry, respectively. Poly(A)⁺-RNA was purified from pooled total RNA (50 μg from each fish) for each experimental group (n = 5-6 per group) using the Oligotex mRNA Mini Kit (Qiagen, Hilden, Germany) and analyzed by denaturing agarose gel electrophoresis. Five SuperSAGE libraries were constructed from 1 μg of gill poly(A)⁺-RNA each, corresponding to the groups sampled at 2 h (LowCa and control) and 12 h (LowCa, HighCa and control), following the protocol described by Matsumura et al. 23 77 , but instead of ditag concatenation, cloning and sequencing, 200 ng of amplified ditags from each library (evaluated using a 2100 Bioanalyzer, Agilent Technologies, Palo Alto, California) were directly sequenced by massively parallel pyrosequencing using a GS20 sequencer (1/2 run; 454 Life Sciences, Branford, CT, USA) at the Max Planck Institute for Molecular Genetics, Berlin.

SuperSAGE tags were extracted from sequenced ditags using the "SuperSAGE_tag_extract_pipe" program pipeline, which included stringent quality-control steps for removal of duplicated ditags, ditags containing incomplete library-identifying linkers and ditags with contaminating linker sequences in interior or ambiguous bases 24 . Comparison of tag frequencies between libraries was carried out by the "SuperSAGE_tag_freq_comp" script 24 . The SuperSAGE tags (transcripts) were deposited in NCBI's Gene Expression Omnibus 78 and are accessible through GEO Series accession number GSE19854.

The relative frequency of each tag across different SuperSAGE libraries was compared by combinations of G-tests for multiple library comparison 41 79 and Z-tests for pairwise comparisons 79 80 using the G_test and SAGEstat software (supplied by the authors) and significance levels (p) set at 0.05. Four different G-tests were performed using different subsets of libraries and decision rules 41 79 : 1) G-Test1 = Control vs. LowCa, 2) G-Test2 = pooled Control and LowCa at 2 h vs. pooled Control and LowCa at 12 h, 3) G-Test3 = Control vs. LowCa and Control vs. HighCa at 12 h, and 4) G-Test4 = Control vs. pooled LowCa and HighCa at 12 h. As a result, four different classes of differentially expressed tags were identified: 1) tags "differentially expressed between control and LowCa water independently of exposure time" (global effect) passed all decision rules in G-Test1 and the pairwise Z-test control vs. LowCa (p < 0.05); 2) tags "regulated by LowCa water over time" (time-dependent effect) passed G_intrinsic in both G-Test1 and G-Test2, failed any of the other decision rules in G-Test1, and the control vs. LowCa effect at 2 or 12 h was confirmed by passing Z-test and failing any of the homogeneity rules in G-Test2; 3) tags "differentially expressed between control and HighCa at 12 h" passed all decision rules in G-Test3 and the Z-test for this comparison; 4) "tags with a global effect of altered water [Ca²⁺] at 12 h" passed all decision rules in G-Test4 and in the Z-test of control vs. Low or HighCa at 12 h.

In order to avoid losing biological information from low abundance transcripts by the exclusion of tags with zero counts before tag relative frequency calculation, the G-test software estimated a "zero-substitution value" of 0.345 ± 0.03 (mean ± standard error) 41 79 . Only tags whose differential expression passed statistical significance in all defined G- and Z-tests and with expression ratio R (tag counts of water [Ca²⁺]-altered libraries/tag counts of control libraries at each time point) ≥2-fold were subjected to further analysis. When tag count was zero the "zero substitution value" 0.345 was assigned before calculation of R.

The annotation of SAGE tags (assignment to the corresponding mRNA and protein) is frequently made by interrogating databases of virtual tags extracted from unigene entries 26 27 42 43 44 . However, since no virtual tag database or unigene entries are available for T. nigroviridis an alternative 2-step annotation procedure was used. In the first step, SuperSAGE tags were mapped to T. nigroviridis genomic/cDNA sequences (available July 2008), running the stand alone Formatdb and BlastN scripts (blast-2.2.14-ia32-linux.tar.gz package available at NCBI 81 ) against each of three DNA datasets: 1) all T. nigroviridis mRNA GenBank 82 ("NCBI cDNAs" dataset; 99,852 GenBank sequences); 2) transcripts resulting from known, novel and pseudo gene predictions from the T. nigroviridis genome 83 ("Ensembl cDNAs" dataset, 23,289 sequences); and 3) the latest assembly of the T. nigroviridis genome 84 ("Tetraodon genome" dataset, 21 chromosomes). The "SuperSAGE_tag_BLAST" suite of programs 24 were used to parse BlastN results under stringent conditions (only perfect tag-to-DNA matches, with 26/26 identical nucleotides, were accepted) and to extract DNA fragments 1000-bp upstream from tag-genomic matches.

In the second step, top BlastN cDNA/genomic matches for each tag were compared to the Swiss-Prot protein database (downloaded March 2009 85 ) using stand-alone BlastX. To assign protein identities to each tag, the expect value (E) was set at <10^-5 and the hierarchical preference order for BlastX matches was: NCBI cDNAs > Ensembl cDNAs > Tetraodon genome 1000-bp fragments. This process was verified manually in cases of multiple BlastX matches using queries from different DNA datasets and for the genes analyzed by qPCR. The alignment directions of BlastN (tag to DNA) and BlastX (DNA to Swiss-Prot protein) were compared to identify tags potentially derived from antisense transcripts (those presenting an inverted match to a sense cDNA or a direct match to an antisense cDNA). Tags matching members from the mammalian claudin family of proteins were assigned to their fish claudin isoform by stand-alone BlastX against the 56 Takifugu rubripes claudin proteins (GenBank accession numbers AAT64039-AAT64094).

The 1,339 tags identified as differentially expressed were grouped according to their expression patterns among different libraries by K-means clustering analysis, using the SAGE data analysis tool and the TransChisq algorithm 86 . Two clustering analyses were performed: the "LowCa analysis" clustered 1,208 tags according to their pattern of expression among control and LowCa libraries (all differentially expressed tags excluding those with significant effect only for HighCa challenge and those with zero tag count in all the control and LowCa libraries), while the "12 h analysis" clustered 1,098 tags (all differentially expressed tags except those with zero tag count in all the 12 h libraries) according to their pattern of expression among 12 h libraries (control, LowCa and HighCa). The optimal number of clusters K was obtained by a combination of the TransChisq algorithm and Gap statistics run in MatLab (MathWorks Ltd., Cambridge, UK) as described in 86 .

Gene ontology mapping was performed on differentially expressed tags by retrieving the GO terms 28 associated with all significant BlastX hits (E value < 10^-10) of each tag-matching DNA against the Swiss-Prot protein database. GO enrichment was then determined for each list of annotated differentially expressed tags (all LowCa-regulated tags, each of the six LowCa clusters, all 12 h-regulated tags or each of the six 12 h clusters) using a one-tailed proportion test 87 at a false discovery rate adjusted p-value of 0.01, between the number of tags containing a given GO in each list compared to the number of tags containing that GO in the normal gill transcriptome at the corresponding time point (all tags contained in control libraries at 2 and 12 h, for the LowCa analysis, or at 12 h only for the 12 h analysis).

Eight differentially expressed genes identified by SuperSAGE (ecac, cldn3a, cldn28a, cldn28b, prvb, ileu, at2a1, hbb and ckm1) were selected and analyzed by quantitative real-time PCR (qPCR) in the same individual RNAs used for SuperSAGE library construction. Specific primers for each gene were designed based on available genomic or cDNA sequences and on the alignment with related genes to avoid cross-amplification, using Beacon Design software (Premier Biosoft Int., Palo Alto, CA). Primer sequences, amplicon sizes and optimized temperatures used in qPCR for each pair are shown in Additional file 8. For cDNA synthesis, 4 μg total RNA was treated with DNA-free Kit (Ambion, UK) and cDNA synthesis was carried out using 500 ng of DNAse-treated total RNA, 200 ng of random hexamers (GE Healthcare, Little Chalfont, UK), 40U of MMLV reverse transcriptase (RT) (Promega, Southampton, UK) and 5U of RNAguard RNase inhibitor (GE Healthcare) in a total volume of 20 μl. Quantification was performed in duplicate reactions using SYBRgreen chemistry (Power SYBR^® Green PCR Master Mix, Applied Biosystems, UK) and the relative standard curve method 88 , using a Bio-Rad iClycler iQ5 qPCR thermocycler and software (Bio-Rad Laboratories). PCR cycling conditions were 10 min at 95°C, followed by 55 cycles of 10 sec at 95°C, 20 sec at 60°C and 30 seconds at 72°C. A final melting curve was carried out between 60 and 98°C for all genes and showed single product/dissociation curves. All amplicons were sequenced to confirm specificity.

Standard curves relating initial template quantity to amplification cycle were generated using serial dilutions of linearized plasmid DNA containing the gene of interest, and efficiency ranged between 75-100% with R² > 0.99. No cross-amplification was detected between closely related genes (cldn28a and cldn28b) by testing the PCR over the heterologous template. Absence of DNA contamination was confirmed by testing the amplification of all genes on a sample that did not receive reverse transcriptase. The expression of three candidate reference genes (18s, g3p and rps18) was quantified in the same conditions and no statistically significant differences were found between experimental groups for any of them. rps18 was chosen as endogenous reference gene to normalize qPCR data as it had the lowest variation and a level of expression of the same order of magnitude to target genes. Statistical significance of relative gene expression between groups was analyzed by two-way ANOVA. Pearson correlations between the qPCR data (average expression for each experimental group) and the cumulative signal obtained by SuperSAGE (sum of tag counts annotated to the same gene) were calculated for each gene after log2 transformation of both variables. A Pearson correlation was also calculated between qPCR and SuperSAGE log2 transformed gene expression fold changes between treated groups and control at each time point. The significance level was 5%.