Embryonic temperature affects muscle fibre recruitment in adult zebrafish: genome-wide changes in gene and microRNA expression associated with the transition from hyperplastic to hypertrophic growth phenotypes

In the myotomes of teleost fish, myotube production and the production of fast muscle fibres (hyperplasia) continues until around 40% of the maximum body length (Weatherley et al.,1988; Johnston et al.,2004). Once the final fibre number (FFN) is established, myotube formation is switched off and growth only involves nuclear accretion and an increase in fibre length and diameter (hypertrophy)(Weatherley et al., 1988; Johnston et al., 2004). Remarkably, in some species, such as the Atlantic salmon (Salmo salarL.), embryonic temperature treatment has been shown to have persistent effects on muscle fibre recruitment, affecting the FFN in adult stages (reviewed by Johnston, 2006). Little is known about changes in gene expression accompanying the transition between growth phenotypes in fish (Fernandes et al., 2005), and nothing about the pathways regulating muscle fibre number.

MicroRNAs (miRNAs) are an important class of 18–24 nucleotide non-coding RNAs that are repressive post-transcriptional regulators of gene expression (Bartel, 2004)involved in most if not all physiological processes, including stem cell differentiation, cell lineage specification, haematopoiesis, neurogenesis,myogenesis, immune responses, insulin secretion and cholesterol metabolism(reviewed by Williams, 2008). Computational studies suggest that around one-third of the protein-coding genes in the human (Homo sapiens) genome are subject to miRNA regulation (Lewis et al.,2005). miRNAs are derived from precursor transcripts containing hairpin structures. The ribonuclease III, Drosha, cleaves the primary transcript (pri-miRNA), releasing ∼60–80 nucleotide precursor miRNA(pre-miRNA) hairpins (Lee et al.,2003). The pre-miRNA is transported to the cytoplasm by Exportin-5 and the endonuclease Dicer processes the stem–loop to a ∼21 bp RNA duplex (Hutvagner and Zamore,2002). Although the two strands of the duplex are initially present in equal amounts, their accumulation is asymmetric at steady state,with the more abundant product referred to as the miRNA and the other stand as a miRNA^* species (Okamura et al., 2008). The miRNA is preferentially incorporated into the RNA-induced silencing (RISC) complex. Recent studies indicate that more than 40% of miRNA^* species are well conserved across Drosophilaspecies, can associate with Argonaute proteins and also have regulatory activity, adding to the richness of miRNA regulation(Okamura et al., 2008).

miRNAs are thought to block translation because the RISC component Ago 2 precludes the binding of the transcription factor eIF4e to the 7-methylguanosine cap of the target mRNA(Kiriakidou et al., 2007). Animal miRNAs show imperfect base pairing to sequences within the 3′untranslated region (UTR) of target mRNAs, although complementarity is higher for the so-called seed region (nucleotides 2–8 from the 5′ end)(Doench and Sharp, 2004). The 3′ UTRs of mRNAs form complex secondary and tertiary structures in vivo, influenced by the intracellular environment and unknown interactions with RNA transcripts and/or proteins, making it difficult to predict the single stranded regions that are accessible to miRNA binding(Zhao and Srivastava, 2007). In addition to their effects on inhibiting translation, miRNAs can affect the stability of mRNAs and mediate their degradation(Bagga et al., 2005). There is evidence that components of RISC, miRNAs and their targets are co-localised to cytoplasmic P-bodies which are thought to act as sites of programmed mRNA degradation (Chan and Slack,2006). miRNAs may therefore have an important role in phenotypic transitions by removing mRNAs that have become inappropriate for the current physiological and/or ontological state.

There are a number of miRNAs that are strongly expressed in muscle,including miR-1, miR-133 and miR-206, which interact with evolutionarily conserved and well characterised transcriptional networks involved in regulating myogenesis (Rao et al.,2006). Experimentally verified targets for miR-1 include histone deacetylase 4 (HDAC4), which represses the transcription factor MEF2C and inhibits muscle differentiation (Lu et al., 2000). miR-206 was also stimulated by the master transcription factor myoD (Rosenberg et al., 2006), although it was shown that miR-206 is mainly induced by myf5 (Sweetman et al.,2008). Switching C2C12 myoblast cultures from media promoting proliferation to media inducing differentiation resulted in a marked up-regulation of miR-1, miR-133 and miR-206(Chen et al., 2006). In vivo, functional overload leading to fibre hypertrophy in the mouse resulted in decreased expression of miR-1 and miR-133(McCarthy and Esser, 2007). Not all the miRNAs involved in the regulation of myogenesis are specifically expressed in muscle. For example, the broadly expressed miR-181 targets Hox-A11, a repressor of myoD expression. Depletion of miR-181 in C2C12 cultures reduced myoD expression and inhibited the differentiation of myoblasts to myotubes (Naguibneva et al.,2006).

The first aim of this study was to determine the effect of embryonic temperature treatment on the life-time recruitment of fast muscle fibres in the zebrafish, a tractable model species with a sequenced genome. The second aim was to use microarrays to identify global changes in mRNA and miRNA expression associated with the transition from hyperplastic to hypertrophic muscle growth phenotypes. In order to identify putative regulatory networks,we then followed a similar approach to that of Tian and colleagues(Tian et al., 2008) involving examining the reciprocal expression of a miRNA and computationally predicted target within a defined physiological context. However, whereas Tian and colleagues matched changes in miRNA and protein expression we identified regulated mRNAs that were predicted targets for reciprocally expressed miRNAs.

