Ethics statement.

This study was approved by the Institutional Ethical Committee of the KU Leuven (project number P102/2012) and experiments were executed in strict accordance with the European Council Directive (2010/63/EC). According to the EU Directive on the protection of animals used for scientific purposes (2010/63/EU) and the Commission Implementing Decision (2012/707/EU), fish are non-protected animals until the stage of free feeding; this limit was set at 120 hpf for zebrafish. The experiments in this study lasted only up to 96 hpf.

Microdissection and RNA isolation.

For each condition, embryos were reared together in a Petri dish (PS, 90 x 16.2 mm, MLS NV, Belgium) with a maximum number of 40 embryos. Using a stereomicroscope (M4A, Wild Heerbruggbino) and a pair of #55 Dumont forceps (Fine Tools, Berlin), 72 hpf-old zebrafish larvae were manually dissected in ice-cold RNAse-free water to separate the head, abdomen and tail regions. The head region contained the brain, sacculi/otolith and jaw parts, the abdomen region contained the yolk, swim bladder, liver and the entire gastrointestinal tract, and the tail region consisted of the remainder reaching to the tail tip. Four biological replicates of each knockdown condition (SCMO, D1D2MO, D3MO) and each body region were collected. A total of 20 head, 40 abdomen and 40 tail fragments were pooled per replicate in an RNAse-free Eppendorf tube and immediately snap frozen on dry ice. Samples were kept at -80°C until processing for RNA extraction.

Total RNA was extracted from the tissue pools using Qiazol Lysis Reagent (Qiagen, Venlo) according to the manufacturer’s instructions. RNA quality and quantity were determined using a Nanodrop spectrophotometer ND-1000 (Thermo Scientific, Rockland) while the RNA integrity was verified using the Qiaxcel qiacube (Qiagen, Venlo). Isolated intact and high quality RNA was further processed for microarray analysis.

cRNA labelling and hybridization.

Total RNA was linearly amplified and labelled using the Low Input Quick Amplification Labelling Kit (LIQA, Agilent Technologies, Santa Clara, CA, USA) according to the manufacturer’s instructions. In brief, 100 ng RNA was reverse transcribed into cDNA using oligo dT primers. Subsequently, cDNA was converted into cRNA and amplified and the resulting cRNA of each sample was labelled once with Cy3-CTP and once with Cy5-CTP. The RNeasy mini spin column kit (Qiagen, Venlo) was then used to purify the cRNA samples. cRNA yield, quality and dye incorporation efficiency were verified using a Nanodrop spectrophotometer.

We used Agilent’s Zebrafish Gene Expression Microarray V3 in a 4x44k format, which is a full genome microarray developed by Agilent and containing 43,803 Danio rerio probes. 825 ng Cy3 and 825 ng Cy5 labelled and purified cRNA were hybridized onto these microarrays for 17 h at 65°C in a rotating (10 rpm) hybridization oven (Agilent Technologies). Four biological replicates of each condition were used and an A-optimal design was used as the hybridization design for each body region separately [19].

After hybridization, slides were immersed in Agilent wash buffers, acetonitrile and in stabilization and drying solution (Agilent Technologies) to wash and protect against ozone-induced Cy5-degradation. Microarray slides were scanned using a Genepix 4100A confocal scanner (Axon Instruments, Union City, CA, USA) at a resolution of 5 μm and excitation wavelengths of 635 nm and 532 nm in an ozone-free environment (NoZone scanner enclosure, SciGene, Sunnyvale, CA, USA). Images were analyzed for spot identification and for quantification of the fluorescent signal intensities using the Genepix Pro software 6.1 (Axon Instruments).

Statistical processing of raw microarray data.

