Skip to main content
  • Research article
  • Open access
  • Published:

Gill transcriptome response to changes in environmental calcium in the green spotted puffer fish



Calcium ion is tightly regulated in body fluids and for euryhaline fish, which are exposed to rapid changes in environmental [Ca2+], homeostasis is especially challenging. The gill is the main organ of active calcium uptake and therefore plays a crucial role in the maintenance of calcium ion homeostasis. To study the molecular basis of the short-term responses to changing calcium availability, the whole gill transcriptome obtained by Super Serial Analysis of Gene Expression (SuperSAGE) of the euryhaline teleost green spotted puffer fish, Tetraodon nigroviridis, exposed to water with altered [Ca2+] was analysed.


Transfer of T. nigroviridis from 10 ppt water salinity containing 2.9 mM Ca2+ to high (10 mM Ca2+ ) and low (0.01 mM Ca2+) calcium water of similar salinity for 2-12 h resulted in 1,339 differentially expressed SuperSAGE tags (26-bp transcript identifiers) in gills. Of these 869 tags (65%) were mapped to T. nigroviridis cDNAs or genomic DNA and 497 (57%) were assigned to known proteins. Thirteen percent of the genes matched multiple tags indicating alternative RNA transcripts. The main enriched gene ontology groups belong to Ca2+ signaling/homeostasis but also muscle contraction, cytoskeleton, energy production/homeostasis and tissue remodeling. K-means clustering identified co-expressed transcripts with distinct patterns in response to water [Ca2+] and exposure time.


The generated transcript expression patterns provide a framework of novel water calcium-responsive genes in the gill during the initial response after transfer to different [Ca2+]. This molecular response entails initial perception of alterations, activation of signaling networks and effectors and suggests active remodeling of cytoskeletal proteins during the initial acclimation process. Genes related to energy production and energy homeostasis are also up-regulated, probably reflecting the increased energetic needs of the acclimation response. This study is the first genome-wide transcriptome analysis of fish gills and is an important resource for future research on the short-term mechanisms involved in the gill acclimation responses to environmental Ca2+ changes and osmoregulation.


Calcium (Ca2+) is a major component of the skeleton in vertebrates and the exoskeleton of some invertebrates (e.g. in bones, scales or shellfish shells) and plays key roles in a wide range of physiological processes, such as muscular contraction, modulation of permeability and excitability of plasma membranes, nerve signal transduction or intracellular signaling [1].

Terrestrial vertebrates obtain calcium mainly through their diet and whole-body calcium homeostasis is mainly achieved by intestinal absorption and kidney reabsorption [2]. The range of calcium in water can range from trace levels in soft fresh water to 10 mM in seawater and fish that live in those waters are able to retrieve most or all the calcium they need from the surrounding water. Furthermore, many species (e.g. estuarine) have the capacity to rapidly adjust to a wide range of environmental calcium concentrations [3]. The ability of teleost fishes to maintain circulating calcium [Ca2+] within narrow limits is well documented [46] and is mediated by rapid rate changes of water and calcium exchange in gill, intestine, kidney and skin epithelia [reviewed in [3]].

The gill epithelium is the main site of active Ca2+ uptake from the water at low [Ca2+], while at higher [Ca2+], such as in seawater (SW), the intestine acquires also an important role [3, 7]. The gill, together with the kidney, also appears to play a role in Ca2+ excretion but little is known about this process [3, 8]. Branchial Ca2+ uptake appears to occur mainly through specialized mitochondrion-rich cells (MR cells) or "chloride cells" following a 3-step process similar to that proposed for Ca2+ reabsorption in the mammalian kidney: passive entry of Ca2+ through apical epithelial Ca2+ channels (e.g. ECaC); transcellular transport in the cytoplasm or sequestration in organelles through binding to Ca2+-binding proteins (CaBPs); and active basolateral extrusion to the blood via a plasma membrane Ca2+-ATPase (PMCA) or a Na+/Ca2+-exchanger (NCX) [3, 9]. When fish are exposed to low or high calcium environments, they respond by modifying gill Ca2+ uptake; altering the number, area and morphology of MR cells, gene expression of epithelial calcium channel (ecac) and other osmoregulatory mediators [1013]. However, there is still relatively little knowledge of the gill molecular machinery involved in the adjustments to changes in [Ca2+], in particular the early events.

Finally, the internal and external Ca2+ sensing mechanisms (e.g. via a membrane calcium-sensing receptor, CaSR), intracellular signaling and the endocrine control of Ca2+ balance in fish, which differs from that of terrestrial animals, have been the subject of several physiological studies and recent reviews [3, 1417], but much remains to be discovered about their mechanisms of action and target genes in Ca2+-transporting epithelia.

This study sought to detail the cellular mechanisms and determine the genes involved in the rapid perception, signal transduction and effector responses of teleost fish gills to changes in environmental Ca2+. To this end, we analyzed the gill transcriptome of the green spotted puffer fish (Tetraodon nigroviridis) soon after transfer from brackish water (10 ppt salinity, 2.9 mM Ca2+) to brackish water containing 0.01 and 10 mM Ca2+. The [Ca2+] tested cover the range teleost fish are generally exposed to in nature: from 0.01 mM in ion-poor fresh water to 10 mM in sea water [3]. The choice of time points analyzed after transfer (2 and 12 h) sought to capture the early transcriptomic responses of the acclimation process. T. nigroviridis was chosen as the experimental model because of two main reasons: 1) it is euryhaline and able to tolerate wide variations in water [Ca2+] (can be found in freshwater (FW) streams and rivers, estuaries and coastal waters [18] and in experimental conditions tolerate direct transfers from FW to SW and vice-versa [19]), and 2) its genome and many cDNAs have been sequenced [20], facilitating gene identification from sequenced data. The gill transcriptomic responses were obtained using a variation of serial analysis of gene expression (SAGE) [21] designated SuperSAGE, a global and quantitative transcriptomic technique based on high throughput sequencing and counting of small transcript-identifier tags, 26-bp instead of 14-bp in SAGE. The larger SuperSAGE tags results in improvements in tag-to-gene annotation (gene identification) and characterization of identified tags [22, 23], and has been adapted to massively parallel pyrosequencing using a bar-code system which simplifies the procedure and permits higher throughput [24].