The F2 generation of a strain of zebrafish originally sourced from a Singapore fish farm were used in all experiments. Broodstock comprised 15 females and 8 males maintained in a freshwater re-circulation system at 26–27°C (12 h dark:12 h light photoperiod). Fish were fed twice daily with bloodworm to bring them into peak breeding condition. The 1–4 cell stage embryos from a minimum of three spawnings (∼600–1000 per spawning) were incubated at embryonic temperature treatments of 22, 26 or 31°C (±0.5°C). After hatching, the larvae were transferred to duplicate tanks maintained at a common temperature of 26–27°C (12 h light:12 h dark). Fish were initially fed 80–200 μm size fry food(ZM, Winchester, UK) and later proprietary flakes supplemented with bloodworm. Zebrafish were killed by an overdose of MS222 and pithing. Morphological measurements were made at each experimental temperature whereas expression studies were only carried out on fish reared at 26°C throughout the whole life-cycle.

Two stages were used for expression studies based on morphological measurements: (1) adults that were still actively recruiting myotubes(M⁺ stage) and (2) adults that had ceased myotube production(M^– stage). The fish from M⁺ and M^– stages were 10–12 mm TL and 28–31 mm TL,respectively. For each RNA extraction for microarray experiments, the pooled dorsal epaxial fast muscle from 10 individuals was used for the M⁺stage and from two individuals for the M^– stage. Total RNA was isolated using the mirVana miRNA isolation kit (Ambion, Austin, TX, USA)and quantified using a Nanodrop spectrophotometer (Thermo Fisher Scientific,Loughborough, UK).

