Dynamic transcriptomic profiles of zebrafish gills in response to zinc depletion
© Zheng et al; licensee BioMed Central Ltd. 2010
Received: 10 February 2010
Accepted: 8 October 2010
Published: 8 October 2010
Skip to main content
© Zheng et al; licensee BioMed Central Ltd. 2010
Received: 10 February 2010
Accepted: 8 October 2010
Published: 8 October 2010
Zinc deficiency is detrimental to organisms, highlighting its role as an essential micronutrient contributing to numerous biological processes. To investigate the underlying molecular events invoked by zinc depletion we performed a temporal analysis of transcriptome changes observed within the zebrafish gill. This tissue represents a model system for studying ion absorption across polarised epithelial cells as it provides a major pathway for fish to acquire zinc directly from water whilst sharing a conserved zinc transporting system with mammals.
Zebrafish were treated with either zinc-depleted (water = 2.61 μg L-1; diet = 26 mg kg-1) or zinc-adequate (water = 16.3 μg L-1; diet = 233 mg kg-1) conditions for two weeks. Gill samples were collected at five time points and transcriptome changes analysed in quintuplicate using a 16K oligonucleotide array. Of the genes represented the expression of a total of 333 transcripts showed differential regulation by zinc depletion (having a fold-change greater than 1.8 and an adjusted P-value less than 0.1, controlling for a 10% False Discovery Rate). Down-regulation was dominant at most time points and distinct sets of genes were regulated at different stages. Annotation enrichment analysis revealed that 'Developmental Process' was the most significantly overrepresented Biological Process GO term (P = 0.0006), involving 26% of all regulated genes. There was also significant bias for annotations relating to development, cell cycle, cell differentiation, gene regulation, butanoate metabolism, lysine degradation, protein tyrosin phosphatases, nucleobase, nucleoside and nucleotide metabolism, and cellular metabolic processes. Within these groupings genes associated with diabetes, bone/cartilage development, and ionocyte proliferation were especially notable. Network analysis of the temporal expression profile indicated that transcription factors foxl1, wt1, nr5a1, nr6a1, and especially, hnf4a may be key coordinators of the homeostatic response to zinc depletion.
The study revealed the complex regulatory pathways that allow the organism to subtly respond to the low-zinc condition. Many of the processes affected reflected a fundamental restructuring of the gill epithelium through reactivation of developmental programs leading to stem cell differentiation. The specific regulation of genes known to be involved in development of diabetes provides new molecular links between zinc deficiency and this disease. The present study demonstrates the importance of including the time-dimension in microarray studies.
Zinc is an essential trace element for all organisms, and is involved in a variety of biological functions [1–3]. It has been recognised as a cofactor for more than 300 catalytic enzymes, and is required for structural and functional integrity of more than 2000 transcription factors involved in the expression of various genes [3, 4]. Therefore, almost every signalling and metabolic pathway is in some way dependent on at least one, and often several, zinc-requiring proteins. In humans zinc deficiency has become a world-wide problem and causes a variety of symptoms including retarded growth, diarrhoea, anorexia, impaired immunity, skin lesions, and abnormal development . Zinc has also been implicated in many diseases including immune system defects , neurodegeneration , diabetes , and cancer . Therefore, elucidation of the molecular targets of zinc becomes exceedingly important.
Intracellular accumulation of zinc is the outcome of influx and efflux processes via zinc transporter proteins. These are mainly from the Slc39 (ZIP) transporter family, which transports zinc into the cytosol, and the Slc30 (ZnT) transporter family, responsible for the flux of zinc away from the cytosol, either into organelles or out of the cell . In addition, metallothioneins (MTs) have high binding affinity for zinc and play a very important role in maintaining stable intracellular zinc availability through the binding or releasing zinc . Levels of MT proteins depend on zinc availability and are, at least partially, regulated at the mRNA level. A well-known mechanism is through the metal-responsive transcription factor 1 (Mtf1), a recognised zinc-sensory transcriptional activator. Upon binding to zinc Mtf1 is able to bind metal-response elements (MREs), and further induces transcription of key target genes such as metallothioneins (MTs) and zinc transporter-1 Slc30a1 (Znt1) [12, 13]. Mtf1 may also function as a transcriptional repressor as exemplified by its down-regulation of zip10 [14, 15]. Mtf1-regulated genes also include several genes not related to zinc in an obvious way [14, 16, 17].
Changes in mRNA expression patterns due to zinc deficiency in some mammalian cells and tissues, such as intestine, liver, hepatocytes and leukocytes, have been investigated using microarray technology [18–23]. Some genes identified in these studies are involved in growth and energy metabolism, hepatic lipid metabolism, and signal transduction pathways that control immune responses. These data have provided a starting point to investigate molecular mechanisms mediating zinc deficiency-derived metabolic disturbances. However, although distinct sets of genes appear to be regulated by zinc deficiency, the response seems to vary from tissue to tissue, between cell types, and treatment conditions. We argue that zinc regulation is not only tissue-specific but also time-dependent, such that groups of genes with distinct functionalities may be regulated at different time-points following a change in zinc availability. In addition, as zinc is mainly obtained from diet in mammals it takes time to reach a state of zinc deficiency. It can also be difficult to distinguish the effect of zinc deficiency from other disturbances associated with poor nutrition.
In addition to absorbing zinc from the diet, fish take up zinc across the gills directly from the water. It is evident from previous studies that zinc regulation in the fish gill is indeed highly dynamic and efficient [15, 24, 25] and that zinc transporters in fish play similar functions as in mammals [26–29]. With the genomic resources available for zebrafish this is an experimental organism of choice for molecular studies on zinc. The aims of this study were to further elucidate the effects of zinc deficiency and underlying regulation mechanisms, particularly the temporal components of gene expression. Our approach therefore employs a comprehensive microarray (in-house 16K oligonucleotide arrays) to explore the time course of transcription profiles in the zebrafish gill. Zebrafish were maintained under both normal and zinc depleted conditions, and transcriptomic analysis phenotypically anchored by physiological measurements of whole body nutrient composition, zinc status and uptake. Transcriptional profiles from gill samples collected at five time points were obtained and analysed in quintuplicate. Expression of a total of 333 transcripts showed differential regulation by zinc deficiency. These were enriched with genes involved in biological processes such as development, cell cycle, cell differentiation, gene regulation, butanoate metabolism, lysine degradation, protein dephosphorylation, nucleobase, nucleoside and nucleotide metabolism, and cellular metabolic processes. DNA-binding proteins was one of the largest molecular function groups among regulated genes, and bioinformatic analysis of the temporal gene expression profile indicated that transcription factors foxl1, wt1, nr5a1, nr6a1, and especially, hnf4a may be key coordinators of the homeostatic response to zinc depletion.
The fish gill is a major pathway for zinc acquisition directly from water. To test the effect of zinc depletion on the zinc up-take capacity across the gill, unidirectional zinc influx was measured at an external zinc concentration of 5 μM Zn2+ for eight fish from each group at each time point. Fish on the zinc-deficiency treatment showed a 30-40% lower zinc uptake compared to the normal controls at all time points except Day 4 (Figure 1B).
Primers used in qPCR validation and the result from qPCR and microarray measurements.
Gene ID (symbol) a
1.21 ± 0.5
0.99 ± 0.2
1.98 ± 0.2
1.52 ± 0.1
2.42 ± 1.0
0.97 ± 0.2
1.47 ± 0.4
0.84 ± 0.1
1.36 ± 0.1
1.31 ± 0.4
1.57 ± 0.8
1.46 ± 0.8
1.36 ± 0.5
1.04 ± 0.8
2.25 ± 0.6
2.06 ± 1.0
0.72 ± 0.4
2.72 ± 1.0
0.53 ± 0.2
0.41 ± 0.1
0.38 ± 0.1
0.56 ± 0.2
0.56 ± 0.3
0.96 ± 0.5
1.08 ± 0.6
1.20 ± 0.5
1.87 ± 0.6
2.28 ± 1.2
1.26 ± 0.8
0.75 ± 0.1
2.07 ± 0.3
2.99 ± 1.5
0.32 ± 0.3
0.30 ± 0.1
1.91 ± 1.5
1.87 ± 1.2
0.21 ± 0.2
0.81 ± 0.4
3.67 ± 0.6
2.72 ± 0.9
0.53 ± 0.2
0.44 ± 0.2
4.90 ± 1.7
4.91 ± 2.9
0.15 ± 0.1
0.38 ± 0.1
7.48 ± 1.6
5.65 ± 1.4
0.91 ± 0.7
0.28 ± 0.1
13.82 ± 3.3
12.52 ± 4.4
0.74 ± 0.2
0.82 ± 0.1
1.11 ± 0.0
1.30 ± 0.4
0.60 ± 0.4
0.94 ± 0.6
1.42 ± 0.4
1.11 ± 0.4
0.89 ± 0.0
0.78 ± 0.2
0.88 ± 0.2
1.22 ± 0.2
0.98 ± 0.0
0.55 ± 0.2
2.01 ± 0.6
1.57 ± 0.4
0.92 ± 0.3
0.96 ± 0.4
1.35 ± 0.6
2.71 ± 0.9
Five fish were withdrawn from either the zinc depleted or adequate conditions at each of five time points: 8 hours, 1, 4, 7, and 14 days. Gills were dissected and RNA extracted from each pair of gills. A total of 16,331 gene targets on the zebrafish arrays were evaluated for a total of 49 RNA samples (25 from the zinc depleted fish and 24 from the controls (one control sample from Day 4 was removed because of failure to pass quality control, N = 4 or 5).
A total of 12,095 probe sets showed significant hybridisation signals (73%), which covered approximately 9,610 genes with Unigene IDs. Genes that showed a fold-change greater than 1.8 and an adjusted P-value less than 0.1, controlling for a 10% false discovery rate (FDR) with the Benjamini-Hochberg method, were regarded as differentially regulated by zinc depletion compared to the control group. We found that a total of 376 reporters, 3.1% of all significant signals in gill, were significantly regulated at one or more time points. Of these, 335 transcripts corresponding to 330 Unigene IDs were analysed further (Additional file 1, Table S1). Treating each time point as a separate event 36.6% of all significant changes in gene expression represented up-regulation and 63.4% down-regulation.
Annotation enrichment among regulated genes at all time-points using zebrafish Unigene gene identifiers.
Protein-tyrosine phosphatase, Tyr-specific/dual-specificity type
cell migration involved in gastrulation
anatomical structure development
multicellular organismal development
transferase activity, transferring one-carbon groups
phosphoprotein phosphatase activity
protein amino acid dephosphorylation
To investigate functional systems and pathways regulated by zinc deficiency in more details we took advantage of the similarities between human and zebrafish genes and used the zebrafish Unigene clusters associated with our reporters to retrieve the corresponding human orthologs. Of the 330 regulated genes with unique Unigene IDs 245 human orthologs were matched (74%, Additional file 1, Table S1) and submitted to DAVID. A total of 227 entries were returned with annotations. Genes regulated at each time point were also classified and analysed in the same manner.
'Developmental Process' was the most significantly overrepresented Biological Process GO term (P = 0.0006), involving 27% of all regulated genes (with annotation in DAVID). The 'gene expression' cluster grouped genes annotated for 'nucleus' (61 genes) and 'regulation of biological processes' (74 genes). Although not shown in Figure 3, these included 29 genes involved in 'transcription regulation', a Swiss-Prot keyword that was itself overrepresented (P = 0.048). The 'butanoate metabolism and lysine degradation' cluster comprised five genes which mapped onto the KEGG pathway for butanoate metabolism, four of which also appear in the KEGG pathway for lysine degradation. Four annotations from separate databases (InterPro domain, SwissProt Protein Information Resource, UniProt sequence feature, UniProt keyword, Gene Ontology Consortium) relating to 'protein-tyrosine phosphatases' appeared at a disproportionately higher frequency among regulated genes than expected to occur by chance (P = 0.0038 - 0.025). Four hierarchically linked GO terms relating to 'nucleobase, nucleoside and nucleotide metabolic process' were overrepresented among the regulated genes. A significantly disproportionate 75% of the regulated genes were labelled with the GO term 'cellular process' and related children terms. Among the listed top 30 most significant categories there were also four annotation terms that did not appear to be directly related to each other or any of the aforementioned clusters. These included terms such as 'regulation of neurotransmitter levels' (5 genes, P = 0.016) and 'cellular lipid metabolic process' (15 genes, P = 0.021).
Transcription factors regulated (>1.8-fold and FDR = 0.1) by zinc depletion for 8 hours to 14 days.
Forkhead box L1 (foxl1)
Iroquois homeobox protein 1, b (irx1b)
Hepatocyte nuclear factor 4, alpha (hnf4a)
Neurogenic differentiation (neurod)
Friend leukemia integration 1a (fli1a)
Forkhead box L1 (foxl1)
Hypothetical protein LOC100005466 (LOC100005466)
Distal-less homeobox gene 4a (dlx4a)
Core promoter element binding protein (copeb)
Homeo box A3a (hoxa3a)
Hypothetical protein LOC792149 (wu:fc23f06)
Nuclear receptor subfamily 5, group A, member 1a (nr5a1a)
Enhancer of zeste homolog 2 (Drosophila) (ezh2)
RING1 and YY1 binding protein (rybp)
Interleukin enhancer binding factor 2 (ilf2)
General transcription factor IIA, 1 (gtf2a1)
Similar to hypoxia-inducible factor 2 alpha (LOC566886)
NK2 homeobox 1a (nkx2.1a)
Estrogen-related receptor alpha-like (esrral)
Distal-less homeobox gene 5a (dlx5a)
Wilms tumor 1a (wt1a)
Hypothetical LOC568537 (LOC568537)
Similar to PHD finger protein 12 (LOC560119)
Notch homolog 1a (notch1a)
Hypothetical LOC573287 (LOC573287)
Hypothetical LOC568537 (LOC568537)
The fish gill provides an elegant experimental system for studying mineral absorption in vivo since the concentration and speciation of metal ions on the apical side can be carefully controlled, and this organ provides the major uptake pathway for ions, including zinc . The teleost fish gill epithelium is comprised of at least five types of cells, including respiratory pavement cells, ion transporting ionocytes (also termed 'chloride cells'), mucous cells, and neural epithelial cells . These are mounted on a cartilage scaffold forming the gill filaments and secondary lamellae where gas and ion exchange with water takes place . Gill zinc uptake is transcellular and primarily takes place in the discrete ionocytes which transport zinc apically from the water and extrude it basolaterally to the opposite side of the epithelium . Furthermore, the gill is highly dynamic, responding to changes in the external environment to which it is in intimate contact, by remodelling its structure and cellular metabolism in order to maintain the homeostasis of essential micronutrients [30, 38, 39]. In the present study the genomic resources of the zebrafish have been exploited to investigate the molecular sequences of events invoked under zinc depleted conditions. This aims to reveal the temporal sequence of events that follow treatment with zinc depletion and provides an opportunity to identify possible global regulators through bioinformatic reverse engineering.
We have established that zinc deficiency was efficiently induced in zebrafish after a 14-day depletion of zinc in both the water and diet. This indicates that any compensatory changes invoked were insufficient to fully maintain whole body zinc concentrations, which were significantly reduced at the end of the experimental period (Figure 1). Intriguingly the unidirectional gill zinc influx was found to be suppressed throughout the experimental period. This suppression of zinc uptake is counterintuitive since it may be expected that under conditions of reduced zinc availability proteins involved in zinc acquisition, such as zinc importers, would be up-regulated in an attempt to maintain zinc status. Indeed, we have shown previously data from the same experiment demonstrating increased expression of the zinc importers, zip3, zip10 and znt5, in response to zinc depletion . However, there was also down-regulation of the basolateral zinc exporter znt1, which would limit zinc transfer into the blood stream. Overall our data would suggest that zinc uptake was rate-limited by the basolateral transfer. However, further detailed physiological studies are required to confirm this interpretation.
The low zinc supply also had an effect on whole body copper status, which interestingly was reduced. The use of zinc to treat copper hyper-accumulation in Wilson's disease patients provides evidence that zinc inhibits copper uptake , an observation confirmed with our own zinc supplementation experiments . This may seem to contradict the observation that zinc deficiency also reduced whole body copper status. One possible explanation is that the reduction in zinc stores probably reduced expression of Mt, which binds copper as well as zinc. Mt binds about one third of the copper in a typical fish liver, but zinc is a much better inducer of mt expression .
The exploitation of microarray technology identified a total of 331 genes representing the molecular response of the gill epithelium to zinc depletion. These gene expression profiles from gill tissue provide an integrated image of zinc regulation occurring in a collection of cell types. Despite Zn2+ influx and whole body zinc content remaining low by the end of the 14-day experimental period, the time-course of overall gene expression indicated establishment of a new steady-state, only 21 genes remained regulated on Day 14 compared with 156 genes on Day 7. The genes regulated during the acclimation to zinc depletion are involved in a variety of biological processes, reflecting the multifaceted biochemical functions of zinc. The transitory nature of the acclimation process was characterised by a succession of functional gene categories that were regulated; transcripts for protein tyrosine phosphatases, anatomical structure morphogenesis, and cellular metabolic processes being preferentially expressed early in the experiment (Day 1), while enrichment of transcripts for proteins involved in nuclear functions, neuron differentiation, and cellular component assembly occurred much later (Day 7). Other functional categories such as 'purine ribonucleoside metabolism process', 'cardiomyopathy', and 'cellular lipid metabolic process' were significantly enriched over the entire time period but not at any particular time-point. Enrichment analysis using zebrafish gene identifiers was limited by the relatively sparse representation of zebrafish genes in the DAVID database. However, the preferential regulation of genes involved in development could be confirmed and a specific enrichment of genes involved in cell migration during development (i.e. slc39a6, wnt4a, and esrra) was observed. The latter group may be of considerable interest considering recent studies implicating zinc as a regulator of the ectoderm-mesenchyme transition (EMT) during zebrafish gastrulation and human cancer [42, 43].
The majority (72%) of the 29 regulated transcription factors were down-regulated under the zinc depleted condition. It is notable that among these genes there were 11 zinc finger proteins, and it may therefore be speculated that the effect of zinc deficiency is a feed-back response to a rather non-specific depletion of zinc from zinc finger domains. However, eight out of the 11 zinc finger proteins were down-regulated, which means that the percentage of down-regulated genes for zinc finger proteins (73%) was practically the same as that for the general cohort of regulated transcription factors. Nearly half (12) of the regulated TFs are involved in regulation of development, including foxl1, neurod1, fli1, nr5a1a, klf6, hmgb3b (zgc:11207 3), epas1, nkx2.1a, irx1, dlx4a, dlx5a, and notch1a. This may help explain the fact that developmental impairment is one of the key features of zinc deficiency . Furthermore, knockout of Mtf1 is lethal to mouse embryos but not at post-embryonic stages . It is possible that zinc signalling during development may normally modulate the expression of these TFs and their downstream targets. In terms of the significance to the fish in our experiment, it is likely that we are observing the signalling for stem-cells in the gill to differentiate and proliferate to structurally modify the gill in an attempt to regain homeostasis. The fish gill is an extremely dynamic structure shown to be reprogrammed and structurally remodelled when faced with environmental changes such as hypercapnia  and hypothermia . We propose that similar changes are occurring during zinc depletion.
Four key regulators of development containing the homeobox sequences were identified as regulated by zinc deficiency. Among these were dlx4a and dlx5a (closest human homologues DLX6 and DLX5, respectively), homologous to the Drosophila distal-less gene (dll). Zinc deficiency caused down-regulation of dlx4a and dlx5a (Figure 4) while zinc supplementation has been found to induce expression of dlx2a (closest human homologue, DLX2) another closely related gene in the same family . All members of the dlx family may be involved in formation of cartilage and/or bone and contribute to the differentiation of mineralising cell types, including osteoblasts and chondrocytes [46, 47]. This is of interest because zinc is abundant in bone and cartilage and recent studies have demonstrated the involvement of zinc transporters in the development of such tissues [48–50]. Specifically, Znt5 -null mice develop hunched backs, show impaired growth, and decrease in bone density due to osteopenia . Foetuses of mice lacking one of the essential zinc importers, slc39a4 (zip4), become critically zinc deficient and show severely underdeveloped and deformed craniofacial and limb features . Knock-out for another zinc transporter of the same family, slc39a13 (zip13), in mice results in a phenotype with defects in bone, teeth, and connective tissue, linked to changes in the BMP signalling pathway [48, 49]. This pathway involves both DLX5 and DLX6, and targeted inactivation of Dlx5 and Dlx6 genes in mice results in abnormalities very similar to those observed in Slc39a4 null mice . The direct effect of zinc on osteogenic differentiation is strongly supported by another recent study which showed that zinc deficiency suppresses matrix mineralisation and delays osteogenic activity in mouse osteoblastic MC3T3-E1 cells through downregulation of Runx2 (runt-related transcription factor 2) , while Dlx5 has been shown to induce expression of Runx2 and initiate osteogenic differenciation in chick calvaria osteoblasts in vitro . Together these data suggest that zinc deficiency, caused by reduced zinc uptake or disruption of zinc transporters, leads to dysregulation of mineralising cells in a process that involves suppression of DLX5/6 expression and downstream effects on osteoblast differentiation.
A major objective of this study was to identify possible key master regulators responsible for orchestrating the response to zinc depletion. In the present study Hnf4a (hepatocyte nuclear factor 4, alpha) was repeatedly implicated as being of potential importance for coordination of gene expression during the early stages of acclimation to zinc restriction. First, hnf4a was one of four transcription factors up-regulated upon zinc depletion at the first sampling point 8 hours (Table 3 and Figure 4C). Second, in a 'Direct Interaction Network' assembled by PathwayArchitect™ from all regulated genes during the experiment, hnf4a was positioned as a hub interacting with 55 nodes, including genes for the zinc-regulatory proteins Slc30a1 (Zip1), Slc30a6 (Zip6) and Mtf1. Third, analysis of promoter sequences of the regulated genes revealed an overrepresentation of genes with TFBSs associated with Hnf4a. Hnf4a does not require binding of a ligand to be activated, instead, binding to its cognate regulatory DNA elements (DR-1) is modified by serine/threonine and tyrosine phosphorylation, which can be catalysed by several protein kinases . Spatial and temporal specificity of target gene transactivation is believed to be controlled by expression of cofactors interacting with Hnf4a. Transcription of hnf4a can be repressed by Snai1a (Snail, snail homolog 1a (Drosophila)) , which has been shown in zebrafish to respond an increase in cellular labile Zn2+ concentration by nuclearisation and repression of target genes . Hence, the early increase in hnf4a expression is plausibly a result of a reduced inhibition by Snail, as it is likely that the sudden depletion of zinc from the water would reduce the concentration of labile Zn2+ in the gill cells. However, it was not investigated if during the first 8 hours of the experiment there was an increase in Hnf4a protein level large enough to impact expression of Hnf4a targets. Because of the profound influence of Zn2+ on protein phosphatase activities  it is also possible that a change in Hnf4a phosphorylation state contributed to the high representation of Hnf4a regulated genes. Hnf4a is known as a key regulator in the maintenance of hepatocyte differentiation and the control of lipid homeostasis . A recent study has shown that zinc supplementation can partially attenuate alcohol-induced lipid accumulation in mice by reversing alcohol-mediated inactivation of Hnf4a and Ppara . Zinc deprivation inactivated the DNA binding ability of Hnf4a and Ppara, thus regulating its target genes involved in lipid metabolism. These data not only provide further support for the key function of hnf4a in zinc deficiency found in the present study, but are also consistent with our findings of significant regulation of genes involved in 'cellular lipid metabolic process' (P = 0.021).
One of the principal mechanisms of non-genomic zinc signalling may be through its ability to inhibit protein phosphatases at nanomolar concentrations, leading to increased phosphorylation and activity of several MAP and tyrosine kinases . There is increasing evidence that this may be a common mechanism for zinc involvement in a variety of conditions including diabetes, neuronal damage during post-transfusional ischemia , and cancer . In the present study tyrosine-specific/dual-specificity phosphatases were overrepresented among genes regulated 24 hours after introduction of zinc depletion. Effects of zinc depletion on phosphatase activities in gill cells would be expected to occur within minutes . Hence, the changes in abundance of mRNA for phosphatases, primarily at the 24 hour time point, could represent a delayed transcriptional response to an earlier change in regulation of phosphatases by zinc at the protein level. Such a component of zinc signalling would be of significance to understanding of the overall effects of zinc on phosphatases and consequently MAP kinase and tyrosine kinase activities.
In conclusion, the experimental design to generate transcription profiles for gills at five time points over a 14-day period allowed us to study the changes in expression of individual genes during the period of acclimation to a zinc depleted environment. A striking feature of the transcription profiles was that although there was a continuum of functional and structural categories of genes regulated during the experiment, there was little continuity in regulation of individual genes except mt2. This oscillating pattern of gene regulation may reflect a homeostatic control mechanism that is gradually adjusting the gill to the new zinc depleted condition. We know from qPCR analysis that zinc importers were maximally up-regulated at day 7, coinciding with maximum down-regulation of the basolateral zinc exporter znt1 . Network analysis of transcriptomics data indicated that effects on phosphatases and a few key transcriptional factors, such as Hnf4a, coordinated early (up to Day 4) responses and may be involved in direct regulation of slc39a1 (zip1) and slc39a6 (zip6) and other zinc effector molecules. Subsequent events, which culminated on day 7, appeared to entail stem cell proliferation and differentiation in attempt to regain zinc homeostasis.
It is common to see microarray studies with only a single time point represented. While such studies often provide very useful information the present study, which contained the resolution of 49 microarray samples spread over two conditions and five time-points, clearly demonstrates that genes may be regulated in the opposite direction at different time-points and that sets of genes regulated at different time points may have similar ontologies but largely different identities.
The zebrafish husbandry and treatment were described previously . Briefly, juvenile zebrafish, Danio rerio (0.44 ± 0.06 g), were obtained from Neil Hardy Aquatica Ltd. (Surrey, UK). The fish were divided into three experimental groups with each group held in four identical tanks (40 fish per tank), and supplied with a continuous flow of aerated reconstituted reverse osmosis water at 26-28°C. The reconstituted water was composed of 0.6 mM NaCl, 41 μM Na2SO4, 13.6 μM KCl, 150 μM CaCl2, 3.4 μM NaHCO3, 78 μM MgCl2. In addition, each tank was equipped with a dosing system which added Zn, as ZnSO4.7H2O (BDH Chemicals), to provide a nominal zinc concentration in the tanks of 0.25 μM (16.3 μg L-1). The fish were fed a purified mash diet (Fish Nutrition Unit, University of Plymouth, UK) containing an analysed zinc concentration of 233 mg kg-1 (3.56 mmol kg-1) at a rate of 4% of their body mass per day. After one-week acclimation zinc deficiency was induced to a group of fish spread over four tanks by changing to reconstituted water without additional zinc supplementation and providing a low zinc (26 mg kg-1) diet. The water and diet conditions for one group of four tanks remained as before, serving as the control. A third group was cultured under zinc supplementation conditions, and the results of this experiment are presented in Zheng et al. . The zinc concentrations in the water were monitored daily using Inductively Coupled Plasma Mass Spectroscopy (PerkinElmer Elan DRC ICP-MS). The measured zinc concentrations for zinc depletion and control groups were 0.04 ± 0.05 μM (2.61 μg L-1) and 0.25 ± 0.09 μM (16.3 μg L-1), respectively. The experiment continued for 14 days and no mortalities or necropsies were observed.
At the end of two weeks nine fish from each group were killed by overdose of benzocaine (Sigma, USA) and analysed for whole body electrolyte and trace elements according to Handy et al.  with modifications (gills were dissected for gene expression analysis). Briefly, the individual fish was over dried at 100°C to a constant weight (for 48 h) and then digested in 2 ml of concentrated nitric acid at 70°C for 8 h. The digested samples were diluted to 6 ml with ultrapure water and then analysed for Ca, Cu, Fe, K, Na, and Zn using inductively coupled plasma optical emission spectrophotometry (ICP-OES, Varian 725-ES) against matrix matched standards. Nutritional parameters for moisture, protein, ash, and lipid were also determined, in triplicate where possible, according to Baker and Davies  with modifications. A micro-Kjeldahl digestion system was used with a 100 ml digestion tube according to the manufacturer's specification. Two fish were pooled for each protein determination and samples prepared using a Gerhardt KB40S digestion block and a Vapodest 40 distillation unit.
Unidirectional whole body influx of Zn2+ was measured as described previously . Under the assay conditions the whole body uptake of Zn2+ from the water represents almost exclusively influx across the gills. Each of the two flux bags was filled with 2 L of the reconstituted water with 5 μM Zn2+. The flux bags were equipped with an airline and placed in a thermostatically controlled water bath set at 28°C. For each time point a total of eight fish from each treatment group (two from each of the four tanks) were transferred into a flux bag. Half an hour later, 0.25 MBq of carrier-free65Zn2+ was added to each flux bag. Water samples were withdrawn from the flux-bags before and after the introduction of65Zn2+. After three hours incubation with65Zn2+ the fish were killed by benzocaine, rinsed in a 50 μM solution of zinc to remove surface-bound65Zn2+, blotted dry and weighed. The65Zn2+ radioactivity of each fish and water sample were measured in an LKB1282 CompuGamma counter, and the actual zinc concentration from the non-radioactive water sample measured by ICP-MS. The unidirectional influx of Zn2+ was calculated according to the formula, Jin= cpm (SA bw t)-1, where SA is the specific activity of65Zn2+ in the water, calculated as [65Zn2+] divided by the total [Zn2+] (cpm pmol-1), bw is the individual fish body weight, and t is the duration of flux (h).
Total RNA was extracted from dissected gill samples of individual fish using TRIzol Reagent (Invitrogen) according to the manufacturer's protocol. The RNA samples were further subject to DNase I digestion using a DNA-free kit (Ambion). The quality of the RNA samples was assessed using the Agilent 2100 Bioanalyser (Agilent Technologies). The reference RNA was purified from whole bodies of untreated zebrafish. The amino-allyl indirect labelling method was used to obtain Cy3 (reference) or Cy5 (samples)-labelled cDNA. The Cy3- or Cy5- cDNA sample was individually purified using QIAquick PCR purification columns (Qiagen).
The zebrafish oligonucleotides arrays were spotted on UltraGAPS™ coated slides (Corning Life Sciences, Promega), using a Qarray2 robot (Genetix Ltd) at the King's College London Genomics Centre, UK. The 16,399 oligonucleotides were designed and synthesised by Compugen and Sigma Genesys as a Zebrafish OligoLibrary ready set, which represent 15,806 LEADS™ clusters plus 171 controls. In addition, 331 customised oligonucleotides and 23 Cl Scorecard™ were added to the array set.
The arrayed slides were pre-hybridised in a buffer containing 25% formamide, 5× SSC, 0.1% SDS, and 1 mg/ml BSA at 42°C for 45 min, washed in double-distilled H2O three times and then air-dried. Equal amounts of Cy3 and Cy5-labelled cDNA were combined and hybridised in 1× hybridisation buffer, containing 25% formamide, 5× SSC, 0.5× Denhardt's, 0.1 μg/μl yeast tRNA, and 0.1% SDS, at 42°C for 16-20 hours in a sealed humid hybridisation chamber (Camlab). The slides were washed for 5 min sequentially in 2× SSC/0.1% SDS buffer (at 42°C), 0.1× SSC/0.1% SDS buffer and 0.1× SSC buffer, air-dried and scanned using a ScanArray (Perkin Elmer).
A total of 49 microarrays, with RNA from one fish per array, were analysed across two treatment groups and five time points. Spot intensities from each scanned image were quantified using Bluefuse software (BlueGnome) which is based on the Bayesian statistical method. Lowess normalisation was applied to the quantified data set. Probe intensities with low confidence (P < 0.05) were filtered out and the normalised data uploaded onto GeneSpring 7.3X software (Agilent Technologies). The dataset is available from Gene Expression Omnibus (GEO) under accession numbers GSE21894 and GSE21914. Profiles generated from five fish of each group at each time point were processed and regarded as biological replicates, with the exception of the day 4 time-point for the control fish, which had four biological replicates. Only genes with expression data from at least three replicate samples were considered expressed and used for subsequent differential expression analysis. Gene expression profiles of the zinc deficient group were compared to those from the control group at the corresponding time point and the significant difference determined by t -test with Benjamini-Hochberg multiple testing correction (false discovery rate, FDR). The genes were deemed as differentially expressed with a fold change (FC) greater than 1.8 and P-value less than 0.1, controlling for a 10% FDR between two groups at each time point.
The oligonucleotide reporters were mapped at the sequence level (using megablast allowing 2 bp mismatches) to Ensembl (Zv5), Refseq, and Genebank. Genebank matches were used to retrieve the relevant Unigene Cluster ID which provided linkage to appropriate human orthologues. Due to the greater proportion of human genes with association to Gene Ontology (GO) terms, functional analysis of genes of interest was performed with the ontology associated with the corresponding human orthologues (with sequence similarity higher than 45%), together with those associated with the Zebrafish Unigene IDs. The enrichment of differentially expressed genes in GO terms and other functional and structural descriptors was carried out using the online Functional Annotation Tool DAVID (Database for Annotation, Visualisation and Integrated Discovery) http://david.abcc.ncifcrf.gov/, with the zebrafish or human genome used as the background population. DAVID interprets the biological functions of listed gene and provides statistical methods for discovering enriched annotations within gene lists. Significance is calculated by a modified Fisher's Exact test with Benjamini correction.
The interactions between zinc and regulated genes were evaluated using the PathwayArchitect® software (Stratagene). The closest human orthologue to each of the genes designated as significantly regulated at one or more time points were imported into the software. Expression of zinc transporters in zebrafish gills from the same animal experiment have been analysed by qPCR and published . Therefore, the significantly regulated zinc transporters, slc39a (zip) 1, 3, and 10, and slc30a (znt) 1 and 5 were added to the list of objects. Because of its established role in zinc regulation metal-responsive transcription factor-1 (Mtf1) was also added to the component list. A Direct Interaction Network was generated by adding interactions (edges) between the entered components (nodes) from the PathwayArchitecht® database.
Upsteam DNA sequences of genes that were differentially regulated at each sampling time were obtained from Ensembl (Zv5). Consequently, only those oligonucleotide entries mapped to Ensembl genes could be used for subsequent analysis of promoter and binding motifs for transcription factors (TFs). The number of analysed sequences from the 8 hour, 1, 4, 7, and 14 day time points were 26, 46, 82, 83, and 6, respectively. In addition to the regulated genes, upstream sequence was also obtained for a cohort of 198 randomly selected genes not significantly regulated at any sampling point during the experiment. DNA sequence corresponding to 2000 bp upstream of the translation start site were downloaded to a local server and used for analysis. Differentially regulated TFs (See Table 3) were identified by their GO annotation. The FindPatterns routine of the Wisconsin Package was used to locate the transcription factor binding site (TFBS) motifs (patterns and regular expressions), as defined by the Transfac (2005) database .
A perl script was used to parse the results, counting the occurrence of each TFBS found among genes in each list (Both the perl script and the GCG output file are available on request). The TFBSs that corresponded to differentially regulated TFs were recorded. The TFBS to Mtf1 was also included. At each time point the percentage of differentially expressed genes with TFBSs corresponding to each regulated transcription factor was calculated. To assess the general frequency of occurrence of specific TFBSs in the genome, the same procedure was repeated for 198 randomly selected, non-regulated genes. The frequency of genes with TFBS to a given TF was then compared between the lists of regulated and non-regulated genes to statistically assess any deviation from a random distribution (P < 0.05, z-test).
Total RNA (2 μg) from nine gill samples collected from each group at each time point was reverse transcribed into cDNA and subjected to qPCR analysis using the SYBR® GreenER™ qPCR SuperMix kit (Invitrogen) according to the manufacturer's protocol except that 20 μl reaction volume was used. The qPCR validation was performed for 10 genes: mt2, slc30a1 (znt1), ppara, foxl1, hmgcs1, apob, cyp11a1, dlx5, dlx2, and rfx2 for which the primer sets are listed in Table 1. As a further validation of the microarray method the same set of genes were measured by qPCR on gill samples from zinc supplemented fish also analysed with the same microarray assay as that used in the present study . The qPCR assays were performed on an ABI prism 7000 with cycling conditions as follows: 5 minutes of denaturation at 95°C and then 40 cycles of 95°C for 30s, 60°C for 1 minute. A standard curve was generated for each gene using serial dilutions of a concentrated cDNA mixture to assess the amplification efficiency. The relative copy number was deduced from the corresponding standard curve using the C t value. To correct for variation in the input RNA concentration the relative gene copies were further normalised to the expression of the 18s rRNA gene by dividing the number of gene copies with the relative copies of 18s. The statistical significance of differences between experimental groups was determined with two-tailed unpaired Student's t -test (P < 0.05).
metal-responsive transcription factor 1
transcription factor binding sites
metal response element.
The authors wish to thank Dr Matthew Arno at the King's Genomics Centre for help with the microarrays and Ms Carmel at University College London for valuable comments on the manuscript. The present study was supported by funding from BBSRC grant S19512.
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.