As a result of transfer of T. nigroviridis from water containing 2.9 mM Ca2+ to water containing 0.01 and 10 mM Ca2+ 1,339 differentially expressed gill transcripts were identified and clustered into sub-groups of genes with distinct patterns of expression. Gene ontology (GO) enrichment analysis of the annotated tags allowed the identification of biological processes, functions and cell compartments that were most affected. Relative expression changes in 8 genes were also compared by quantitative PCR (qPCR). The data obtained in this study provide the first comprehensive catalogue of genes implicated in the rapid acclimation responses of fish gill cells to changes in water Ca2+ availability.


Total plasma calcium and Ca2+ channel gene expression

Transfer of T. nigroviridis from water containing 2.9 mM Ca2+ to water containing 0.01 (LowCa) and 10 mM Ca2+ (HighCa) resulted, respectively, in a small decrease and increase in blood plasma total [Ca2+], 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 [Ca2+] were detected (not shown). These results are in agreement with previous studies [1013, 25] and confirmed that the experimental approach induced a challenge to calcium balance mechanisms in T. nigroviridis gills.

Figure 1
figure 1

Blood plasma total calcium and ecac gill mRNA expression upon exposure to different water [Ca2+]. Each bar is the mean ± S.E.M of (A) the plasma total calcium levels (mM) and (B) the relative expression of the epithelial calcium channel mRNA in the gills (ecac expression quantified by qPCR and normalized to the reference gene rps18) of T. nigroviridis individuals (n = 5-6/group) after 2 or 12 h transfer to 10 ppt artificial water containing 2.9 mM Ca2+ (control water, black bars), 10 mM Ca2+ (HighCa water, light grey bars) or 0.01 mM Ca2+ (LowCa water, dark grey bars). Different letters indicate statistically significant differences (p < 0.05) between groups, evaluated by two-way ANOVA.

Effects on global gene expression

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).

Table 1 Summary of tag extraction from the T. nigroviridis gill SuperSAGE libraries

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 Ca2+ 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 [Ca2+] treatments were identified of which 1,339 were up- or down-regulated at least 2-fold and selected for further analyses.

Tag annotation

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.

Table 2 Summary of tag annotation of the 1,339 differentially expressed 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 Nla III 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.

Most up- or down-regulated tags

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 [Ca2+]. 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.

Table 3 Top 10 up- and down-regulated annotated tags for each challenge

Clustering of expression patterns and enriched gene ontology

K-means clustering was used to group transcripts according to their patterns of expression between control and low Ca2+ water over time (LowCa analysis) and across the three Ca2+ availability conditions at 12 h (12 h analysis, to compare the effects of high vs. low Ca2+ 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.

Figure 2
figure 2

Results from K -means clustering of differentially expressed SuperSAGE tags according to their expression patterns. Pattern of expression of all tags for each cluster (mean ± S.E.M of expression, normalized by making the sum of tag counts for all the analyzed libraries equal to 100%), (A) for the LowCa analysis (effects of transfer to 0.01 mM Ca2+ at 2 or 12 h), (B) for the 12 h analysis, describing tag expression at 12 h exposure to LowCa (0.01 mM Ca2+), control (2.9 mM Ca2+) or HighCa (10 mM Ca2+) water. The number of tags grouped in each cluster is shown in each graph. C- Bubble graph showing the relation between clusters of the two analyses. Bubble area is proportional to the number of tags in common by each pair of cluster number (LowCa and 12 h analysis). Tags falling in clusters numbered "zero" were those not analyzed (due to failure of inclusion criteria, see methods) for a given analysis (LowCa or 12 h). Solid line - control; dashed line - LowCa.

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.

Enriched categories in the "LowCa analysis"

When detailing the expression patterns of transcript tags whose expression was altered by water with low [Ca2+], 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 Ca2+ 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 Ca2+ changes.

Enriched categories in the "12 h analysis"

The most interesting clusters for the comparison between the effects of low and high water [Ca2+] 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 Ca2+, 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).

Differential transcript expression of proteins from selected functional groups

To focus on the effects of water Ca2+ 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 Ca2+-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.

Table 4 Summary list of proteins/genes belonging to selected functional groups to which differentially expressed tags were annotated

Fourteen genes for transcription factors were regulated by Ca2+ 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.

Quantitative PCR of differentially expressed genes

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.

Figure 3
figure 3

Relationship between relative change of gene expression measured by qPCR and SuperSAGE. Relative change of gene expression (fold change) between Ca2+-challenged (Low or HighCa water) and control gills for the same response period (2 or 12 h) measured by qPCR (target gene/rps18) and SuperSAGE (sum of tag counts annotating to the same gene). r is the Pearson correlation coefficient.


To our knowledge, this is the first genome-wide transcriptome profile of gill tissue to water Ca2+ availability, and the deepest gill transcriptome study to date, identifying differential expression in 1,339 out of 79,367 different transcripts analyzed. Other studies have analyzed gill transcriptomic changes in response to different challenges such as osmotic stress, heat stress, hypoxia or infection by parasites [3033], but the depth of analysis was much lower (maximum 16,000 different cDNA clones analyzed) and none was concerned with responses to calcium per se.

Three notes of caution should be considered in the interpretation of these results. First, the gill is a heterogeneous, multifunctional organ composed by several types of cells and tissues (epithelial, neural, muscular/vascular and bony/cartilaginous) [8] and the expression profiles produced can result from any calcium responsive cell type. Confirmation of the hypotheses generated by these data as being specific genes to Ca2+ uptake/homeostasis requires confirmation by specific experiments, with focus on the gill epithelial MR cells, considered the primary site of active Ca2+ uptake in this tissue and where most significant changes are expected to occur [3, 7]. Second, gene ontology enrichment analysis is mainly based on the annotation and comparison of lists of differentially expressed genes with data that is biased towards a few (mainly mammalian) model organisms and the involvement of a particular protein in biological processes in fish may differ from that in those species [34]. Third, in common with other transcriptomic studies the produced gene sets are based on altered levels of mRNA and assume parallel, but not confirmed, changes in protein abundance [34].

Because T. nigroviridis does not normally live in waters containing 0.01 mM Ca2+, it could be questioned whether some of the transcriptomic changes reflect acute stress induced by the low [Ca2+] itself. This is unlikely since the "response to stress" biological function category was not enriched among the differentially expressed genes. Furthermore, of 41 proteins included in "response to stress" in the global data set only five (Serine/threonine-protein kinase Sgk1, Serine/threonine-protein kinase Sgk3, Heat shock cognate 71 kDa protein, Heat shock 70 kDa protein 1 and two isoforms of telethonin were differentially expressed, each in a different cluster and equally represented as up-regulated in LowCa at 2 h and 12 h and HighCa at 12 h (additional file 3).

Plasma total calcium and gill expression of ecac

T. nigroviridis exposed to lower [Ca2+] had lower total calcium plasma levels than those exposed to high [Ca2+] at both 2 h and 12 h post-transfer (Figure 1), indicating a short-term dependence of total circulating levels on external [Ca2+], 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 Ca2+ 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 Ca2+ transport and restore Ca2+ 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 [1013] and acted as a positive control in our experiment.

Global Analysis of SAGE data

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, 4244]. 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.,[4547]].