Statistical analyses of raw microarray data were performed using the R package Limma [20] as described by Vergauwen et al. [21] for each body region separately. Both knockdown conditions were contrasted against SCMO. Spots for which the criterion FG < BG + 2SD (FG: foreground, BG: background, SD: standard deviation of the local background of the entire array [22]) was true for all arrays in the dataset, were excluded from analysis. Background correction was carried out using a normal-exponential convolution model [23]. Within-array adjustment was done by Loess normalization [24], which performs an intensity-dependent normalization of the ratio of red and green. Linear models were fitted to intensity ratios, after which probes were ranked in order of evidence of differential transcription using an empirical Bayes method [20]. Contrasts were fitted to the linear models and cut off at FDR < 0.05 (false discovery rate, multiple testing correction using the Benjamini-Hochberg procedure) and |log2FC| > 0.585 (log2 fold change, corresponding to a fold induction of at least 1.5 or -1.5). Raw and analysed microarray data have been deposited in NCBI’s Gene Expression Omnibus (GEO, http://www.ncbi.nlm.nih.gov/geo) and are accessible through the GEO series accession number GSE61625 [reviewer only link: http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?token=ejojsqegzbmfnaf&acc=GSE61625].

Biological interpretation of microarray data.

Before further data analysis, for distinct probes on the array targeting identical transcripts with identical names, Agilent primary accession numbers and Genbank accession numbers, the mean of the log₂FC values was calculated. Hierarchical clustering analysis of the samples was performed using MultiExperiment Viewer software 4.8.1 (MeV, http://www.tm4.org/mev.html) based on Pearson correlation with average linkage with bootstrapping. The biological function of significantly differentially expressed transcripts in D1D2MO and D3MO was analyzed in Blast2GO (www.blast2go.org). The most prominently affected biological processes were identified using pie charts of the blast2GO score which takes into account the number of affected transcripts within each gene ontology (GO) class and the level of specificity of the respective GO class. This score was cut off at 5 to 10% of the total number of annotated transcripts in their respective condition. Subsequently, these GO classes were grouped manually into broader biological categories (e.g. metabolic process, development) to create a general overview of the data (See also S1 Table). Additionally, KEGG pathway analysis was used to identify strongly affected pathways. Only pathways containing at least 3 differentially expressed enzymes were considered (See also S1 Table).

Assessment of morphological effects.

An early life stage (ELS) test was performed in three experimental replicates for each condition, resulting in a total sample size of 120 embryos per knockdown condition. The test was performed in 48-well plates with one embryo per well. The first replicate contained 24 embryos per knockdown condition. The second and third replicates each contained 48 embryos per knockdown condition. For the morphological analysis we also included an uninjected control (UC) group and the assessing researcher was blinded with respect to the knockdown condition of the embryos. From 24 hpf up to 96 hpf, lethal and sublethal morphological and physiological endpoints were evaluated every 24 h using a stereomicroscope. The observed endpoints per time point are listed in Table 1. The heart rate of the embryos was recorded at 24 hpf and at 48 hpf by counting the heartbeats for 10 seconds (n = 96, from two experimental replicates). For morphological analysis and recording digital images at 96 hpf, embryos were lightly anaesthetized in 0.1 g/L (pH 7.3 ± 0.1) MS-222 (Tricaine, Sigma-Aldrich, Cas nr: 886-86-2). To determine the effects on growth, the length of larvae along the body axis from head to the tail tip was measured at 96 hpf on digital images (acquired using a Canon EOS 600D, 18 Mpix) using the Image J software (version 1.47).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Observed endpoints at specific time points during the zebrafish early life stage test.

https://doi.org/10.1371/journal.pone.0123285.t001

Assessment of larval motility and energy stores.

Based on the transcriptional patterns indicating effects on muscle function and energy metabolism, we performed additional analyses of motility and energy stores. For analysis of motility SCMO, D1D2MO and D3MO embryos were transferred to individual wells of a 48-well plate and were incubated up to 96 hpf. In each plate, embryos from each condition were included and the total sample size for each condition was 24. At 96 hpf, larval movements were recorded and monitored with a Zebrabox (ViewPoint, Lyon, France) for 75 minutes in a 100% full light period (1160.7 lux), using the video tracking mode with the detection threshold set at 126. The first 30 minutes were used as an acclimation period and were therefore not included in the analysis. The total distance (mm) travelled during the 45 remaining minutes was calculated.

For measurement of energy stores, SCMO, D1D2MO and D3MO embryos were raised in Petri dishes and incubated up to 72 hpf. Four biological replicates of each condition, each replicate containing 48 pooled 72 hpf-old whole larvae, were collected. Carbohydrates, lipids and proteins were extracted and measured as described earlier [21]. In brief, larvae homogenized in deionized water were used for total lipid analysis. Total lipid extraction was performed based on the chloroform-methanol protocol according to Bligh and Dyer [25] and absorption was measured at 375 nm, relative to a calibration curve of cholesterol (Acros Organics, Geel, Belgium) in chloroform. From a separate homogenate, proteins were precipitated using perchloric acid and the total protein content was determined using the Bradford method [26] by measuring absorption at 592 nm, relative to a calibration curve of bovine serum albumin (Sigma, Bornem, Belgium) in water. From the same larval homogenate, the total carbohydrate content was measured in the supernatant, using Anthrone reagent prepared according to Roe and Dailey [27] and absorption was measured at 630 nm, relative to a calibration curve of glycogen (Sigma, Bornem, Belgium) in perchloric acid. Total caloric content was calculated using the appropriate energy combustion values (39.5 kJ/g lipid, 24 kJ/g protein and 17.5 kJ/g carbohydrates [28]).

Statistical analysis of morphological and physiological data.

*Statistical tests were performed using IBM SPSS statistics (version 20) and Graphpad Prism (version 6.0). A p-value below 0.05 was considered significantly different. One-way analysis of variance (ANOVA) with Bonferroni corrected post hoc tests was used to determine significant effects between the MO groups in larval length and heart rate. A Kaplan-Meier survival curve analysis was conducted with log-rank Mantel-Cox test to analyze the hatching results in a time-dependent manner. A chi square test followed by pairwise chi square tests with Bonferroni correction was performed to identify significant differences in frequencies of malformations such as non-inflated swim bladder and pericardial oedema among all conditions. Kruskal-Wallis tests followed by Bonferroni corrected Mann-Whitney U tests were performed to determine significant differences in the travelled distance, the carbohydrate, protein and lipid contents, and the total caloric content. Knockdowns were initially compared to SCMO. If further differences were observed, additional comparisons were made.

D1D2MO and D3MO have different effects on carbohydrate and protein metabolism.

While the transcriptional response related to metabolism in D3MO was mostly limited to protein turnover (Fig 4), the response of D1D2MO was characterized by consistent upregulation of a remarkable number of transcripts coding for key proteins in carbohydrate metabolism, including glycogen breakdown, almost all steps of gluconeogenesis/glycolysis, the citric acid cycle and all five complexes of the electron transport chain (Fig 3). High TH levels are known to cause increased mitochondrial respiratory function, increased oxygen consumption, increased basal metabolic rate, and weight loss [36]. The consistent upregulation in D1D2MO tail could reflect a compensatory response to a decreased metabolic rate due to low T₃ levels. The liver is the central organ responsible for regulating carbohydrate metabolism and so far localized expression of D1 has been shown in liver of 48, 72, 96 and 120 hpf zebrafish embryos [17,31,37]. Although changes in energy metabolism in the tail are mainly attributed to muscle tissue, they are possibly indirectly affected by the liver.

In D3MO larvae, the transcriptional response involved in carbohydrate metabolism was limited, but carbohydrate stores were increased. The impact of TH on glucose metabolism has been studied extensively in humans and rats, where T₃ is known to stimulate gluconeogenesis and cause hyperglycemia [38–40]. Among other effects, T₃ stimulates hepatic glucose production and increases GLUT4-mediated transport of glucose into skeletal muscle and adipose cells [41–43]. This could explain the increased carbohydrate levels in D3MO zebrafish larvae. Although this is often related to increased expression of enzymes involved in gluconeogenesis [44], this was not observed in our study. This gluconeogenesis stimulating effect has also been attributed to an excessive availability of alanine as a substrate for gluconeogenesis, resulting from a thyrotoxicosis-associated catabolic state in isolated rat livers [45]. Correspondingly, we found decreased protein levels in D3MO larvae (Fig 7B). Increased proteolysis in skeletal muscle of hyperthyroid rats was found to be mediated mainly via the ubiquitin-proteasome pathway [46]. Among the transcriptional regulations involved in metabolism in D3MO, we mostly found transcripts involved in protein turnover including both protein degradation (proteasome subunits) and protein synthesis (Fig 4). Such conversion of proteins to carbohydrates could be partly responsible for the observed reduced growth of D3MO (Fig 6B), because structural proteins play an important part in growth (e.g. muscle) and because such conversions are energy expensive. Additionally, hyperthyroidism has been related to a decrease in metabolic efficiency by controlling the coupling of mitochondrial oxidative phosphorylation and the cycling of extramitochondrial substrate/futile cycles [47,48]. In contrast to free-feeding organisms, zebrafish embryos depend on a limited energy resource in the yolk up until the start of free-feeding around 5 days post fertilization. Therefore, the larvae analysed in this study could not compensate for a potentially increased energy loss. Although a potential increase of the metabolic rate and decrease of metabolic efficiency in D3MO did not result in a significant decrease of total available energy compared to SCMO, it could have contributed to the reduced growth. In contrast, in the case of D1D2MO, hypothyroidism is known to be related to reduced gluconeogenesis and resulting episodes of hypoglycemia [49]. Therefore, the observed upregulation of almost all enzymes involved in gluconeogenesis in D1D2MO larvae could have been a compensating mechanism facilitating maintenance of normal carbohydrate levels in D1D2MO larvae. D1D2MO also showed many upregulated transcripts involved in protein degradation and concurrently protein levels were lower than in SCMO. This could be related to decreased growth. Methimazole, an inhibitor of TH synthesis, has also been shown to cause growth retardation in 5-day-old zebrafish larvae [50] and growth arrest in 3-week-old zebrafish larvae [16]. These effects could be due to the absence of the TH stimulus needed to drive growth and development. We conclude that, even though both knockdowns resulted in decreased growth, the transcriptional pattern indicates that this was caused by different mechanisms.

TH is also known to stimulate both lipolysis and lipogenesis [51]. Since lipid levels were similar to control levels in both knockdown conditions (Fig 7C), and total energy stores were not significantly reduced compared to SCMO (Fig 7D), we conclude that altered triglyceride cycling was not an important contributor to the observed metabolic state.

Both knockdowns affect muscle development and motility.

TH is known to be involved in muscle development and muscle maturation, which is associated with a switch between myosin isoforms during postnatal development in mammals [52–54] and during metamorphosis in tadpoles [55]. In the tail of both deiodinase knockdown conditions, several transcripts coding for proteins involved in muscle development were upregulated (Fig 4). These genes included hsp90 and unc-45, involved in sarcomere formation during muscle development; mef2ca, a transcription factor regulating muscle gene expression; tropomyosin, myosin heavy polypeptide and myosin regulatory light polypeptide which are sarcomere components and desmin, which connects the sarcomere to the cytoskeleton.

We also observed decreased motility of both knockdowns (Fig 6F). Since total energy stores were not significantly decreased compared to SCMO, factors other than energy availability could be related to the observed decreased motility. We did not observe an increased incidence of malformations of the pectoral fins, malformations or irregularities of somites, nor of malformations or curvatures of the tail (Table 1). A first important contributor could be impaired muscle development, indicated by the differentially expressed transcripts involved in muscle development together with the decreased protein levels. A second factor could be impaired visual perception, for which we found evidence at the transcriptional level (See 4.1). A third factor could be impaired nervous system development, which has been linked to maternal thyroid hormone contribution in zebrafish [56]. We observed transcriptional alterations related to nervous system development in D3MO but not in D1D2MO larvae (Fig 5). Motility is needed for hatching and subsequent emergence to the water surface to inflate the swim bladder [57]. We also observed impaired swim bladder inflation in both knockdowns and delayed hatching in D3MO (Fig 6C and 6E), but we must note that it is difficult to distinguish between cause and effect in this context. For example, Liu and Chan [50] reported that T₄-supplemented larvae showed less compact myotome arrangement and that they showed no free-swimming behaviour. The authors suggested that this was due to the defective swim bladder. (See 4.3 for more detailed discussion of swim bladder inflation)

D3MO affects heart function.

TH balance is known to be essential for heart morphology and function, and T₃ has previously been shown to increase heart rate [58,59]. Furthermore, shortening of the action potential and the resulting tachycardia increases susceptibility to heart failure, known as thyrotoxic heart and is thought to be related to sympathetic adrenalin signalling [58,60]. We observed an increased heart rate at 24 hpf in D3MO embryos, while at 48 hpf the heart rate was comparable to that of the controls (Fig 6A). At 72 hpf, several transcripts involved in muscle development and contraction, including some specifically involved in heart development (e.g. cardiac myosin light chain-1), were differentially transcribed (mostly upregulated) in abdomen of D3MO in particular (Fig 5). Moreover, at 96 hpf we observed an increased occurrence of pericardial oedema, possibly a sign of heart failure (Fig 6D). Contrastingly, D1D2MO heart rate was comparable to controls combined with only little transcriptional response and no pericardial oedema.

Both knockdowns impair swim bladder inflation.

T₃ is known to drive growth and development, including the development of the swim bladder. During early development, zebrafish undergo an embryonic-to-larval transition phase marking an important switch from yolk sac- to exogenous feeding larvae around 120 hpf. This transition includes swim bladder inflation, structural and functional maturation of the mouth and gastrointestinal tract, and resorption of the yolk sac [50]. Dong et al. [17] and Thisse et al. [31] showed localized expression of D1 and D2 in the swim bladder tissue of 96 and 120 hpf zebrafish larvae. Liu and Chan [50] showed that TH plays a major role in the transition phase and reported that inhibition of the TH-TR (thyroid hormone receptor) axis retarded swim bladder development. Furthermore, the swim bladder in fish is postulated as a homolog of the tetrapod lung [61]. Holt et al. [62] showed that hypothyroidism in neonatal rats also retarded the growth of the lung. The observed reduced growth and swim bladder inflation (Fig 6B and 6C) in D1D2MO zebrafish larvae was thus expected based on the previous evidence of hypothyroidism in D1D2MO since it could be rescued by T₃ supplementation [14].

Our observation of similar (and even more prominent) effects on growth, development and swim bladder inflation in D3MO seems counterintuitive because these animals are expected to suffer from hyperthyroidism—confirmed by the previous finding that T₃ supplementation partially mimicked the effects observed in D3MO [10]. Contradicting findings have been reported on the effects of T₃ supplementation. While T₃ supplementation was shown to improve swim bladder inflation, promote growth, differentiation and survival in some studies [63–65], it was shown in some other studies to impair swim bladder morphogenesis [10], delay hatching and inhibit growth and development [10,66] similar to what we observed. It can even cause swim bladder deflation possibly by suppressing surfactant production [50]. The reason for these contradicting results may be that the final effect on growth and development depends on the dose and the thyroid status of the organism, with a beneficial effect of a small increase in T₃ levels but an adverse effect of T₃ levels above physiological boundaries.