The differential expression of miRNAs between M⁺ and M^– phenotypes was studied using a microarray with probes designed from the miRBase sequence database version 8.2 (Sanger Institute,Cambridge, UK; http://www.microRNA.sanger.ac.uk/sequences). Total RNA (10 μg) was analysed by LC Sciences (Houston, TX, USA) microRNA microarray service(http://www.lcsciences.com). The samples were enriched for small RNAs and labelled with fluorescent dyes:Cy3 for M⁺ and Cy5 for M^– phenotypes. A pair of labelled samples (six replicates) was hybridised to μParaFlo^Rmicrofluidic chips. Each chip includes multiple redundant regions for each miRNA. Multiple control probes were included in each chip. The stringency was estimated from the intensity ratio (>30) of perfect match and single-based match detection probes. The data were filtered and log₂ transformed and significant differences (probability values) between M⁺ and M^– phenotypes calculated using paired t-tests(Minitab, State College, PA, USA). The European Bioinformatics Institute accession number for this experiment is E-TABM-526.

As a first step to validating targets and to confirm the microarray data,we performed northern blot analysis using total RNA from independent samples comprising three individuals per phenotype. Approximately 30 μg of total RNA for each sample was separated in a denaturing 15% (m/v) polyacrylamide gel by electrophoresis. RNA was transferred to zeta probe membranes (BioRad, Hemel Hempstead, UK) using a semi-dry electro-blotting unit (Fisher Brand,Loughborough, UK) and membranes were UV cross-linked. The selected miRNA probe sequences (supplementary material Table S1) were cross-checked with the sequences listed in the miRNA registry(http://www.sanger.ac.uk/software/Rfam/mirna/). Oligonucleotides of the reverse complement to the mature miRNA were used as probes. Probes (10 μmol l^–1) were prepared by T4 polynucleotide kinase labelling of antisense oligonucleotides withγ ³²P dATP (1.5 MBq). Pre-hybridisations and hybridisations were carried out using Ultra-Hyb oligo hybridisation buffer (Ambion) at 37°C and the blots were washed twice with 0.2× SSC/0.1% (w/v) SDS at 37°C. The membranes were exposed to a phosphor storage screen, incubated at room temperature for 4–5 days and analysed using a Personal Molecular Imager FX (BioRad).

Dot blots of control probes corresponding to specific miRNAs were used to verify transfer and hybridisation, using primers with identical sequences to mature miRNAs. RNA oligonucleotides 19 and 24 nucleotides long were used as size markers. Equal loading of the gels was confirmed by re-probing the filters with a ³²P-labelled U6 RNA probe. Filters were stripped and re-probed up to 5 times. Loss of the probe was confirmed by phosphorimaging of the membrane before re-probing. Quantitative analysis of radiolabelled probes hybridising to blots was performed by auto-radiography using an Instant Imager(Canberra Packard, Meriden, CT, USA). Signals appearing on the northern blots were normalised to the corresponding U6 signal. Variations in the amount of total RNA present on the blot were calculated using U6 and used to adjust the final radioactive signals obtained to values per microgram of total RNA from the hybridisation. Statistical analysis was performed using StatView 4.01 software (Abacus Concepts, Berkeley, CA, USA), using ANOVA followed by Scheffe's F post-analyses of significance.

Hybridisations were performed using eight RNA extractions per phenotype and microarray experiments were performed using a two-colour-based gene expression system at an Agilent (Palo Alto, CA, USA) certified microarray service provider (University Health Network, Toronto, ON, Canada). Arrays were scanned using Gene pix 4000A/B scanners. Evaluation of data for microarray analysis was performed using Gene Spring software (Agilent Technologies, Mississauga,ON, Canada). Signal intensities reflected overall expression level and a detection confidence score. Signals were log₂ transformed and those that were at or below background level were discarded. Genes were filtered based both on the fold-change of ±2 and on confidence with a P-value ≤0.05. The differentially expressed genes were then clustered using the gene tree function with a Pearson correlation and average linkage. Fold-change in expression was calculated from the average signal intensity of each group and mRNAs with a fold-change ≥2 were selected for further consideration. The European Bioinformatics Institute accession number for this experiment is E-TABM-552.

The expression patterns of genes differentially expressed between M⁺ and M^– stages in the microarray experiment were examined by qPCR. A set of fast muscle RNA samples independent from those in the microarray was used. For M⁺ stages, dissections pooled from three zebrafish of 9–12 mm TL were used for each of five samples. For M^– stages, fast muscles of five individual zebrafish of 28–34 mm TL were used. Total RNA was obtained using a standard phenol:chloroform extraction method. RNA concentration and contamination carry-over were analysed using a NanoDropTM 1000 spectrophotometer (Thermo Scientific). All RNA had 260 nm:280 nm absorbance ratios between 1.9 and 2.3 and 260 nm:230 nm absorbance ratios of >2.2. The integrity of RNA was confirmed by analysing ∼1 μg of RNA by agarose gel electrophoresis and each sample had clear 28S and 18S ribosomal RNA bands with no visible RNA degradation. cDNA synthesis with 850 ng of total RNA was carried out using the QuantiTect reverse transcription kit (Qiagen, Hilden, Germany) following the manufacturer's instructions and including a genomic DNA removal step. To reduce carried-over `poisons' that might reduce qPCR efficiency, cDNA was diluted 100 times in nuclease-free water. The expression of 16 candidate genes was examined using the primers listed in supplementary material Table S2. Primers were designed to distinguish between all potential paralogues that could be identified in Ensembl (release 51; WTSI/EBI, Cambridge, UK) and NCBI zebrafish databases and so that at least one in each pair spanned an exon boundary. In three cases, this was not possible due to the limited regions available to distinguish certain highly homologous gene paralogues, although these primers were still positioned within different exons. A 4 μl sample of each cDNA was used as a template for qPCR, using 10 μl Brilliant SYBR Green QPCR master mix (Stratagene) and a Mx30005P qPCR thermocycler(Stratagene, La Jolla, CA, USA), in 20 μl reactions, performed in duplicate and containing 200 nmol of primer. Cycling parameters were as follows: one cycle of 15 min at 95°C, 40 cycles of 30 s at 95°C, 30 s at 60°C and 30 s at 72°C, followed by a DNA dissociation analysis. Sybergreen fluorescence was recorded during the extension phase of cycling. Each qPCR plate contained all sample cDNAs to avoid plate-to-plate heterogeneity. To estimate the amplification efficiency of each primer set, a cDNA dilution series was created from a pool of all cDNAs. Raw data were analysed using Mx30005P qPCR software (Stratagene) and the threshold fluorescence of dRn values was adjusted to be in the exponential phase of amplification. Cycle threshold values of samples and the dilution series were manually exported into REST 2008 (Pfaffl et al.,2002) (downloaded from http://www.gene-quantification.de/rest-2008.html),which was used to calculate reaction efficiencies and relative expression levels of M⁺ and M^– samples, normalised to the expression of two housekeeping genes (β-actin and rpl13) that were stably expressed across samples. To assess statistical differences in relative expression values, a non-parametric randomisation test was performed,using 5000 bootstrap replicates to resample the expression differences.

Briefly, a standard RT-PCR reaction was used with primers specific to myhz1(2) (supplementary material Table S2), to amplify a double-stranded cDNA product that was ligated into pCR4-TOPO T/A vector(Invitrogen, Paisley, UK) and transformed into competent Escherichia coli (Invitrogen). Cloned products were sequenced using T3/T7 primers to confirm the expected sequence and determine strand orientation. This demonstrated that the primers amplified the most divergent region of myhz(2) (the extreme 3′ UTR) which shares no more than 85%sequence identity with other myhz genes on the tandem (not shown). T3 and T7 RNA polymerases (Roche Diagnostics, Burgess Hill, West Sussex, UK) were used to synthesise RNA probes in sense and antisense directions with concurrent incorporation of digoxigenin (Roche) following the manufacturer's instructions. In situ hybridisation of probes to small bundles of fast muscle stripped from the epaxial myotomes of M⁺ and M^– zebrafish was performed using a modified standard protocol(Thisse and Thisse, 2008). Hybridised probes were detected with an alkaline phosphatase-conjugated anti-digoxigenin antibody (Roche) using NBT/BCIP (Roche). Cryosectioning of muscle bundles was performed using a cryostat (Leica Microsystems, CM1850,Nussloch, Germany) after first freezing tissues in isopentane cooled to–159°C with liquid N₂.

Computationally predicted targets associated with the significantly differentially expressed miRNAs were obtained from the miRBase Targets Version 5 database(http://microrna.sanger.ac.uk/targets/v5/). Down-regulated miRNAs were matched with the predicted up-regulated mRNA targets and vice versa. The miRBase Targets database uses the miRanda algorithm to identify potential binding sites for a given miRNA in genomic sequences. In the current version of the program, alignments require no more than one base in the `seed region' at the 5′ end of the miRNA to be non-complementary for the target to be discarded. Targets selected in this manner are further screened for thermodynamic stability of RNA folding and for conservation of alignment in the 3′ UTR of orthologous genes in at least two species (see http://microrna.sanger.ac.uk/targets/v5/info.htmlfor further details).

The evolutionary relationships of fast muscle myosin heavy chain (MyHC)genes orientated in tandem in several vertebrate genomes were reconstructed using phylogenetic analysis. Full-length amino acid sequences of 28 MyHC genes were obtained from release 51 Ensembl genome databases of zebrafish,stickleback (Gasterosteus aculeatus), tiger pufferfish (Takifugu rubripes), green-spotted pufferfish (Tetraodon nigroviridis) and human. Sequences were aligned using promals(Pei and Grishin, 2007)followed by manual alignment quality checking and removal of indels. The sequence alignment is available on request to I.A.J. Maximum likelihood was performed with Phyml (Guindon and Gascueal, 2003), using the LG model, with concurrent estimation of the γ-distribution of among-site rate variation, the number of invariable sites and employing 1000 bootstrap replicates. Neighbour joining(NJ) and maximum parsimony (MP) analyses were performed using Mega 4.0(Tamura et al., 2007),resampling the data with 1000 bootstrap iterations as a measure of branch confidence.

The relationships between the number of fast muscle fibres per myotomal cross-section and fish TL (mm) for different embryonic temperature treatments are shown in Fig. 1A. Two patterns of post-embryonic muscle fibre production were sequentially observed in zebrafish as evidenced by the spatial distribution of fibre diameters. Stratified hyperplasia involves the production of myotubes in discrete zones,resulting in layers of fibres of increasing diameter with increasing distance from the site(s) of fibre production (arrow in Fig. 1B). In the smallest fish,∼7–8 mm TL, stratified hyperplasia was the only mechanism of fibre expansion apparent, except in a few individuals from the 31°C ET treatment. In larger fish, mosaic hyperplasia was the predominant means of expansion of fibre number, with fibres of the smallest size class occurring on the surface of existing muscle fibres throughout all regions of the myotomal cross-section (illustrated by asterisks in Fig. 1C). The final number of muscle fibres (FFN) produced was estimated from the asymptote of a Gompertz curve fitted to values of fibre number and TL(Table 1). Values of FFN obtained from the Gompertz model were in good agreement with average values calculated from fish that had no fibres in the smallest size class (7–10μm; Table 1). ET treatment resulted in significant differences in FFN (P<0.01). From the model, FFN was 18.8% higher at 26°C than at 22°C (P<0.01)and 13.7% higher at 26°C than at 31°C (P<0.05). The fish length at which the recruitment of fast muscle fibres stopped, estimated from the Gompertz model, was 23.0 mm at 31°C increasing to 27.8 mm at 22°C and 29.8 mm at 26°C (Table 1). Thus, ET treatment also had a significant effect on the body length at which the transition between M⁺ and M^–phenotypes occurred. In contrast, the relationship between the cross-sectional area of fast myotomal muscle and TL was similar for all ET treatments (not shown).

One-hundred and sixty-eight miRNAs were expressed in the fast myotomal muscle of adult zebrafish reared at 26–27°C over the whole life-cycle out of 219 miRNAs on the microarray, although only 75 were consistently expressed in all individuals. We assessed the relative expression of miRNAs on the basis of their background-subtracted and normalised fluorescence intensity signals and classified them as high(20,000–55,000 units), moderate (10,000–19,999 units) or low(2000–9900) abundance. The most abundant miRNAs were miR-1, let-7a,let-7c, let-7f, miR-17a, miR20a,b, miR-126, miR-133c, miR-181a, miR-203b,miR-206, miR-214 and miR-738. As absolute expression is influenced by variation in probe concentration and the efficiency of printing pins, this result should be treated with caution. We identified 14 miRNAs that were up-regulated (Fig. 2) and 15 miRNAs that were down-regulated in the M^– phenotype in 6/6 individuals (Fig. 3). The expression of a selection of the differentially regulated miRNAs was successfully validated by northern blotting (supplementary material Table S1; Fig. 4). The signal intensity of northerns was quantified and normalised to U6, and found to correlate well with the microarray data (Fig. 4). Five members of the dre-let-7 family of miRNA (let-7b, e, g,h, j) were consistently up-regulated in the M^– phenotype(Fig. 2). Three members of the dre-miR-19 family (miR-19b, c, d) and two members of the dre-miR-130 family(miR-130b, c) were significantly down-regulated in the M^–phenotype (Fig. 3). The most down-regulated miRNA in the M^– phenotype by 47/34-fold(microarray/northerns) was dre-miR-9^*, the expression of which was relatively low (Figs 3 and 4).

A whole genome zebrafish array was used to identify genes differentially expressed by at least 2-fold between 8/8 M⁺ and M^–samples with a significance level of P<0.05. The identity of the genes was confirmed by blasting the probe sequence on the array against the zebrafish genome assembly version 7 (Zv7; http://ensembl.genomics.org.cn/Danio_rerio/index.html)and non-redundant nucleotide and EST databases(http://blast.ncbi.nlm.nih.gov/). When closely related putative paralogues were identified for sequences on the array, the identity shared by probe sequences and corresponding regions in their paralogues was calculated to check for potential spurious cross-hybridisation artifacts. Following this analysis there were 34 down-regulated genes (Table 2)and 30 up-regulated genes (Table 3) in the M^– phenotype(Fig. 5). The three most down-regulated genes (by 10- to 20-fold) in the M^– phenotype together with the sixth (by ∼9-fold) and eighth (by 7-fold) most down-regulated genes corresponded to MyHC genes(Table 2). These MyHC genes were part of a tandem repeat on chromosome 5(Fig. 6). Expression of all five down-regulated MyHC genes was validated by qPCR, using highly specific primers and fold down-regulation was generally even greater than that observed on the microarray (Fig. 6). The final member of the cluster not detected by the microarray analysis was shown by qPCR to be significantly up-regulated in the M^– phenotype(Fig. 6). The spatial expression of one of the MyHC genes [myhz1(2)] was investigated by in situ hybridisation. myhz1(2) was highly expressed in nascent muscle fibres, but not in larger diameter fibres(Fig. 6). The discovery of these MyHC genes using the approach adopted is highly encouraging because the M^– phenotype was defined in terms of the absence of myotubes and the smallest size class of muscle fibre.

The genomes of several other vertebrates have similar tandems of fast skeletal muscle MyHC genes [e.g. stickleback (Gasterosteus aculeatus), human, Fig. 6;medaka (Oryzias latipes) (Liang et al., 2007)]. A maximum likelihood phylogenetic analysis was performed for complete amino acid sequences of these genes for zebrafish,stickleback and human, within a framework containing in total 13 human MyHC genes (Fig. 6). To test the sensitivity of the tree topology to the reconstruction method, neighbour joining and maximum parsimony analyses were also performed on the same data,producing highly comparable well-supported trees (not shown). All the included vertebrate fast skeletal MyHC sequences found in tandems branched as a clade internal to other MyHC types, many of which are conserved in both teleosts and mammals, and including cardiac and other non-fast skeletal muscle isoforms(Fig. 6). Zebrafish and stickleback fast skeletal MyHC cluster sequences branched internally to myh13 of the human cluster. Therefore, it is possible that the tandem zebrafish/stickleback fast skeletal MyHC genes and their single orthologues in pufferfish are actually co-orthologues of myh13 and that other human MyHC genes on the tandem were derived independently. Further, zebrafish and stickleback clusters form separate clades(Fig. 6), suggesting that these tandem arrangements were independently derived after the speciation event separating these lineages.

The relatively small numbers of differentially expressed genes were classified manually according to their function based on a search of the literature primarily using the PubMed, GoogleScholar, iHop(http://www.ihop-net.org/)and Kegg(http://www.genome.jp/kegg/)databases. The most abundant categories of down-regulated genes in the M^– phenotype were contractile proteins (23.5%) and sarcomeric structural/cytoskeletal proteins (14.7%, Table 2; Fig. 5). The next most abundant category of down-regulated genes in the M^– phenotype was involved with either tyrosine metabolism or amino acid transport. Two genes encoding transcription factors were significantly down-regulated in the M^– phenotype on the array, the myogenic regulatory factor myf5 and sox11a (Table 2). Cystathionine γ-lyase (cth), which catalyses the production of gaseous H₂S from cysteine and functions as a neuromodulator and physiological vasodilator involved in the regulation of blood pressure (Yang et al.,2008), was down-regulated 3.2-fold in the M^–phenotype (Table 2). There were two genes significantly down-regulated in the M^– phenotype that were involved in tyrosine metabolism (an orthologue of 4-hydroxy-phenyl pyruvate dioxygenase and fumarylacetoacetate hydrolase, fah) and melanin biosynthesis (tyrosine-related protein 1b, tyrp1b), perhaps reflecting some inadvertent contamination with pigment cells in the M⁺ phenotype samples.

The significantly up-regulated genes in the M^– phenotype showed a very different profile with only contractile proteins and energy metabolism genes represented in the functional categories observed for the down-regulated genes (Fig. 5). pvalb4, which is a member of the parvalbumin gene family that code for sarcoplasmic Ca²⁺ binding proteins involved in muscle relaxation (Jiang et al.,1996), was up-regulated 8.3-fold on the microarrays(Table 3) and 16.9-fold by qPCR(Table 4). Immune-related genes related to cell–cell interactions and cytokine pathways comprised around 15% of the up-regulated genes and included β-2-microglobulin precursor(b2m, 7.3-fold on array and 5.1-fold by qPCR), CD9 antigen-like(cd9l, 2.6-fold on array and 2.4-fold by qPCR), invariant chain-like protein 1 (iclp1) and tyrosine kinase, non-receptor 2 (tnk2, Table 3). Enolase 3 (eno3), a myoblast-specific enhancer was up-regulated 2.2-fold on the array and hypoxia-inducible factor 1, α-subunit inhibitor (hif1an) was up-regulated 2.8-fold on the array and 4.2-fold by qPCR (Tables 3 and 4). This latter protein functions as part of an oxygen-sensing system in muscle and under normoxic conditions hydroxylation of the C-terminal transactivation domain of HIF-1α by hif1an represses its transcription(Semenza, 1999). Two components of G-protein signalling were significantly up-regulated in the M^– phenotype, a regulator of G protein signalling 5(zgc:64006) and ADP-ribosylation factor-like 6 interacting protein 1(Tables 3 and 4).

Seven of the down-regulated genes (20.5%) in the M^–phenotype were predicted targets for significantly up-regulated miRNAs(Table 2). Two of these genes were predicted targets for miRNAs (three miRNAs for fah and four miRNAs for aspn; Table 2). Dre-miR365 was predicted to target three mRNAs (fzd8awhich is a WNT inhibitor, aspn and fkbp1b; Table 2). Seven of the up-regulated genes in the M^– phenotype were predicted targets for significantly down-regulated miRNAs(Table 3). Those mRNAs predicted to be targets for more than two miRNAs were pvalb4 with four, glyceraldehyde 3-phosphate dehydrogenase (gapdh) with four, and regulator of G-protein signalling (zgc:64006) with three(Table 3). Dre-miR-181c was predicted to bind to the 3′ UTR of pvalb4, gapdh, cd9l,zgc:64006 and slc25a4 (Table 3).

Traditionally the expansion of muscle fibre number in teleosts has been described in morphological terms as the formation of myotubes at discrete germinal zones (stratified hyperplasia) or myotube production throughout the myotome (mosaic hyperplasia) (Rowlerson and Veggetti, 2001). Mosaic hyperplasia is absent in some fish species that only attain a small ultimate body length, e.g. the guppy(Poecilia reticulata) (Veggetti et al., 1993). It has been reported in the literature that mosaic hyperplasia is also absent in zebrafish restricting its use as a general model for teleost post-embryonic myogenesis (see Biga and Goetz, 2006), although this study examined only relatively small fish. Recently mosaic hyperplasia was observed in zebrafish larva at 6 mm standard length (∼7.5 mm TL)(Patterson et al., 2008). The present study considerably extends these observations and demonstrates that mosaic hyperplasia is the main mechanism of muscle fibre expansion in zebrafish, as is the case for larger fish of commercial importance such as Atlantic salmon (Salmo salar L.) (see Johnston, 2006). This is an important finding because, compared with most farmed species, zebrafish has considerable advantages for investigation of the genetic mechanisms controlling growth, including the impact of selection and domestication. For example, in addition to a sequenced genome(http://www.ensembl.org/Danio_rerio/index.html)the small size and short generation time (∼4 months) of zebrafish are ideal for maintaining replicated selected and non-selected lines at relatively low cost.

Embryonic temperature has been shown to alter the number and diameter of fast and slow myotomal muscle fibres in a phylogentically diverse range of teleost species (reviewed by Johnston,2006). However, only a few studies have determined the long-term consequences of embryonic temperature for muscle growth in adult stages(Johnston et al., 2003; Macqueen et al., 2008; López-Albors et al.,2008). Macqueen and colleagues(Macqueen et al., 2008)incubated Atlantic salmon embryos at 2, 5, 8 or 10°C until the completion of eye pigmentation and then transferred them to common rearing conditions. Fish at lower temperatures remained smaller until smoltification 18 months later, but showed substantial compensatory catch-up growth in seawater over the next 18 months. The final number of fast muscle fibres was highest for the 5°C treatment and reduced at higher and lower treatments(Macqueen et al., 2008). ET treatment was also shown to alter the number of myonuclei per centimetre of fibre length in isolated single muscle fibres in this species(Johnston et al., 2003). In the present study, we found that zebrafish showed an optimal embryonic temperature for FFN of 26°C, which resulted in 18.8% more fast fibres than at 22°C and 13.7% more fibres than at 31°C(Table 1). Therefore zebrafish provides a good model for developmental plasticity to temperature in commercial species such as Atlantic salmon(Fig. 1), but has the advantage that the outcome of embryonic treatment on FFN can be established in less than 3 months. The present study showed that embryonic temperature affects both the intensity of myotube production (Fig. 1) and the body length at which the transition between M⁺ and M^– phenotypes is completed and the FFN established (Table 1),extending previous studies. These observations require direct temperature effects on embryonic tissues such as the myogenic stem cell containing external cell layer. Cell lineage and vital dye tracking studies in zebrafish have shown that during mid-segmentation the somites undergo a 90 deg. rotation from their starting positions (Hollway et al., 2007). The cells in the posterior somite domain differentiate into the primary embryonic fast muscle fibres whereas those in the anterior compartment form the external cell layer on the outside of the embryonic slow muscle layer (Hollway et al.,2007; Stellabotte et al.,2007). Cells derived from the Pax3/7-expressing external cells migrate through the somite to form additional fast muscle fibres in the late embryo and larval stages (Hollway et al.,2007; Stellabotte et al.,2007). As the external cell layer persists in later stages it is a strong candidate for providing some or all of the myogenic progenitor cells required for juvenile and adult growth(Hollway et al., 2007; Stellabotte et al., 2007). We next used microarrays to obtain genome-wide information on changes in miRNA and mRNA expression between the M⁺ and M^–phenotypes.

We have chosen to investigate changes in gene and regulatory RNA expression between two complex growth phenotypes delineated by the active production of myotubes in fast myotomal muscle. Growth involves a population(s) of myogenic progenitor cells (MPCs) or myoblasts that remain capable of proliferation and are regulated by signalling pathways responsive to both nutritional status and environmental conditions. Myoblast fusion involves several processes,including the recognition and adhesion of myoblasts, the breakdown of muscle membranes and the remodelling of the actin cytoskeleton(Richardson et al., 2008). The primary event in myotube formation is myoblast–myoblast fusion giving rise to a syncytial structure with several nuclei. The secondary events of myotube elongation involve the accretion of a large number of additional nuclei and appear to involve distinct myoblast–myotube fusion events and separate regulatory pathways (Horsley and Pavlath, 2004). Genetic analysis involving Drosophila has identified a set of genes that have conserved functions in myotube formation across the metazoans (Richardson et al.,2008). Myoblast fusion and muscle formation were disrupted in zebrafish embryos lacking a functional Rac gene, which encodes a small GTPase that regulates the actin cytoskeleton (Pajinici et al., 2008). In mammals, two conserved orthologues of Drosophila Bag2 and Dock180 respectively activated the GTPases AFR6 (ADP ribosylation factor 6) and Rac and were required for myoblast fusion and myotube differentiation (Pajcini et al., 2008). Knockdowns of Dock1 and Dock5orthologues of the Drosophila gene myoblast city (Mbc) also result in the failure of myoblasts to fuse(Moore et al., 2007). Several genes and pathways have been discovered that are associated with the secondary events of myotube formation including the cytokine coding interleukin 4, IL-4(Horseley et al., 2003) and myoferlin(Doherty et al., 2005). For example, the transcription factor NFATC2 regulates secretion of IL-4, which is essential for nuclear accretion during myotube elongation. Myoblasts derived from NFATC2^–/– mice still form thin syncitial structures with a few nuclei associated with primary myotube formation, but fail to recruit additional nuclei and increase in diameter in the same way as cultures from wild-type animals (Horsley et al., 2003). In mammals another cytokine, interleukin-6 (IL-6),which is locally and transiently produced by growing myofibres and associated satellite cells is involved in muscle fibre hypertrophy(Serrano et al., 2008). There is evidence that IL-6 deficiency impairs myoblast proliferation and myonuclear accretion in growing muscle by impairing STAT3 activation and expression of its target gene cyclin D1(Serrano et al., 2008). It is therefore to be expected that myotube formation (specific to the M⁺phenotype) would share some common mechanisms and patterns of gene expression with muscle fibre growth involving fibre hypertrophy (present in the M⁺ and M^– phenotypes) as well as having its own distinct features.

Genes associated with sarcomere structural proteins and the cytoskeleton comprised ∼15% of the significantly down-regulated genes in the M^– phenotype (Fig. 5). Thymosin β4, the seventh most down-regulated gene in the M^– phenotype, is an actin monomer-sequestering protein regulating unpolymerised actin to control the assembly of microfilaments(Dedova et al., 2006), which has been implicated in promoting cell migration, angiogenesis, cell survival and wound healing. Studies with C2C12 myoblasts have shown that promyogenic members of the Ig superfamily bind to each other in a cis fashion,forming complexes with N- and M-cadherin. These complexes containβ-catenin and are enriched at sites of cell–cell contact between myoblasts (Kang et al., 2003). In the M^– phenotype of zebrafish, ctnna2, an orthologue of human catenin, was down-regulated 2.2-fold(Table 2). Another of the down-regulated genes in the zebrafish M^– phenotype was aspn, which is one of the class I members of the small leucine-rich repeat proteoglycans (SLRPs) which also include decorin and biglycan(Henry et al., 2001). In the mouse, asporin is strongly expressed in the skeleton and more weakly in the fascia surrounding muscle fibres (Henry et al., 2001). Asporin inhibits TGF-B/Smad signalling by colocalising with TGFB-1 on the cell surface and inhibiting its binding to the TGFB type II receptor (Nakajima et al.,2007). Knockdown of asporin by small interfering RNA (siRNA)inhibits TGF-B1-induced gene expression and blocks chondrogenesis(Nakajima et al., 2007)whereas targeted mutations of class I SRLPs result in abnormal collagen fibril formation (Henry et al.,2001).

The second largest category of up-regulated genes in the M^– phenotype were immune-related genes, including several that function in cytokine pathways and are therefore candidates for involvement with myoblast fusion and/or muscle hypertrophy. Several transmembrane proteins containing immunoglobulin domains function in the recognition and adhesion of myoblasts(Richardson et al., 2008). The second most highly up-regulated gene in the M^– phenotype was b2m which functions in the folding, peptide binding and surface display of class I antigens (Yu et al.,2009). CD9 antigen-like (cd9l), which was down-regulated∼2.5 times (Tables 3 and 4), is a cell surface molecule that interacts with integrins and other membrane proteins. Most cytokine receptors are capable of recruiting and/or activating non-receptor protein kinases that induce downstream signalling pathways(Taniguchi, 1995). Tyrosine kinase, non-receptor 2 (tnk2) was up-regulated 2.2-fold in the M^– phenotype (Table 3). There is evidence that macrophages are involved in muscle regeneration and can stimulate myogenic cell growth in vitropromoting myoblast fusion into myotubes and myogenin expression leading to differentiation (Arnold et al.,2007). Invariant chain-like protein (CD74 antigen) has MHC2 interacting and thyroglobulin domains.

Our gene expression analysis has successfully identified five paralogues of fast skeletal myosin heavy chain organised in a tandem repeat on chromosome 5 of zebrafish that are very highly down-regulated in the M^–phenotype and specifically expressed in very small diameter muscle fibres(Fig. 6). Several of these genes are also highly expressed in the somites of zebrafish embryos(Xu et al., 2000). Myosin heavy chain isoforms specific to small diameter muscle fibres have previously been reported in the common carp (Cyprinus carpio L.)(Ennion et al., 1999). Interestingly, phylogenetic analysis suggests that these tandem copies and similar clusters of orthologous myosin heavy chain genes found in other vertebrates including stickleback and human(Fig. 6) as well as medaka(Liang et al., 2007) are not synapomorphies and were derived separately in each of these lineages. This suggests that some selective advantage exists, at least in some vertebrates,for having multiple tandem copies of MyHC genes. However, considering that one of the zebrafish tandem copies was not down-regulated in M^–zebrafish, it is possible that complex lineage-specific patterns of myosin heavy chain gene regulation has occurred, which might contribute to species-specific differences in fast-twitch myotube formation patterns.

Numerous aspects of muscle phenotype can be correlated with changes in body size. For example, the maximum tail-beat frequency (contraction duration per cycle) and aerobic metabolic capacity are known to decrease with increasing body length (James et al.,1998; Davies and Moyes,2007). We found a huge fold increase in pvalb4 expression in the M^– phenotype (Tables 3 and 4). Pvalb4 is a cytoplasmic Ca²⁺ binding protein involved in muscle relaxation. In the rainbow trout, the content of parvalbumin isoform 1 was shown to decrease along the trunk and was associated with a slowing of muscle relaxation rate(Coughlin et al., 2007), as occurs with increasing body size. An orthologue of the myozenin gene(myoz1; calsarcin–calcineurin binding protein) was down-regulated 2.3-fold in the M^– phenotype on the microarray(Table 2), although no significant difference in expression was observed by qPCR when independent M⁺ and M^– samples were used(Table 4). Myoz1knockout mice are deficient in calsarcin-2 and show enhanced NFAT activity and calcineurin signalling leading to a slower oxidative phenotype(Frey et al., 2008). The higher mRNA levels of myoz1 observed in the M⁺ phenotype on the array may be related to the increased aerobic character of fast muscle observed in small compared with large fish(Davies and Moyes, 2007).

Fourteen up-regulated (Fig. 2) and 15 down-regulated miRNAs(Fig. 3) were identified in the M^– phenotype providing evidence for the involvement of miRNAs in muscle growth transitions; 57% of the down-regulated mRNAs and 73% of the up-regulated mRNAs were predicted targets for one or more differentially expressed miRNAs (Figs 2 and 3; Tables 2 and 3). Bioinformatic approaches to identify mRNA targets for miRNAs have involved assessing Watson–Crick base-pairing at nucleotides 2–7 at the 5′ end of the miRNA (the so-called `seed match') (Brennecke et al.,2005). However, many computationally predicted targets have failed to be confirmed experimentally (Didiano and Hobert, 2006) and some validated miRNAs were not identified by the current algorithms (Nicolas et al.,2008). miRNA target site interactions may also involve local accessibility of the binding site. Validated miRNA target sites often have destabilising elements or high free energy in regions flanking the 5′ or 3′ ends of the target site (Xiao et al., 2009). The popular computational miRNA–mRNA prediction algorithm miRanda detects potential target sites based on the alignment score and minimum free energy (MFE) of the miRNA bound to the potential target site,and is likely to overestimate the number of true targets(Moxon et al., 2008). Thus the strongest candidate miRNAs to have a role in this growth transition are those with multiple targets. Of particular interest was dre-miR-181c which was expressed at a 2-fold lower level in the M^– phenotype and was predicted to bind to the 3′ UTR of five of the up-regulated genes pvalb4, gapdh, cd9l, zgc:64006 and slc25a4(Table 3). It is of interest that miR-181 has been shown to be up-regulated before or at the same time as muscle differentiation markers such as creatine kinase in cell culture(Naguibeva et al., 2006). In vivo miR-181 was weakly expressed in adult mouse tibial muscle, but was strongly up-regulated following repair from injury (Naguibeva et al., 2006). miR-181 was also shown to repress the translation of Hox-A11 a repressor of the differentiation program (Naguibeva et al., 2006).

MiR-1 and miR-206 are known from functional and expression studies to interact with conserved transcriptional networks regulating myogenesis(Callas et al., 2008) and were significantly differentially expressed between phenotypes, but at less than 2-fold (not shown). For example, miR-206 was 31% higher in the M⁺phenotype, which contained actively differentiating muscle fibres, whereas miR-133c expression was not significantly different between phenotypes. These two miRNAs have been shown to be induced by transferring C2C12 myoblasts to differentiation medium (Kim et al.,2006). Although miR-206 transfection advanced myosin heavy chain expression after changing to differentiation medium, miR-133 transfection did not. Inhibition of miR-206 by antisense oligonucleotide inhibited cell cycle withdrawal and differentiation and evidence was presented that mRNA for the p180 subunit of DNA polymerase was degraded by miR-206(Kim et al., 2006). miR-206 also regulates the expression of connexin43, a component of gap junctions required for the fusion of myoblasts and muscle differentiation in vitro (Anderson et al.,2006). However, in this case regulation occurs by inhibiting translation without targeting the mRNA for degradation(Anderson et al., 2006).