Muscle contraction, cytoskeleton organization and cytokinesis

The main biological processes affected by short-term changes of water [Ca2+] 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 [Ca2+]. Some of the same genes were down-regulated 12 h after exposure to high [Ca2+].

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 [Ca2+]. 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 [Ca2+] [4852]. 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 [Ca2+], [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., [4547]]. Additionally, other components of the cytoskeleton were also regulated by low or high Ca2+ 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 Ca2+ 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 Ca2+ in muscle contraction. Indeed, the concomitant up-regulation of cytoplasmatic CaBP PVB, creatine kinase, Ca2+ 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.

Calcium ion homeostasis

One of the most enriched GO categories induced by low water Ca2+ availability was, as expected, calcium-ion homeostasis. A rapid up-regulation after 2 h exposure to LowCa water was detected for AT2A1 (SERCA1), the Ca2+ pump responsible for re-sequestration of Ca2+ into the SR/ER that, together with SR/ER Ca2+-release channels and cytoplasmatic/ER luminal CABPs, allow the tight regulation of cytosolic [Ca2+] in numerous cell types [61]. In addition, several genes for Ca2+-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 Ca2+ 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 Ca2+ uptake pathway to maintain body calcium balance remains to be determined. However, the current accepted model for transepithelial Ca2+ uptake across fish gill MR cells includes a transcellular transport step of Ca2+ bound to CaBPs, such as Ca2+-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 Ca2+-buffers supports their possible role in Ca2+ uptake across fish gills, a hypothesis that may be validated by confirming these transcriptional changes specifically in the gill MR epithelial cells responsible for Ca2+ uptake.

Ca2+ sensing and signal transduction

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 [Ca2+] 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 [Ca2+], 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 Ca2+ sensor network [65, 66] in fish gills involved in the initial perception of Ca2+ availability and in the consequent activation of signaling pathways and target proteins.

Ion and water homeostasis

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 Ca2+ transporter regulated by environmental [Ca2+] in gills, while the steady state expressions of PMCA and NCX may accommodate the needs for transepithelial Ca2+ transport under different [Ca2+] [10]. However, it is also possible that the exchange activities of these and other transporters (and not their mRNA expression levels) are modulated by [Ca2+] through Ca2+ 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 Ca2+-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 Ca2+ 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, 6973]. 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 [Ca2+] water compared to control water, while the expression of cldn3a was rapidly up-regulated by low [Ca2+] at 2 h and unchanged in both low and high [Ca2+] 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 [Ca2+] 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 HCO3- for apical Cl-/HCO3- exchanger or protons for apical Na+/H+ exchanger [12]. Their up-regulation suggests an impact of Ca2+ 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.

Responses to low versus high water [Ca2+]

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 [Ca2+] at 2 h. Whether there was an equivalent pattern of response to high [Ca2+] 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 [Ca2+] 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 Ca2+ 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.

Figure 4
figure 4

Main biological processes affected by transfer to water at different Ca2+ concentrations in Tetraodon gills. Each triangle represent the differential expression of genes contained in the most represented biological processes (represented in the lower panel and grouped into broad categories in the upper panel) between LowCa (0.01 mM Ca2+) or HighCa (10 mM Ca2+) water and control water (2.9 mM Ca2+, blue triangle). Triangle width is inversely proportional to the adjusted p values obtained for biological process GO categories in our GO enrichment/clustering analysis (wider for most enriched processes); shown categories were selected among the most enriched GOs (adjusted p values > 10-4).


We have used SuperSAGE to carry out a comprehensive analysis of the gill transcriptome during the rapid responses to changes in environmental Ca2+, which identified 1,339 differentially expressed transcripts. The generated transcript expression patterns provide a framework of water calcium-responsive genes in the gill during the initial response after transfer to different [Ca2+]. This molecular response entails initial perception of alterations, activation of signaling networks and effectors and suggests active remodeling of cytoskeletal proteins during the initial acclimation process (Figure 4). It is also possible that some alterations of gene expression may represent disruption of calcium sensitive pathways as a result of the transfer shock which we cannot distinguish at present from those of acclimation. Nevertheless, this data allows the generation of new hypotheses about the mechanisms of acclimation to environmental calcium changes, which can be tested in specific experiments (including localization studies, expression analysis at different time-frames coupled to morphological examinations of the gills, confirmation at the protein level, functional characterization, etc.) in order to refine current models of calcium transport in the gill. In addition, this study also provides a valuable diversity of novel transcripts and different transcript types which will be of great interest for further exploration.

In summary, our results indicate that when T. nigroviridis are transferred to a low calcium environment, their gills respond rapidly by up regulating genes related to Ca2+ signaling/homeostasis but also actin cytoskeleton reorganization, cytokinesis and muscle contraction, which may be involved in previously described morphological adjustments in the cells where transepithelial Ca2+ transport occurs, the MR cells, in response to altered Ca2+ levels. Genes related to energy production and energy homeostasis are also up-regulated, probably reflecting the increased energetic needs of the acclimation response. The responses to high calcium availability appear to affect the same biological processes as the low calcium challenge, although with opposite effects.


Experimental procedure for environmental calcium exposure

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 Ca2+ concentrations: 2.9 mM (control group), 0.01 mM (as in ion-poor fresh water, LowCa group) or 10 mM Ca2+ (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 MgCl2; 7.1 mM MgSO4 and 0.01, 2.9 or 10 mM CaCl2, for the Control, LowCa or HighCa ABW, respectively. LowCa ABW was supplemented with 5.7 mM choline chloride (C5H14NOCl) 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 (NH3/NH4+) and nitrites (NO2-) 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 [Ca2+] 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).

RNA extraction, construction of SuperSAGE libraries and sequencing

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 data analysis

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 Gintrinsic 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 [Ca2+] 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 [Ca2+]-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.

Tag-to-gene annotation

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, 4244]. 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).

Clustering analysis

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 (GO) annotation and enrichment analysis

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).

Quantitative polymerase chain reaction

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 R2 > 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%.



18S ribosomal RNA


artificial brackish water


