- Research article
- Open Access
Cross-tissue and cross-species analysis of gene expression in skeletal muscle and electric organ of African weakly-electric fish (Teleostei; Mormyridae)
BMC Genomicsvolume 16, Article number: 668 (2015)
African weakly-electric fishes of the family Mormyridae are able to produce and perceive weak electric signals (typically less than one volt in amplitude) owing to the presence of a specialized, muscle-derived electric organ (EO) in their tail region. Such electric signals, also known as Electric Organ Discharges (EODs), are used for objects/prey localization, for the identification of conspecifics, and in social and reproductive behaviour. This feature might have promoted the adaptive radiation of this family by acting as an effective pre-zygotic isolation mechanism. Despite the physiological and evolutionary importance of this trait, the investigation of the genetic basis of its function and modification has so far remained limited. In this study, we aim at: i) identifying constitutive differences in terms of gene expression between electric organ and skeletal muscle (SM) in two mormyrid species of the genus Campylomormyrus: C. compressirostris and C. tshokwe, and ii) exploring cross-specific patterns of gene expression within the two tissues among C. compressirostris, C. tshokwe, and the outgroup species Gnathonemus petersii.
Twelve paired-end (100 bp) strand-specific RNA-seq Illumina libraries were sequenced, producing circa 330 M quality-filtered short read pairs. The obtained reads were assembled de novo into four reference transcriptomes. In silico cross-tissue DE-analysis allowed us to identify 271 shared differentially expressed genes between EO and SM in C. compressirostris and C.tshokwe. Many of these genes correspond to myogenic factors, ion channels and pumps, and genes involved in several metabolic pathways. Cross-species analysis has revealed that the electric organ transcriptome is more variable in terms of gene expression levels across species than the skeletal muscle transcriptome.
The data obtained indicate that: i) the loss of contractile activity and the decoupling of the excitation-contraction processes are reflected by the down-regulation of the corresponding genes in the electric organ’s transcriptome; ii) the metabolic activity of the EO might be specialized towards the production and turn-over of membrane structures; iii) several ion channels are highly expressed in the EO in order to increase excitability; iv) several myogenic factors might be down-regulated by transcription repressors in the EO.
Bioelectrogenesis (i.e., the ability to produce strong or weak electric signals by specialized organs) has evolved several times independently in aquatic vertebrates . In fact, it can be observed in the marine electric rays (Torpediniformes) and skates (Rajiformes), in the African freshwater Mormyridae and Gymnarchidae (Osteoglossiformes; Mormyroidea), in the South American knifefishes (Gymnotiformes), in several catfish species (Siluriformes), and in few marine stargazers (Perciformes; Uranoscopidae). In all the above-mentioned groups, electric organs originate from myogenic tissue; the only exception are members of the family Apteronotidae (Gymnotiformes), where the electric organs are formed by modified spinal motor neurons . The amount of excitable cells within each electric organ determines the electric potential of an EOD, which can range from few millivolts to several hundreds of volts (e.g., in the electric eel Electrophorus electricus) .
African weakly-electric fishes of the family Mormyridae constitute a group of teleost fishes, formed by approximately 200 species , all endemic to African riverine and, partially, lacustrine systems. As their name suggests, they are able to produce only weak electric fields, in the order of millivolts to few volts, which are not used for predation or defence. The cells forming their electric organ are compressed disk-like cells commonly called electrocytes. In many species they are longitudinally stacked behind each other in order to form columns of cells embedded within tubes of isolating connective tissue. The synchronous activity of each electrocyte defines the output of the electric organ, known as Electric Organ Discharge (EOD) . Such weak pulses are mainly employed for the localization and discrimination of objects in water (active electrolocation) , for the recognition of conspecific individuals [7, 8], and in social and reproductive behaviour [9, 10].
Besides these functional roles, the EODs of Mormyrids display remarkable levels of differentiation in terms of shape and length across different species . Such differences, which can be dramatic even among closely related species, are considered to have acted as effective prezygotic isolation mechanisms, promoting thus the adaptive radiation observed in this family, facilitated by either ecological speciation [12–14] or speciation driven by sexual selection .
In all mormyrids, the adult electric organ is located in the caudal peduncle and is formed by four columns of electrocytes, two dorsal and two ventral ones. Each electrocyte is innervated by electromotoneurons originating in the spinal cord . Electric organs arise in juvenile fishes from several myomeres of the deep lateral muscle; their myogenic origin is confirmed by the presence of disorganized myofibrils within the electrocytes [16, 17].
During the last decade, several studies have investigated the genetic basis of EOD function and evolution [18–22]. Some of them have stressed the importance of duplication in a class of sodium channel genes for the origin of EOD production and diversification in weakly-electric fishes, by revealing the presence of important functional substitutions across paralogs and by discovering their differential patterns of expression between electric organ (EO) and skeletal muscle (SM) [18, 19, 22]. More recent studies however [20, 21], based on high throughput genomic technologies (e.g., SSH; RNA-Seq) have identified many differentially expressed genes, belonging to multiple functional classes (e.g., transcription factors; ion channels; sarcomeric proteins), between skeletal muscle and electric organ in several weakly-electric species (including representatives of Mormyridae). All studies conducted so far have been focusing on gene expression differences between two tissues –i.e., skeletal muscle and electric organ- within a species (cross-tissue approach), overlooking possible differences in the same tissue across different species (cross-species approach). However, tissue-specific gene expression differences across several species might underlie important phenotypic differences  which, in the case of the electric organ of mormyrid fishes, could explain the species-specific variability of EODs.
The aim of the present work is twofold; first we aim at exploring the differential patterns of gene expression between skeletal muscle and electric organ (cross-tissue comparison) in adult specimens belonging to the mormyrid genus Campylomormyrus (C. compressirostris and C. tshokwe). We focus then, on the identification of the differentially expressed genes, that are in common between the two species, and that might be responsible for the functional differences between the two tissues, and compare them to the results obtained by previous studies. The second task is to find differences in gene expression among three mormyrid species (C.compressirostris, C.tshokwe, and the outgroup species Gnathonemus petersii; Fig. 1) for the skeletal muscle and electric organ separately (cross-species comparison), and identify genes potentially related to phenotypic differences in EOD shape and duration.
Transcriptome sequencing and assembly
Sequencing of the twelve cDNA libraries produced a total amount of 371,043,357raw read pairs, resulting in 330,595,546 quality-filtered read pairs (89.1 %); see Additional file 1 for per library sequencing statistics. Trinity assembly resulted in 260,598 and 369,030 contigs for C. compressirostris and C.tshokwe cross-tissue transcriptomes respectively (Table 1); 357,832 and 399,878 contigs were obtained for the SM and EO cross-species assemblies respectively (Table 2). Contigs were then compared to the Danio rerio proteome, retrieving 18,458 and 19,363 unique proteins for C. compressirostris and C.tshokwe respectively; of these retrieved matches, 7971 (43.1 %) and 8993 (46.4 %) hits corresponded to full or nearly full-length coding sequences (Fig. 2a). For the cross-species assemblies, 20,023 and 20,352 contigs, for the SM and EO respectively, matched unique proteins in the D. rerio proteome, with 8662 (43.3 %) and 8768 (43.1 %) hits corresponding to full or nearly full-length coding sequences (Fig. 2b).
Differential Expression (DE) analysis
After transcript quantification with RSEM and DE-analysis with edgeR, 1313 transcripts resulted to be differentially expressed between EO and SM in C. compressirostris (356 up-regulated in EO and 957 down-regulated in EO) and 1002 in C. tshokwe (594 up-regulated in EO and 408 down-regulated in EO). Of all differentially expressed transcripts, 271 resulted to be shared between the two species (97 up-regulated in EO and 174 down-regulated in EO) (Fig. 3).
In order to obtain an initial overview of transcriptome-wide gene expression patterns, we performed a principal component analysis on the expression levels of all assessed transcripts with non-zero levels in both assemblies (see methods for details). The results clearly separate the data according to tissue rather than species (Fig. 4). A distance matrix was then obtained from the same dataset and a neighbour-joining gene expression tree was built (see methods) in order to analyse global evolutionary trends in more detail. The obtained tree (Fig. 5) shows a clustering pattern for the EO data where species characterized by similar EODs (G. petersii, C. compressirostris) are grouped together, whereas the species with a rather different EOD (C. tshokwe) forms an isolated cluster. SM data, on the other hand, do not seem to form any particular clustering scheme.
After DE analysis, 166 and 950 genes resulted to be differentially expressed across the three analysed species for SM and EO, respectively. Within the skeletal muscle 78 genes are up-regulated in G. petersii, 16 in C. compressirostris, 17 in C. tshokwe, and 55 in C. compressirostris and C. tshokwe together (Fig. 6a). As far as the electric organ is concerned, 232 genes are up-regulated in G. petersii, 87 in C. compressirostris, and 631 in C. tshokwe (Fig. 7a).
In order to identify over-represented functional categories and pathways within the C. compressirostris and C. tshokwe transcriptomes, the sets of up- and down-regulated genes in the electric organ were subject to an enrichment analysis involving Gene Ontology categories, as well as Reactome and KEGG biological pathways (Fig. 8). The number of retrieved categories/pathways is proportional to the amount of differentially expressed genes present in each analysed set; with fewer categories for the EO up-regulated genes compared to the EO down-regulated genes in C. compressirostris, and a similar number of categories identified between the two sets in C. tshokwe. For both species the most represented categories within the EO up-regulated genes are related to ion channel transport (e.g., sodium ion transport), multicellular organismal development, development and patterning of nerves (branching morphogenesis of a nerve; semaphorin interactions; axon guidance). The genes found to be down-regulated in the EO regard mainly functional classes like: metabolic pathways specific to muscle tissue (oxidative phosphorylation, pyruvate metabolism, calcium signalling pathway), and muscle specific categories (muscle cell differentiation, striated muscle contraction, cardiac muscle contraction).
Given the information provided by the category-based functional annotation, and to better understand the functional differences in terms of gene expression between EO and SM, independently from the species analysed, a literature search was performed on the shared set of differentially expressed genes between C. compressirostris and C. tshokwe. For each gene, phenotypic information consequent to its mis-expression (e.g., via knockdowns or non-sense mutations) in D. rerio was retrieved from the “Zebrafish Model Organism Database” (ZFIN; http://zfin.org/). All shared genes were divided into five “general” functional classes, which synthesize the categories reported in Fig. 8. The chosen categories are: “electrical activity” (genes responsible for the differential accumulation and transfer of ions across the plasma membrane), “muscular activity” (genes important for keeping a functional muscle phenotype), “metabolism” (genes involved in metabolic pathways), “transcription factors” (genes regulating gene expression) and “signal transduction” (molecules involved in signalling pathways) (Tables 2, 3, 4, 5 and 6). Many of the genes present in the category “electrical activity” are up-regulated in the EO (Table 2), they include genes coding for Na+/K+ pumps (atp1a2a), voltage-gated sodium (scn4aa) and potassium channels (kcnq5a) and cholinergic receptors (chrna7). However, other voltage-gated ion channels result to be down-regulated in the EO (kcna3, cacna2d2). There are then two members of the subfamily J of inwardly-rectifying potassium channels that show distinct patterns of expression, with one member (kcnj9) up-regulated in EO and the other (kcnj12) down-regulated. All the genes included in the class “muscular activity” are down-regulated in the EO (Table 3). As far as the “metabolism” genes are concerned (Table 4), most of the EO up-regulated transcripts are involved in the metabolism of fatty acids, glycerol, and phospholipids (e.g., acsl3b, gdpd4a, cds1), whereas the down-regulated transcripts are more involved in muscle-specific, energy production processes, like glycolysis (aldoab) and gluconeogenesis (gpib). Among transcription factors (Table 5), two of the four known myogenic factors (transcription factors that activate the expression of sarcomeric proteins), are down-regulated in the EO (myog, myf6), while the other two (myoD, myf5) do not show significant differences in expression between the two tissues. Two basic helix-loop-helix (bHLH) transcription factors (hey1, hes6) and one co-factor (her6) are up-regulated in the electric organ. Two myocyte enhancer factors (mef2aa, mef2b) show high levels of expression in the EO, whereas two regulators of SM cell proliferation (six1b, six4b) are lowly expressed in the EO. Most of the EO up-regulated genes involved in signal transduction (Additional file 9) belong to the G-protein coupled receptor (GPCR) signalling pathway (e.g., arhgef7a, arhgef7b, gpr22) and to the fibroblast growth-factor receptor (FGFR) signalling pathway (e.g., fgf8a, kal1b)
The two sets of differentially expressed genes identified for SM and EO across the three analysed species were each partitioned into four sub-clusters with related expression patterns (see the methods section for details); each sub-cluster was then subjected to an enrichment analysis like the one described in the previous paragraph and in the methods section. Of the analysed sub-clusters for the SM dataset, one out of four showed significantly enriched terms, all related to nucleotides metabolic processes (Fig. 6b). Conversely, three out of four sub-clusters were significantly enriched in functional categories for the EO dataset (Fig. 7b). The most representative enriched functional categories are: glutamate receptor activity (sub-cluster 1); TCA cycle and fatty acid metabolism (sub-cluster 2); ion transport, neuronal system, and striated muscle contraction (sub-cluster 4).
For each of the analysed sub-clusters genes with known phenotypic effect in D. rerio or H. sapiens are reported in Table 7.
Functional annotation of the 271 differentially expressed genes that are shared between C. compressirostris and C. tshokwe has revealed marked differences within several functional categories, which are probably critical in determining the observed phenotypic differences between the electric organ and the skeletal muscle. Below, the functional implications of the differentially expressed genes, in the light of what is known from other fish models, are discussed.
The up-regulation of the atp1a2a gene is explained by the fact that its product, the Na+/K+ ATP-ase, is fundamental for keeping the electrochemical gradient across the plasma membrane. Over-expression of this gene was already observed in the mormyrid fish Brienomyrus brachyistius , as well as in several species of south-American weakly-electric fishes (Gymnotiformes) . Voltage-gated ion channels, on the other hand, are important for dissipating the electric potential generated by the ATP-ases and therefore for producing an EOD in response to an action potential. In the electric organ of the analyzed species, one gene coding for a voltage-gated sodium channel (scn4aa) is highly expressed in the electric organ; differential expression of this gene and of its paralog (scn4ab) between EO and SM was demonstrated by Zakon et al.  for mormyrid and gymnotiform fishes, and suggest the role of gene duplication followed by neo-functionalization as a main driver for the evolution of electric communication . Other over-expressed genes that increase cell excitability are the potassium channels kcnq5a and kcnj9. The latter belongs to the family of inwardly rectifying potassium channels, a class of ion channels that favour the influx of K+ ions in the cell; up-regulation of members of this family was observed in the EO of the Electric eel (E. electricus) .
Repression of muscular phenotype in the EO
Many of the differentially expressed transcription factors retrieved in this study are fundamental for the regulation of myogenic development. In particular, we have found that two bHLH transcription factors: hey1 and hes6, in co-operation with her6, are up-regulated in the EO; these factors are known to negatively regulate the expression of myogenic factors in several model organisms [25, 26], including electric fish . Two of the four known myogenic regulatory factors (MRFs: myog, myf6) are down-regulated in the EO, both genes are fundamental for muscle development and differentiation ; in particular, knock-down experiments on myf6 have shown the degradation of posterior somites in D. rerio , the region where the adult EO originates . Another gene important for muscle development is the myocyte enhancer factor mef2aa. Unlike MRFs, this gene is up-regulated in the electric organ of the two species analysed here, as well as in other electric fish species [20, 21], and it is also important for the correct development of posterior somites in D. rerio . The concerted activity of transcriptional repressors and co-repressors of the myogenic program results in the down-regulation of genes coding for muscle specific proteins (Table 4), which finally determine the non-muscle characteristics of the EO like: i) the presence of few, non-contractile, myofibrils  (e.g., tcap ); ii) loss of calcium compartmentalization activity (e.g., atp2a1, atp2a2, casq1a ); and iii) decoupling of the excitation-contraction process (e.g., ryr1, stac3, jph1 [32–34]).
The observed differences in terms of gene expression between EO and SM suggest that the metabolic machinery of the electric organ could be mainly devoted to the production and turn-over of membrane structures. Indeed, many of the metabolism-related genes up-regulated in the EO are involved in the metabolism of fatty acids (acsbg2, acsl3b, cpt2), glycerophospholipids (cds1, gdpd4a, gdpd5a, gdpd5b), and cholesterol (idi1). On the other hand, most of the SM up-regulated genes are involved in typical processes of muscle metabolism like: glycolysis (aldoab, glo1, me3); gluconeogenesis (gpib, pgam2); and aminoacids metabolism (acy1, ckma, ckmt2a).
The grouping pattern emerging from the principal component analysis, where expression levels tend to group in an “organ-wise” rather than a “species-wise” fashion, is putatively due to the fact that expression levels are conserved for the same organ across different species for functional reasons. A similar pattern was already observed for more tissues across broader phylogenetic distances .
The clustering scheme obtained from the neighbour-joining analysis for the EO data might be indicative of the observed differences in terms of EOD among the three species, which may be reflected in the expression levels of a conspicuous part of the EO transcriptome. Previous studies [23, 35] have revealed that, for most tissues, gene expression levels tend to accumulate over evolutionary time, such that more closely related species have more similar expression levels. However, for tissues characterized by increased levels of adaptation (e.g., testis and liver in mammals), expression trees tend to group according to phenotypic similarity .
The results of the enrichment analysis conducted on the expression clusters for the EO have revealed interesting results. In particular, sub-cluster 2 and sub-cluster 4 (Fig. 7b) are enriched in terms which might underlie the observed EOD differences across the three species; both sub-clusters are characterized by genes which are mainly up-regulated in the EO of C. tshokwe. The terms relative to sub-cluster 2 are all related to metabolic pathways, the metabolism of fatty acids in particular. Many of the observed genes are involved in the production and turnover of cell membranes (e.g., smpd1, oxct1a, mlycd, cpt2). Sub-cluster 4 is mainly characterized by genes involved in ion transport and neuronal function; of particular importance here are sodium/potassium ATPases (atp1a1a, atp1a1b, atp1b2a), as their over-expression in the EO of C. tshokwe might explain the higher amplitude observed in its EOD . Other genes which may potentially influence EOD features are potassium channels (kcnq5b, kcnma1a, kcnk2a, kcnj11).
The cross-tissue analysis of differentially expressed genes between skeletal muscle and electric organ in two species of African weakly-fishes suggests that: i) the loss of contractile activity and the decoupling of the excitation-contraction processes are reflected by the down-regulation of the corresponding genes in the electric organ; ii) the metabolic activity of the EO might be specialized towards the production and turn-over of membrane structures; iii) several ion channels are highly expressed in the EO in order to increase excitability; iv) several myogenic factors might be down-regulated by transcription repressors in the EO.
The cross-species analysis has revealed that the EO transcriptome is more variable in terms of gene expression levels across species than the SM transcriptome. The functional annotation indicates that the most diverging functional classes across species in the EO include “metabolism of fatty acids” and “ion transport”.
In order to better understand the role played by the differentially expressed gene identified in this study, the onset of transgenic experiments (e.g., knockdown) will be necessary either in fully established model organisms (D. rerio), or in emerging models for electric fish (E. electricus).
Several specimens of C. compressirostris, C. tshokwe and G. petersii were collected in the wild during a sampling campaign conducted at the Congo River rapids south of Brazzaville (Republic of the Congo, August/September 2012). For the present study, adult female specimens for each species were selected, kept in captivity for a maximum period of 2 weeks and then euthanized for tissue sample collection. Gender and sexual maturity were assessed after dissection by checking for the presence of mature ovaries in the selected specimens. Electric organ and skeletal muscle tissue samples were dissected from the caudal peduncle and the posterior trunk musculature respectively and immediately transferred into RNAlater® (Life Technologies). The research followed internationally recognized guidelines and applicable national legislation. We received ethical approval from the deputy of animal welfare of the University of Potsdam.
RNA extraction and cDNA library preparation
The dissected tissues were processed at the University of Potsdam for RNA extraction: they were first removed from the RNAlater®-containing vials; shock frozen in liquid nitrogen and then homogenized into a buffer containing guanidine isothiocyanate and β-mercaptoethanol using a Mini-beadbeater-1 (Biospec). Total RNA was extracted using the RNeasy® Mini Kit (Qiagen), RNA quality and concentration was inspected using a Fragment Analyzer™ (Advanced Analytical Technologies, Inc.).
For the present study eight cDNA libraries were selected for sequencing (one library per species per tissue; two biological replicates), each library was obtained by pooling the total RNA from a minimum of two to a maximum of 4 different individuals (see Additional file 1). The paired-end (100 nt), strand-specific cDNA libraries were prepared using the NEXTflex™ Directional RNA-Seq Kit V2 (dUTP based) (Bioo Scientific); preparation was performed in six steps: i) mRNA enrichment from total RNA via polyA selection; ii) fragmentation; iii) first and second strand syntheses; iv) A-tailing; v) adapter and barcode ligation and vi) PCR amplification. Fragment size distribution and quality was estimated using an Agilent 2100 Bioanalyzer with the High Sensitivity DNA Chip.
Transcriptome sequencing, assembly and annotation
Transcriptome sequencing was performed at the Max Delbrück Center for Molecular Medicine; the multiplexed cDNA libraries were sequenced using one lane of an Illumina HiSeq2000 sequencing system. After sequencing, the resulting raw reads were subject to five processing steps using the program Flexbar v2.4  : i) filtering reads with uncalled bases; ii) trimming of reads at 3′-end to get a minimum average Phred quality score of 20; iii) barcode detection, removal and reads separation; iv) adapter detection and removal, and v) filtering of reads shorter than 20 bp after trimming. Quality control of both raw and processed reads was performed with FastQC v0.10.1 (Babraham Bioinformatics).
The processed reads were assembled de novo (i.e., without using a reference genome) with Trinity r20131110  (kmer length = 25). Two reference transcriptomes were produced from C. compressirostris and C. tshokwe respectively, by assembling together the reads obtained from the EO and SM libraries. Combining all reads across all tissues and all biological replicates for each species (cross-tissue assembly), or across all species and all replicates for each tissue (cross-species assembly) (Fig. 9) into a single RNA-seq dataset, allows to correctly compare transcript abundances from the analysed tissues or species by aligning the short reads from each library independently onto the same set of reference transcripts (see below and Additional files 2, 3, 4, 5, 6 and 7 for more details) .
Transcriptome annotation was conducted using the stand-alone version of the blastx algorithm implemented in Blast + v2.2.29  (E-value cutoff = 10−10) against the proteome of Danio rerio (Uniprot ID = UP000000437). Likely coding sequences were extracted from Trinity transcripts using TransDecoder (http://transdecoder.github.io/) and the longest translated Open Reading Frames (ORFs) were reported (Table 1). Protein domains were searched on the PFAM database (Pfam-A.hmm available at http://pfam.xfam.org/) using HMMER v3.1b1. The retrieved ORFs were later annotated by “blasting” them against the SwissProt (http://web.expasy.org/docs/swiss-prot_guideline.html) database using the blastp algorithm. Transcripts’ completeness was assessed by computing the proportion of transcripts and ORFs that matched to full-length top hits in their respective searches using the Perl script “analyze_blastPlus_topHit_coverage.pl” (provided with Trinity) (Fig. 2 and Additional file 8) .
The use of a single reference species for annotation is sub-optimal in terms of number of retrieved orthologs if compared to other methods like iterative BLAST searches to multiple species , however the use of a long-time established model organism (D. rerio) facilitates the enrichment analysis and the identification of experimental evidence for the functional role of a given gene.
Transcript abundance quantification and DE-analysis
Short reads were individually mapped to their respective transcriptome assemblies using Bowtie v1.0.0  with default parameters. Gene expression levels were estimated using RSEM v1.2.12 . Putative transcript artifacts and lowly expressed transcripts were filtered out using the Perl script “filter_fasta_by_rsem_values.pl” (provided with Trinity).
Principal component analyses of cross-species data was performed on a matrix of expression values of 14,436 genes with non-zero values in both assemblies, using the function “prcomp” in R. The same dataset was used for building a distance matrix of Spearman’s correlation coefficients (ρ) which was then subjected to neighbour-joining tree construction with the function “nj” in R (bootstrap = 10,000). Differential expression analysis was performed using the Bioconductor package edgeR  (minimum fold change = 4, p-value cutoff = 0.001 after FDR correction). The differentially expressed transcripts were then subject to an enrichment analysis using the Cytoscape plugin ClueGO v2.1.4 [45, 46], in order to identify over-represented functional categories from the Gene Ontology (GO) database  (http://geneontology.org), as well as over-represented biological pathways from KEGG (http://www.genome.jp/kegg/) and Reactome (http://www.reactome.org/) [48, 49]. Statistical significance was assessed using a Fisher’s exact test with FDR p-value correction (≤0.05). Cross-species data were partitioned into expression clusters using a k-means algorithm (k = 4), implemented in a perl script provided with Trinity.
All the Illumina reads used for this study are available at the Sequence Read Archive (SRA; http://www.ncbi.nlm.nih.gov/sra), under the accession number SRP050174.
Bass A. Electric organs revisited: evolution of a vertebrate communication and orientation organ. In: Bullock TH, Heiligenberg W, editors. Electroreception. New York: Wiley; 1986. p. 13–70.
Kirschbaum F. Myogenic electric organ precedes the neurogenic organ in apteronotid fish. Naturwissenschaften. 1983;70:205–7.
Moller P. Electric organs. In: Moller P, editor. Electr fishes hist behav. London: Chapman & Hall; 1995. p. 385–402.
Daget J, Gosse JP. Thys van den audenaerde DFE: check list of the freshwater fishes of africa = catalogue Des poissons D’eau douce d’Afrique, vol. 1. Paris/Tervuren: ORSTOM/MRAC; 1984.
Lissmann HW. On the function and evolution of electric organs in fish. J Exp Biol. 1958;35:156–91.
Lissmann HW, Machin KE. The mechanism of object location in gymnarchus niloticus and similar fish. J Exp Biol. 1958;35:451–86.
Feulner PGD, Plath M, Engelmann J, Kirschbaum F, Tiedemann R. Electrifying love: electric fish use species-specific discharge for mate recognition. Biol Lett. 2009;5:225–8.
Kramer B, Kuhn B. Species recognition by the sequence of discharge intervals in weakly electric fishes of the genus Campylomormyrus (Mormyridae, Teleostei). Anim Behav. 1994;48:435–45.
Bratton BO, Kramer B. Patterns of the electric organ discharge during courtship and spawning in the mormyrid fish, Pollimyrus isidori. Behav Ecol Sociobiol. 1989;24:349–68.
Crawford JD. Sex recognition by electric cues in a sound-producing mormyrid fish, pollimyrus isidori. Brain Behav Evol. 1991;38:20–38.
Alves-Gomes J, Hopkins CD. Molecular insights into the phylogeny of mormyriform fishes and the evolution of their electric organs. Brain Behav Evol. 1997;49:343–51.
Tiedemann R, Feulner P, Kirschbaum F. Electric organ discharge divergence promotes ecological speciation in sympatrically occurring African weakly electric fish (campylomormyrus). In: Evol action. Berlin Heidelberg: Springer; 2010. p. 307–21.
Feulner PGD, Kirschbaum F, Mamonekene V, Ketmaier V, Tiedemann R. Adaptive radiation in African weakly electric fish (Teleostei: Mormyridae: Campylomormyrus): a combined molecular and morphological approach. J Evol Biol. 2007;20:403–14.
Feulner P, Kirschbaum F, Tiedemann R. Adaptive radiation in the Congo river: an ecological speciation scenario for African weakly electric fish (teleostei; mormyridae; campylomormyrus). J Physiol. 2008;102:340–6.
Arnegard ME, McIntyre PB, Harmon LJ, Zelditch ML, Crampton WGR, Davis JK, et al. Sexual signal evolution outpaces ecological divergence during electric fish species radiation. Am Nat. 2010;176:335–56.
Szabo T. Development of the electric organ of mormyridae. Nature. 1960;188:760–2.
Denizot JP, Kirschbaum F, Westby GWM, Tsuji S. On the development of the adult electric organ in the mormyrid fish Pollimyrus isidori (with special focus on the innervation). J Neurocytol. 1982;11:913–34.
Zakon HH, Lu Y, Zwickl DJ, Hillis DM. Sodium channel genes and the evolution of diversity in communication signals of electric fishes: convergent molecular evolution. Proc Natl Acad Sci U S A. 2006;103:3675–80.
Arnegard ME, Zwickl DJ, Lu Y, Zakon HH. Old gene duplication facilitates origin and diversification of an innovative communication system--twice. Proc Natl Acad Sci U S A. 2010;107:22172–7.
Gallant JR, Hopkins CD, Deitcher DL. Differential expression of genes and proteins between electric organ and skeletal muscle in the mormyrid electric fish Brienomyrus brachyistius. J Exp Biol. 2012;215:2479–94.
Gallant JR, Traeger LL, Volkening JD, Moffett H, Chen PH, Novina CD, et al. Genomic basis for the convergent evolution of electric organs. Science (80-). 2014;344:1522–5.
Novak AE, Jost MC, Lu Y, Taylor AD, Zakon HH, Ribera AB. Gene duplications and evolution of vertebrate voltage-gated sodium channels. J Mol Evol. 2006;63:208–21.
Brawand D, Soumillon M, Necsulea A, Julien P, Csárdi G, Harrigan P, et al. The evolution of gene expression levels in mammalian organs. Nature. 2011;478:343–8.
Traeger LL, Volkening JD, Moffett H, Gallant JR, Chen P-H, Novina CD, et al. Unique patterns of transcript and miRNA expression in the South American strong voltage electric eel (Electrophorus electricus). BMC Genomics. 2015;16:243.
Fischer A, Gessler M. Delta-Notch--and then? Protein interactions and proposed modes of repression by Hes and Hey bHLH factors. Nucleic Acids Res. 2007;35:4583–96.
Buas MF, Kabak S, Kadesch T. The Notch effector Hey1 associates with myogenic target genes to repress myogenesis. J Biol Chem. 2010;285:1249–58.
Hinits Y, Osborn DPS, Hughes SM. Differential requirements for myogenic regulatory factors distinguish medial and lateral somitic, cranial and fin muscle fibre populations. Development. 2009;136:403–14.
Wang Y-H, Li C-K, Lee G-H, Tsay H-J, Tsai H-J, Chen Y-H. Inactivation of zebrafish mrf4 leads to myofibril misalignment and motor axon growth disorganization. Dev Dyn. 2008;237:1043–50.
Wang Y, Qian L, Dong Y, Jiang Q, Gui Y, Zhong TP, et al. Myocyte-specific enhancer factor 2A is essential for zebrafish posterior somite development. Mech Dev. 2006;123:783–91.
Zhang R, Yang J, Zhu J, Xu X. Depletion of zebrafish Tcap leads to muscular dystrophy via disrupting sarcomere-membrane interaction, not sarcomere assembly. Hum Mol Genet. 2009;18:4130–40.
Ebert AM, Hume GL, Warren KS, Cook NP, Burns CG, Mohideen MA, et al. Calcium extrusion is critical for cardiac morphogenesis and rhythm in embryonic zebrafish hearts. Proc Natl Acad Sci U S A. 2005;102:17705–10.
Jurynec MJ, Xia R, Mackrill JJ, Gunther D, Crawford T, Flanigan KM, et al. Selenoprotein N is required for ryanodine receptor calcium release channel activity in human and zebrafish muscle. Proc Natl Acad Sci U S A. 2008;105:12485–90.
Dowling JJ, Arbogast S, Hur J, Nelson DD, McEvoy A, Waugh T, et al. Oxidative stress and successful antioxidant treatment in models of RYR1-related myopathy. Brain. 2012;135(Pt 4):1115–27.
Horstick EJ, Linsley JW, Dowling JJ, Hauser MA, McDonald KK, Ashley-Koch A, et al. Stac3 is a component of the excitation-contraction coupling machinery and mutated in Native American myopathy. Nat Commun. 2013;4:1952.
Perry GH, Melsted P, Marioni JC, Wang Y, Bainer R, Pickrell JK, et al. Comparative RNA sequencing reveals substantial genetic variation in endangered primates Comparative RNA sequencing reveals substantial genetic variation in endangered primates. Genome Res. 2012;22:602–10.
Paul C, Mamonekene V, Vater M, Feulner PGD, Engelmann J, Tiedemann R, et al. Comparative histology of the adult electric organ among four species of the genus Campylomormyrus (Teleostei: Mormyridae). J Comp Physiol A. 2015;201(4):357–74.
Dodt M, Roehr JT, Ahmed R, Dieterich C. FLEXBAR-flexible barcode and adapter processing for next-generation sequencing platforms. Biology (Basel). 2012;1:895–905.
Grabherr MG, Haas BJ, Yassour M, Levin JZ, Thompson DA, Amit I, et al. Full-length transcriptome assembly from RNA-Seq data without a reference genome. Nat Biotechnol. 2011;29:644–52.
Haas BJ, Papanicolaou A, Yassour M, Grabherr M, Blood PD, Bowden J, et al. De novo transcript sequence reconstruction from RNA-seq using the Trinity platform for reference generation and analysis. Nat Protoc. 2013;8:1494–512.
Camacho C, Coulouris G, Avagyan V, Ma N, Papadopoulos J, Bealer K, et al. BLAST+: architecture and applications. BMC Bioinformatics. 2009;10:421.
Tzika AC, Helaers R, Schramm G, Milinkovitch MC. Reptilian-transcriptome v1.0, a glimpse in the brain transcriptome of five divergent Sauropsida lineages and the phylogenetic position of turtles. Evodevo. 2011;2:19.
Langmead B, Trapnell C, Pop M, Salzberg SL. Ultrafast and memory-efficient alignment of short DNA sequences to the human genome. Genome Biol. 2009;10:R25.
Li B, Dewey CN. RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinformatics. 2011;12:323.
Robinson MD, McCarthy DJ, Smyth GK. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. 2010;26:139–40.
Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, et al. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 2003;13:2498–504.
Bindea G, Mlecnik B, Hackl H, Charoentong P, Tosolini M, Kirilovsky A, et al. ClueGO: a Cytoscape plug-in to decipher functionally grouped gene ontology and pathway annotation networks. Bioinformatics. 2009;25:1091–3.
Ashburner M, Ball CA, Blake JA, Botstein D, Butler H, Cherry JM, et al. Gene ontology: tool for the unification of biology. Gene Ontology Consortium Nat Genet. 2000;25:25–9.
Kanehisa M. KEGG: Kyoto Encyclopedia of Genes and Genomes. Nucleic Acids Res. 2000;28:27–30.
Joshi-Tope G, Gillespie M, Vastrik I, D’Eustachio P, Schmidt E, de Bono B, et al. Reactome: a knowledgebase of biological pathways. Nucleic Acids Res. 2005;33(Database issue):D428–32.
Doganli C, Kjaer-Sorensen K, Knoeckel C, Beck HC, Nyengaard JR, Honoré B, et al. The α2Na+/K + −ATPase is critical for skeletal and heart muscle function in zebrafish. J Cell Sci. 2012;125(Pt 24):6166–75.
Friedrich T, Lambert AM, Masino MA, Downes GB. Mutation of zebrafish dihydrolipoamide branched-chain transacylase E2 results in motor dysfunction and models maple syrup urine disease. Dis Model Mech. 2012;5:248–58.
Chen Y-H, Pai C-W, Huang S-W, Chang S-N, Lin L-Y, Chiang F-T, et al. Inactivation of Myosin binding protein C homolog in zebrafish as a model for human cardiac hypertrophy and diastolic dysfunction. J Am Heart Assoc. 2013;2:e000231.
Seguchi O, Takashima S, Yamazaki S, Asakura M, Asano Y, Shintani Y, et al. A cardiac myosin light chain kinase regulates sarcomere assembly in the vertebrate heart. J Clin Invest. 2007;117:2812–24.
Postel R, Vakeel P, Topczewski J, Knöll R, Bakkers J. Zebrafish integrin-linked kinase is required in skeletal muscles for strengthening the integrin-ECM adhesion complex. Dev Biol. 2008;318:92–101.
Li H, Zhong Y, Wang Z, Gao J, Xu J, Chu W, et al. Smyd1b is required for skeletal and cardiac muscle function in zebrafish. Mol Biol Cell. 2013;24:3511–21.
Pan W, Pham VN, Stratman AN, Castranova D, Kamei M, Kidd KR, et al. CDP-diacylglycerol synthetase-controlled phosphoinositide availability limits VEGFA signaling and vascular morphogenesis. Blood. 2012;120:489–98.
Molt S, Bührdel JB, Yakovlev S, Schein P, Orfanos Z, Kirfel G, et al. Aciculin interacts with filamin C and Xin and is essential for myofibril assembly, remodeling and maintenance. J Cell Sci. 2014;127(Pt 16):3578–92.
Bitomsky N, Conrad E, Moritz C, Polonio-Vallon T, Sombroek D, Schultheiss K, et al. Autophosphorylation and Pin1 binding coordinate DNA damage-induced HIPK2 activation and cell death. Proc Natl Acad Sci U S A. 2013;110:E4203–12.
Gyda M, Wolman M, Lorent K, Granato M. The tumor suppressor gene retinoblastoma-1 is required for retinotectal development and visual function in zebrafish. PLoS Genet. 2012;8:e1003106.
Filippi A, Jainok C, Driever W. Analysis of transcriptional codes for zebrafish dopaminergic neurons reveals essential functions of Arx and Isl1 in prethalamic dopaminergic neuron development. Dev Biol. 2012;369:133–49.
Liu Y, Semina EV. pitx2 Deficiency results in abnormal ocular and craniofacial development in zebrafish. PLoS One. 2012;7:e30896.
Bessarab DA, Chong S-W, Srinivas BP, Korzh V. Six1a is required for the onset of fast muscle differentiation in zebrafish. Dev Biol. 2008;323:216–28.
Funding was provided by the University of Potsdam and the Leibniz-SAW-project GENART. Computational resources were kindly provided by the Unit of Bioinformatics (J. Selbig, S. Hartmann) at the University of Potsdam.
The authors are immensely grateful to Dr. Victor Mamonekene (Université Marien Ngouabi, Brazzaville) for the logistic and scientific support provided during their stay in Congo.
The authors declare that they have no competing interests.
FL participated in specimen sampling, extracted the RNA, carried out the bioinformatic analyses and drafted the manuscript. FK participated in specimen sampling, determined sex and species for the collected samples and participated in manuscript drafting and supervision. IK produced the cDNA libraries and prepared them for sequencing. CD participated in experimental design and provided sequencing resources. RT participated in specimen sampling, conceived and supervised the study. All authors read and provided advice on the manuscript.
Sequencing statistics. For each of the eight produced libraries we report: the number of pooled individuals, the number of pre- and post- processing read pairs, and the percent of retained reads after quality filtering. (DOCX 17 kb)
Blast results of C. compressirostris transcriptome assembly. Tabular blast report showing the results of the blastx comparison between the transcriptome of C. compressirostris and the proteome of D. rerio. (TXT 14648 kb)
Blast results of C. tshokwe transcriptome assembly. Tabular blast report showing the results of the blastx comparison between the transcriptome of C. tshokwe and the proteome of D. rerio. (TXT 19887 kb)
Changes in gene expression between EO and SM for C. compressirostris transcriptome. Results of the edgeR DE-analysis. The table reports the fold-change in the expression levels and relative uncorrected and corrected p-values between EO and SM for all genes in the C. compressirostris transcriptome. (TXT 5878 kb)
Changes in gene expression between EO and SM for C. tshokwe transcriptome. Results of the edgeR DE-analysis. The table reports the fold-change in the expression levels and relative uncorrected and corrected p-values between EO and SM for all genes in the C. tshokwe transcriptome. (TXT 7336 kb)
Transcripts’ length distribution of the C. compressirostris Trinity assembly. (PDF 24 kb)
Transcripts’ length distribution of the C. tshokwe Trinity assembly. (PDF 5 kb)
Distribution of length coverage between retrieved ORFs and reference database (SwissProt). (PDF 44 kb)
List of differentially expressed genes belonging to the class “signal transduction”, obtained from the cross-tissue comparison. (DOCX 39 kb)