Ca2+-transporting ATPase sarcoplasmic reticulum type, fast twitch skeletal muscle isoform or sarcoplasmic/endoplasmatic reticulum Ca2+ ATPase1 (alternative name)

C2 h/C12 h:

exposure to control (2.9 mM Ca2+) water for 2 or 12 h, respectively


calcium ion


calcium-binding protein




calcium-sensing receptor




expressed sequenced tag


fresh water/seawater


glyceraldehyde-3-phosphate dehydrogenase


gene ontology


hemoglobin β

HighCa12 h:

exposure to water with high (10 mM) Ca2+ concentration for 12 h


leukocyte elastase inhibitor


creatine kinase muscle type 1

LowCa2 h/LowCa12 h:

exposure to water with low (0.01 mM) Ca2+ concentration for 2 or 12 h, respectively

MR cells:

mitochondrion-rich cells




plasma membrane Ca2+-ATPase


parvalbumin b


quantitative real-time PCR


ribosomal protein S18


standard error of the mean


serial analysis of gene expression


single nucleotide polymorphism


sarcoplasmic/endoplasmic reticulum


transient receptor potential cation channel subfamily V member 6 or epithelial calcium channel (alternative name).

Gene/protein names and symbols are written:

respectively, in italic or in plain type, and follow the nomenclature defined by Swiss-Prot [89]

in some cases:

an alternative name/symbol is also indicated to facilitate comparison with available bibliography. For abbreviations of genes identified by SuperSAGE, see Additional file 3 or [89].


  1. Wendelaar Bonga SE, Pang PK: Control of calcium regulating hormones in the vertebrates: parathyroid hormone, calcitonin, prolactin, and stanniocalcin. Int Rev Cytol. 1991, 128: 139-213.

    Article  CAS  PubMed  Google Scholar 

  2. Hoenderop JGJ, Nilius B, Bindels RJM: Calcium absorption across epithelia. Physiol Rev. 2005, 85: 373-422. 10.1152/physrev.00003.2004.

    Article  CAS  PubMed  Google Scholar 

  3. Guerreiro PM, Fuentes J: Control of calcium balance in fish. Fish Osmoregulation. Edited by: Baldisserotto B, Mancera Romero JM, Kapoor BG. 2007, Enfield, USA: Science Publishers

    Google Scholar 

  4. Abbink W, Bevelander GS, Hang X, Lu W, Guerreiro PM, Spanings T, Canario AV, Flik G: PTHrP regulation and calcium balance in sea bream (Sparus auratus L.) under calcium constraint. J Exp Biol. 2006, 209: 3550-3557. 10.1242/jeb.02399.

    Article  CAS  PubMed  Google Scholar 

  5. Abbink W, Bevelander GS, Rotllant J, Canario AV, Flik G: Calcium handling in Sparus auratus: effects of water and dietary calcium levels on mineral composition, cortisol and PTHrP levels. J Exp Biol. 2004, 207: 4077-4084. 10.1242/jeb.01254.

    Article  CAS  PubMed  Google Scholar 

  6. Flik G, Fenwick JC, Kolar Z, Mayer-Gostan N, Wendelaabonga SE: Effects of low ambient calcium levels on wholebody Ca2+ flux rates and internal calcium pools in the freshwater cichlid teleost, Oreochromis Mossambicus. J Exp Biol. 1986, 120: 249-264.

    Google Scholar 

  7. Guerreiro PM, Fuentes J, Flik G, Rotllant J, Power DM, Canario AV: Water calcium concentration modifies whole-body calcium uptake in sea bream larvae during short-term adaptation to altered salinities. J Exp Biol. 2004, 207: 645-653. 10.1242/jeb.00765.

    Article  CAS  PubMed  Google Scholar 

  8. Evans DH, Piermarini PM, Choe KP: The multifunctional fish gill: dominant site of gas exchange, osmoregulation, acid-base regulation, and excretion of nitrogenous waste. Physiol Rev. 2005, 85: 97-177. 10.1152/physrev.00050.2003.

    Article  CAS  PubMed  Google Scholar 

  9. Flik G, Verbost PM, Wendelaar Bongar SE: Calcium transport process in fishes. Cellular and Molecular Approaches to Fish Ionic Regulation. Edited by: Wood C, Shuttleworth T. 1995, San Diego, CA: Academic

    Google Scholar 

  10. Liao BK, Deng AN, Chen SC, Chou MY, Hwang PP: Expression and water calcium dependence of calcium transporter isoforms in zebrafish gill mitochondrion-rich cells. BMC Genomics. 2007, 8: 354-10.1186/1471-2164-8-354.

    Article  PubMed Central  PubMed  Google Scholar 

  11. Shahsavarani A, Perry SF: Hormonal and environmental regulation of the epithelial calcium channel (ECaC) in the gill of rainbow trout (Oncorhynchus mykiss). Am J Physiol Regul Integr Comp Physiol. 2006, 291: R1490-1498.

    Article  CAS  PubMed  Google Scholar 

  12. Craig PM, Wood CM, McClelland GB: Gill membrane remodeling with soft-water acclimation in zebrafish (Danio rerio). Physiol Genomics. 2007, 30: 53-60. 10.1152/physiolgenomics.00195.2006.

    Article  CAS  PubMed  Google Scholar 

  13. Pan TC, Liao BK, Huang CJ, Lin LY, Hwang PP: Epithelial Ca2+ channel expression and Ca2+ uptake in developing zebrafish. Am J Physiol Regul Integr Comp Physiol. 2005, 289: R1202-1211.

    Article  CAS  PubMed  Google Scholar 

  14. Fiol DF, Kultz D: Osmotic stress sensing and signaling in fishes. Febs J. 2007, 274: 5790-5798. 10.1111/j.1742-4658.2007.06099.x.

    Article  CAS  PubMed  Google Scholar 

  15. Abbink W, Flik G: Parathyroid hormone-related protein in teleost fish. Gen Comp Endocrinol. 2007, 152: 243-251. 10.1016/j.ygcen.2006.11.010.

    Article  CAS  PubMed  Google Scholar 

  16. Loretz CA: Extracellular calcium-sensing receptors in fishes. Comp Biochem Physiol A Mol Integr Physiol. 2008, 149: 225-245. 10.1016/j.cbpa.2008.01.037.

    Article  PubMed  Google Scholar 

  17. Wagner G, Dimattia G: The stanniocalcin family of proteins. J Exp Zool A Comp Exp Biol. 2006, 305: 769-780.

    Article  PubMed  Google Scholar 

  18. Tetraodon nigroviridis. []

  19. Lin CH, Tsai RS, Lee TH: Expression and distribution of Na, K-ATPase in gill and kidney of the spotted green pufferfish, Tetraodon nigroviridis, in response to salinity challenge. Comp Biochem Physiol A Mol Integr Physiol. 2004, 138: 287-295. 10.1016/j.cbpb.2004.04.005.

    Article  CAS  PubMed  Google Scholar 

  20. Jaillon O, Aury JM, Brunet F, Petit JL, Stange-Thomann N, Mauceli E, Bouneau L, Fischer C, Ozouf-Costaz C, Bernot A, Nicaud S, Jaffe D, Fisher S, Lutfalla G, Dossat C, Segurens B, Dasilva C, Salanoubat M, Levy M, Boudet N, Castellano S, Anthouard V, Jubin C, Castelli V, Katinka M, Vacherie B, Biemont C, Skalli Z, Cattolico L, Poulain J, De Berardinis V, Cruaud C, Duprat S, Brottier P, Coutanceau JP, Gouzy J, Parra G, Lardier G, Chapple C, McKernan KJ, McEwan P, Bosak S, Kellis M, Volff JN, Guigo R, Zody MC, Mesirov J, Lindblad-Toh K, Birren B, Nusbaum C, Kahn D, Robinson-Rechavi M, Laudet V, Schachter V, Quetier F, Saurin W, Scarpelli C, Wincker P, Lander ES, Weissenbach J, Roest Crollius H: Genome duplication in the teleost fish Tetraodon nigroviridis reveals the early vertebrate proto-karyotype. Nature. 2004, 431: 946-957. 10.1038/nature03025.

    Article  PubMed  Google Scholar 

  21. Velculescu VE, Zhang L, Vogelstein B, Kinzler KW: Serial analysis of gene expression. Science. 1995, 270: 484-487. 10.1126/science.270.5235.484.

    Article  CAS  PubMed  Google Scholar 

  22. Matsumura H, Kruger DH, Kahl G, Terauchi R: SuperSAGE: a modern platform for genome-wide quantitative transcript profiling. Curr Pharm Biotechnol. 2008, 9: 368-374. 10.2174/138920108785915157.

    Article  CAS  PubMed  Google Scholar 

  23. Matsumura H, Reich S, Ito A, Saitoh H, Kamoun S, Winter P, Kahl G, Reuter M, Kruger DH, Terauchi R: Gene expression analysis of plant host-pathogen interactions by SuperSAGE. Proc Natl Acad Sci USA. 2003, 100: 15718-15723. 10.1073/pnas.2536670100.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  24. Terauchi R, Matsumura H, Kruger DH, Kahl G: SuperSAGE: the most advanced transcriptome technology for functional genomics. Handbook of Plant Functional Genomics: Concepts and Protocols. Edited by: Kahl G, Meksem K. 2008, New York: John Wiley & Sons, 37-54.

    Chapter  Google Scholar 

  25. Greenwood MP, Flik G, Wagner GF, Balment RJ: The corpuscles of Stannius, calcium-sensing receptor, and stanniocalcin: responses to calcimimetics and physiological challenges. Endocrinology. 2009, 150: 3002-3010. 10.1210/en.2008-1758.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  26. Wang SM: Understanding SAGE data. Trends Genet. 2007, 23: 42-50. 10.1016/j.tig.2006.11.001.

    Article  PubMed  Google Scholar 

  27. Wheeler DL, Church DM, Federhen S, Lash AE, Madden TL, Pontius JU, Schuler GD, Schriml LM, Sequeira E, Tatusova TA, Wagner L: Database resources of the National Center for Biotechnology. Nucleic Acids Res. 2003, 31: 28-33. 10.1093/nar/gkg033.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  28. Ashburner M, Ball CA, Blake JA, Botstein D, Butler H, Cherry JM, Davis AP, Dolinski K, Dwight SS, Eppig JT, Harris MA, Hill DP, Issel-Tarver L, Kasarskis A, Lewis S, Matese JC, Richardson JE, Ringwald M, Rubin GM, Sherlock G: Gene ontology: tool for the unification of biology. Nat Genet. 2000, 25: 25-29. 10.1038/75556.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  29. Loh YH, Christoffels A, Brenner S, Hunziker W, Venkatesh B: Extensive expansion of the claudin gene family in the teleost fish, Fugu rubripes. Genome Res. 2004, 14: 1248-1257. 10.1101/gr.2400004.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  30. Kalujnaia S, McWilliam IS, Zaguinaiko VA, Feilen AL, Nicholson J, Hazon N, Cutler CP, Balment RJ, Cossins AR, Hughes M, Cramb G: Salinity adaptation and gene profiling analysis in the European eel (Anguilla anguilla) using microarray technology. Gen Comp Endocrinol. 2007, 152: 274-280. 10.1016/j.ygcen.2006.12.025.

    Article  CAS  PubMed  Google Scholar 

  31. Evans TG, Somero GN: A microarray-based transcriptomic time-course of hyper- and hypo-osmotic stress signaling events in the euryhaline fish Gillichthys mirabilis: osmosensors to effectors. J Exp Biol. 2008, 211: 3636-3649. 10.1242/jeb.022160.

    Article  CAS  PubMed  Google Scholar 

  32. Kalujnaia S, McWilliam IS, Zaguinaiko VA, Feilen AL, Nicholson J, Hazon N, Cutler CP, Cramb G: Transcriptomic approach to the study of osmoregulation in the European eel Anguilla anguilla. Physiol Genomics. 2007, 31: 385-401. 10.1152/physiolgenomics.00059.2007.

    Article  CAS  PubMed  Google Scholar 

  33. Boutet I, Long Ky CL, Bonhomme F: A transcriptomic approach of salinity response in the euryhaline teleost, Dicentrarchus labrax. Gene. 2006, 379: 40-50. 10.1016/j.gene.2006.04.011.

    Article  CAS  PubMed  Google Scholar 

  34. Kultz D, Fiol D, Valkova N, Gomez-Jimenez S, Chan SY, Lee J: Functional genomics and proteomics of the cellular osmotic stress response in 'non-model' organisms. J Exp Biol. 2007, 210: 1593-1601. 10.1242/jeb.000141.

    Article  CAS  PubMed  Google Scholar 

  35. Hwang P-P, Tung Y-C, Chang M-H: Effect of environmental calcium levels on calcium uptake in tilapia larvae Oreochromis mossambicus. Fish Physiol Biochem. 1996, 15: 363-370. 10.1007/BF01875578.

    Article  CAS  PubMed  Google Scholar 

  36. Chen YY, Lu FI, Hwang PP: Comparisons of calcium regulation in fish larvae. J Exp Zool A Comp Exp Biol. 2003, 295: 127-135. 10.1002/jez.a.10195.

    Article  PubMed  Google Scholar 

  37. Velculescu VE, Madden SL, Zhang L, Lash AE, Yu J, Rago C, Lal A, Wang CJ, Beaudry GA, Ciriello KM, Cook BP, Dufault MR, Ferguson AT, Gao Y, He TC, Hermeking H, Hiraldo SK, Hwang PM, Lopez MA, Luderer HF, Mathews B, Petroziello JM, Polyak K, Zawel L, Kinzler KW: Analysis of human transcriptomes. Nat Genet. 1999, 23: 387-388. 10.1038/70487.

    Article  CAS  PubMed  Google Scholar 

  38. Matsumura H, Ito A, Saitoh H, Winter P, Kahl G, Reuter M, Kruger DH, Terauchi R: SuperSAGE, Technoreview. Cell Microbiol. 2005, 7: 11-18. 10.1111/j.1462-5822.2004.00478.x.

    Article  CAS  PubMed  Google Scholar 

  39. Molina C, Rotter B, Horres R, Udupa SM, Besser B, Bellarmino L, Baum M, Matsumura H, Terauchi R, Kahl G, Winter P: SuperSAGE: the drought stress-responsive transcriptome of chickpea roots. BMC Genomics. 2008, 9: 553-10.1186/1471-2164-9-553.

    Article  PubMed Central  PubMed  Google Scholar 

  40. Zhu J, He F, Wang J, Yu J: Modeling transcriptome based on transcript-sampling data. PLoS One. 2008, 3: e1659-10.1371/journal.pone.0001659.

    Article  PubMed Central  PubMed  Google Scholar 

  41. Schaaf GJ, Ruijter JM, van Ruissen F, Zwijnenburg DA, Waaijer R, Valentijn LJ, Benit-Deekman J, van Kampen AH, Baas F, Kool M: Full transcriptome analysis of rhabdomyosarcoma, normal, and fetal skeletal muscle: statistical comparison of multiple SAGE libraries. FASEB J. 2005, 19: 404-406.

    CAS  PubMed  Google Scholar 

  42. Keime C, Damiola F, Mouchiroud D, Duret L, Gandrillon O: Identitag, a relational database for SAGE tag identification and interspecies comparison of SAGE libraries. BMC Bioinformatics. 2004, 5: 143-10.1186/1471-2105-5-143.

    Article  PubMed Central  PubMed  Google Scholar 

  43. Liang P: SAGE Genie: a suite with panoramic view of gene expression. Proc Natl Acad Sci USA. 2002, 99: 11547-11548. 10.1073/pnas.192436299.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  44. Lash AE, Tolstoshev CM, Wagner L, Schuler GD, Strausberg RL, Riggins GJ, Altschul SF: SAGEmap: a public gene expression resource. Genome Res. 2000, 10: 1051-1060. 10.1101/gr.10.7.1051.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  45. Boon W-M, Beissbarth T, Hyde L, Smyth G, Gunnersen J, Denton DA, Scott H, Tan S-S: A comparative analysis of transcribed genes in the mouse hypothalamus and neocortex reveals chromosomal clustering. Proc Natl Acad Sci USA. 2004, 101: 14972-14977. 10.1073/pnas.0406296101.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  46. Calsa T, Figueira A: Serial analysis of gene expression in sugarcane (Saccharum spp.) leaves revealed alternative C4 metabolism and putative antisense transcripts. Plant Mol Biol. 2007, 63: 745-762. 10.1007/s11103-006-9121-z.

    Article  CAS  PubMed  Google Scholar 

  47. Hoffman BG, Zavaglia B, Witzsche J, Ruiz de Algara T, Beach M, Hoodless PA, Jones SJ, Marra MA, Helgason CD: Identification of transcripts with enriched expression in the developing and adult pancreas. Genome Biol. 2008, 9: R99-10.1186/gb-2008-9-6-r99.

    Article  PubMed Central  PubMed  Google Scholar 

  48. Perry SF, Wood CM: Kinetics of branchial calcium uptake in the rainbow trout: effects of acclimation to various external calcium levels. J Exp Biol. 1985, 116: 411-433.

    Google Scholar 

  49. Katoh F, Kaneko T: Effects of environmental Ca2+ levels on branchial chloride cell morphology in freshwater-adapted killifish Fundulus heteroclitus. Fisheries Sci. 2002, 68: 347-355. 10.1046/j.1444-2906.2002.00432.x.

    Article  CAS  Google Scholar 

  50. Laurent P, Hobe H, Dunel-Erb S: The role of environmental sodium chloride relative to calcium in gill morphology of freshwater salmonid fish. Cell Tissue Res. 1985, 240: 675-692. 10.1007/BF00216356.

    Article  CAS  Google Scholar 

  51. Chang IC, Lee TH, Yang CH, Wei YY, Chou FI, Hwang PP: Morphology and function of gill mitochondria-rich cells in fish acclimated to different environments. Physiol Biochem Zool. 2001, 74: 111-119. 10.1086/319304.

    Article  CAS  PubMed  Google Scholar 

  52. Uchida K, Hasegawa S, Kaneko T: Effects of a low-Ca2+ environment on branchial chloride cell morphology in chum salmon fry and immunolocalization of Ca2+-ATPase in chloride cells. Can J Zool. 2002, 80: 1100-1108. 10.1139/z02-090.

    Article  CAS  Google Scholar 

  53. Daborn K, Cozzi RR, Marshall WS: Dynamics of pavement cell-chloride cell interactions during abrupt salinity change in Fundulus heteroclitus. J Exp Biol. 2001, 204: 1889-1899.

    CAS  PubMed  Google Scholar 

  54. Coemans B, Matsumura H, Terauchi R, Remy S, Swennen R, Sagi L: SuperSAGE combined with PCR walking allows global gene expression profiling of banana (Musa acuminata), a non-model organism. Theor Appl Genet. 2005, 111: 1118-1126. 10.1007/s00122-005-0039-7.

    Article  CAS  PubMed  Google Scholar 

  55. Tang Z, Li Y, Wan P, Li X, Zhao S, Liu B, Fan B, Zhu M, Yu M, Li K: LongSAGE analysis of skeletal muscle at three prenatal stages in Tongcheng and Landrace pigs. Genome Biol. 2007, 8: R115-10.1186/gb-2007-8-6-r115.

    Article  PubMed Central  PubMed  Google Scholar 

  56. Ling K-H, Hewitt C, Beissbarth T, Hyde L, Banerjee K, Cheah P-S, Cannon P, Hahn C, Thomas P, Smyth G, Tan S-S, Thomas T, Scott H: Molecular networks involved in mouse cerebral corticogenesis and spatio-temporal regulation of Sox4 and Sox11 novel antisense transcripts revealed by transcriptome profiling. Genome Biol. 2009, 10: R104-10.1186/gb-2009-10-10-r104.

    Article  PubMed Central  PubMed  Google Scholar 

  57. Tsai JC, Hwang PP: Effects of wheat germ agglutinin and colchicine on microtubules of the mitochondria-rich cells and Ca2+ uptake in tilapia (Oreochromis mossambicus) larvae. J Exp Biol. 1998, 201: 2263-2271.

    CAS  PubMed  Google Scholar 

  58. Nilsson S: Filament position in fish gills is influenced by a smooth muscle innervated by adrenergic nerves. J Exp Biol. 1985, 118: 433-437.

    Google Scholar 

  59. Arif SH: A Ca2+-binding protein with numerous roles and uses: parvalbumin in molecular biology and physiology. Bioessays. 2009, 31: 410-421. 10.1002/bies.200800170.

    Article  CAS  PubMed  Google Scholar 

  60. de Groof AJ, Fransen JA, Errington RJ, Willems PH, Wieringa B, Koopman WJ: The creatine kinase system is essential for optimal refill of the sarcoplasmic reticulum Ca2+ store in skeletal muscle. J Biol Chem. 2002, 277: 5275-5284. 10.1074/jbc.M108157200.

    Article  CAS  PubMed  Google Scholar 

  61. Puzianowska-Kuznicka M, Kuznicki J: The ER and ageing II: calcium homeostasis. Ageing Res Rev. 2009, 8: 160-172. 10.1016/j.arr.2009.05.002.

    Article  CAS  PubMed  Google Scholar 

  62. Perry SF: The chloride cell: structure and function in the gills of freshwater fishes. Annu Rev Physiol. 1997, 59: 325-347. 10.1146/annurev.physiol.59.1.325.

    Article  CAS  PubMed  Google Scholar 

  63. Porta AR, Bettini E, Buiakova OI, Baker H, Danho W, Margolis FL: Molecular cloning of ictacalcin: A novel calcium-binding protein from the channel catfish, Ictalurus punctatus. Mol Brain Res. 1996, 41: 81-89. 10.1016/0169-328X(96)00069-1.

    Article  CAS  PubMed  Google Scholar 

  64. Loretz CA, Pollina C, Hyodo S, Takei Y: Extracellular calcium-sensing receptor distribution in osmoregulatory and endocrine tissues of the tilapia. Gen Comp Endocrinol. 2009, 161: 216-228. 10.1016/j.ygcen.2008.12.020.

    Article  CAS  PubMed  Google Scholar 

  65. Haeseleer F, Imanishi Y, Sokal I, Filipek S, Palczewski K: Calcium-binding proteins: intracellular sensors from the calmodulin superfamily. Biochem Biophys Res Commun. 2002, 290: 615-623. 10.1006/bbrc.2001.6228.

    Article  CAS  PubMed  Google Scholar 

  66. Santamaria-Kisiel L, Rintala-Dempsey AC, Shaw GS: Calcium-dependent and -independent interactions of the S100 protein family. Biochem J. 2006, 396: 201-214. 10.1042/BJ20060195.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  67. Penniston JT, Enyedi A: Modulation of the plasma membrane Ca2+ pump. J Membr Biol. 1998, 165: 101-109. 10.1007/s002329900424.

    Article  CAS  PubMed  Google Scholar 

  68. Vanagas L, Rossi RC, Caride AJ, Filoteo AG, Strehler EE, Rossi JP: Plasma membrane calcium pump activity is affected by the membrane protein concentration: evidence for the involvement of the actin cytoskeleton. Biochim Biophys Acta. 2007, 1768: 1641-1649. 10.1016/j.bbamem.2007.03.012.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  69. Bagherie-Lachidan M, Wright SI, Kelly SP: Claudin-3 tight junction proteins in Tetraodon nigroviridis: cloning, tissue-specific expression, and a role in hydromineral balance. Am J Physiol Regul Integr Comp Physiol. 2008, 294: R1638-1647.

    Article  CAS  PubMed  Google Scholar 

  70. Bagherie-Lachidan M, Wright SI, Kelly SP: Claudin-8 and -27 tight junction proteins in puffer fish Tetraodon nigroviridis acclimated to freshwater and seawater. J Comp Physiol B. 2009, 179: 419-431. 10.1007/s00360-008-0326-0.

    Article  CAS  PubMed  Google Scholar 

  71. Tipsmark CK, Baltzegar DA, Ozden O, Grubb BJ, Borski RJ: Salinity regulates claudin mRNA and protein expression in the teleost gill. Am J Physiol Regul Integr Comp Physiol. 2008, 294: R1004-1014.

    Article  CAS  PubMed  Google Scholar 

  72. Tipsmark CK, Kiilerich P, Nilsen TO, Ebbesson LO, Stefansson SO, Madsen SS: Branchial expression patterns of claudin isoforms in Atlantic salmon during seawater acclimation and smoltification. Am J Physiol Regul Integr Comp Physiol. 2008, 294: R1563-1574.

    Article  CAS  PubMed  Google Scholar 

  73. Tipsmark CK, Luckenbach JA, Madsen SS, Kiilerich P, Borski RJ: Osmoregulation and expression of ion transport proteins and putative claudins in the gill of southern flounder (Paralichthys lethostigma). Comp Biochem Physiol A Mol Integr Physiol. 2008, 150: 265-273. 10.1016/j.cbpa.2008.03.006.

    Article  PubMed  Google Scholar 

  74. Giffard-Mena I, Boulo V, Aujoulat F, Fowden H, Castille R, Charmantier G, Cramb G: Aquaporin molecular characterization in the sea-bass (Dicentrarchus labrax): the effect of salinity on AQP1 and AQP3 expression. Comp Biochem Physiol A Mol Integr Physiol. 2007, 148: 430-444. 10.1016/j.cbpa.2007.06.002.

    Article  PubMed  Google Scholar 

  75. Cutler CP, Martinez AS, Cramb G: The role of aquaporin 3 in teleost fish. Comp Biochem Physiol A Mol Integr Physiol. 2007, 148: 82-91. 10.1016/j.cbpa.2006.09.022.

    Article  PubMed  Google Scholar 

  76. ASAB: Guidelines for the treatment of animals in behavioural research and teaching. Anim Behav. 2003, 65: 249-255. 10.1006/anbe.2003.2068.

    Article  Google Scholar 

  77. Matsumura H, Reuter M, Kruger DH, Winter P, Kahl G, Terauchi R: SuperSAGE. Methods Mol Biol. 2008, 387: 55-70.

    Article  CAS  PubMed  Google Scholar 

  78. The Gene Expression Omnibus (GEO) repository. []

  79. Schaaf GJ, van Ruissen F, van Kampen A, Kool M, Ruijter JM: Statistical comparison of two or more SAGE libraries: one tag at a time. Methods Mol Biol. 2008, 387: 151-168.

    Article  CAS  PubMed  Google Scholar 

  80. Kal AJ, van Zonneveld AJ, Benes V, van den Berg M, Koerkamp MG, Albermann K, Strack N, Ruijter JM, Richter A, Dujon B, Ansorge W, Tabak HF: Dynamics of gene expression revealed by comparison of serial analysis of gene expression transcript profiles from yeast grown on two different carbon sources. Mol Biol Cell. 1999, 10: 1859-1872.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  81. NCBI FTP directory for BLAST executables. []

  82. NCBI GenBank. []

  83. Ensembl anonymous FTP site- Tetraodon nigroviridis cDNAs. []

  84. The Tetraodon nigroviridis genome (Assembly 8, Golden Path v.2, March 2007). []

  85. Uniprot FTP repository. []

  86. Kim K, Zhang S, Jiang K, Cai L, Lee IB, Feldman LJ, Huang H: Measuring similarities between gene expression profiles through new data transformations. BMC Bioinformatics. 2007, 8: 29-10.1186/1471-2105-8-29.

    Article  PubMed Central  PubMed  Google Scholar 

  87. Benjamini Y, Hochberg Y: Controlling the false discovery rate: a practical and powerful approach to multiple testing. J R Stat Soc. 1995, 57: 289-300.

    Google Scholar 

  88. Cikos S, Bukovska A, Koppel J: Relative quantification of mRNA: comparison of methods currently used for real-time PCR data analysis. BMC Mol Biol. 2007, 8: 113-10.1186/1471-2199-8-113.

    Article  PubMed Central  PubMed  Google Scholar 

  89. Uniprot. []

Download references


We thank Dr V Schein (CCMAR) for assistance in carrying out the in vivo calcium availability experiment, Dr H Kuhl (MPIMG Berlin) for performing the 454 sequencing, Dr JM Ruijter and Dr G Schaaf (Academic Medical Centre, Amsterdam) for providing statistical software (G-test and Z-test) and for help in designing the statistical tests to detect differential expression in SuperSAGE data, Dr N Kolmakov (CCMAR) and Mr Arlindo Martins for help with Perl Scripts for tag-to-gene annotation, Dr K Kim and Dr H Huang (Department of Statistics, University of California) for providing the software and Matlab codes for clustering analysis and Gap statistics and Dr R Ben-Hamadou (CCMAR) for running the Matlab codes for Gap statistics. PISP was in receipt of fellowship BPD/25247/2005 from the Portuguese National Science Foundation. This work was partly funded by Marine Genomics Europe, project no.GOCE-CT-2004-505403, through a Technology Platform Access grant to sequence the SuperSAGE libraries at the max Plank Institute of Molecular Genetics (Berlin, Germany).

Author information

Authors and Affiliations


Corresponding author

Correspondence to Adelino VM Canário.

Additional information

Authors' contributions

PISP was involved in all aspects of the project, including experimental design and set-up, SuperSAGE libraries construction and data analysis, qPCR and writing. HM assisted in SuperSAGE libraries construction and data analysis and RT on SuperSAGE data analysis. RR was responsible for the 454 sequencing. MAST performed the GO enrichment analysis. DMP participated in the discussion of results and writing of the manuscript. AVMC devised the study, provided resources, participated in the planning of all experiments, statistical analysis, discussion of results, and writing of the manuscript. All authors read and approved the final manuscript.

Electronic supplementary material


Additional file 1: Venn diagram summarizing the mapping of differentially expressed tag sequences to the three available datasets containing T. nigroviridis DNA sequences. (PDF 48 KB)


Additional file 2: Percentage of tag mapping to different T. nigroviridis DNA datasets using different levels of stringency. (PDF 20 KB)


Additional file 3: Main data file summarizing all the expression (tag counts, expression ratios and clustering) and annotation data (BlastN and BlastX results) for the 1,339 differentially expressed tags identified. Tags are sorted by LowCa cluster number, annotation and gene symbol. Significant BlastN hits (26/26 identical nucleotides) and BlastX hits (E < 10-5) shown in white, non-significant hits in grey. * indicates BlastX hits confirmed by manual annotation. (XLS 504 KB)


Additional file 4: List of enriched GO categories in clusters of the LowCa analysis. Cluster number, GO code, GO type (P = biological process, F = molecular function, C = cellular compartment), description and adjusted p value are shown for significantly enriched categories (adjusted p value < 0.01). (XLS 86 KB)


Additional file 5: List of enriched GO categories in clusters of the 12 h analysis. Cluster number, GO code, GO type (P = biological process, F = molecular function, C = cellular compartment), description and adjusted p value are shown for significantly enriched categories (adjusted p value < 0.01). (XLS 50 KB)


Additional file 6: Detailed information for the differentially expressed tags/genes analyzed by qPCR. Includes gene accession numbers, name and sequence; tag sequence, counts, classification as differentially expressed (p < 0.05) or not, clustering and Pearson correlation coefficient and p value between SuperSAGE and qPCR results. (XLS 44 KB)


Additional file 7: Parameters used in G-test statistical comparison of tag-proportions among Tetraodon gill SuperSAGE libraries(PDF 25 KB)

Additional file 8: Primer sequences and PCR product sizes of genes selected for qPCR. (PDF 22 KB)

Authors’ original submitted files for images

Rights and permissions

Open Access This article is published under license to BioMed Central Ltd. This is an Open Access article is distributed under the terms of the Creative Commons Attribution License ( ), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Reprints and permissions

About this article

Cite this article

Pinto, P.I., Matsumura, H., Thorne, M.A. et al. Gill transcriptome response to changes in environmental calcium in the green spotted puffer fish. BMC Genomics 11, 476 (2010).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: