Skip to main content

Advertisement

The transcriptional correlates of divergent electric organ discharges in Paramormyrops electric fish

Abstract

Background

Understanding the genomic basis of phenotypic diversity can be greatly facilitated by examining adaptive radiations with hypervariable traits. In this study, we focus on a rapidly diverged species group of mormyrid electric fish in the genus Paramormyrops, which are characterized by extensive phenotypic variation in electric organ discharges (EODs). The main components of EOD diversity are waveform duration, complexity and polarity. Using an RNA-sequencing based approach, we sought to identify gene expression correlates for each of these EOD waveform features by comparing 11 specimens of Paramormyrops that exhibit variation in these features.

Results

Patterns of gene expression among Paramormyrops are highly correlated, and 3274 genes (16%) were differentially expressed. Using our most restrictive criteria, we detected 145–183 differentially expressed genes correlated with each EOD feature, with little overlap between them. The predicted functions of several of these genes are related to extracellular matrix, cation homeostasis, lipid metabolism, and cytoskeletal and sarcomeric proteins. These genes are of significant interest given the known morphological differences between electric organs that underlie differences in the EOD waveform features studied.

Conclusions

In this study, we identified plausible candidate genes that may contribute to phenotypic differences in EOD waveforms among a rapidly diverged group of mormyrid electric fish. These genes may be important targets of selection in the evolution of species-specific differences in mate-recognition signals.

Background

Understanding the genomic basis of phenotypic diversity is a major goal of evolutionary biology [1]. Adaptive radiations and explosive diversification of species [2] are frequently characterized by interspecific phenotypic differences in divergence of few, hypervariable phenotypic traits [3,4,5,6]. Such systems offer exceptional advantages to study the genomic bases of phenotypic diversity: they can provide replication under a controlled phylogenetic framework [7], and couple ample phenotypic differentiation with relatively “clean” genomic signals between recently diverged species [8]. Study of the genomic mechanisms underlying hypervariable phenotypic traits has identified, in some cases, relatively simple genetic architectures [9,10,11,12,13]. More often, the genetic architecture underlying such traits can be complex and polygenic [14,15,16,17]. It has long been recognized that changes in gene expression can affect phenotypic differences between species [18], and RNA-seq based approaches have greatly facilitated the study of this relationship [19]. A growing number of studies have examined differences in gene expression in phenotypic evolution (e.g., [19,20,21,22,23,24,25,26,27]). While these studies do not investigate mutational causes, analysis of differential gene expression (DGE) can be a useful approach in examining the genomic basis of divergent phenotypes.

African weakly electric fish (Teleostei: Mormyridae) are among the most rapidly speciating groups of ray-finned fishes [28, 29]. This is partly due to the diversification of the genus Paramormyrops [30, 31] in the watersheds of West-Central Africa, where more than 20 estimated species [32] have evolved within the last 0.5–2 million years [30]. Extensive evidence has demonstrated that electric organ discharges (EODs) exhibit little intraspecific variation, yet differ substantially among mormyrid species [33,34,35]. This pattern is particularly evident in Paramormyrops [30, 36], in which EOD waveforms evolve much faster than morphology, size, and trophic ecology [37].

Mormyrid EODs are a behavior with a dual role in electrolocation [38, 39] and intraspecific communication [40, 41]. EOD waveforms vary between species principally in terms of their complexity, polarity, and duration [30, 42], and all three dimensions of variation are evident among Paramormyrops (Fig. 1). Furthermore, recent discoveries of intraspecific polymorphism in EOD waveform in P. kingsleyae [43] and polarity among P. sp. ‘magnostipes’ [35] present a unique opportunity to study the genomic basis of phenotypic traits within a rapidly diverging species group.

Fig. 1
figure1

Electric organ discharge (EOD) diversity and electric organ anatomy in Paramormyrops. EOD traces from specimens in this study and representative parasagittal sections of the five Paramormyrops operational taxonomic units (OTUs) considered in this study. 200x magnification on P. kingsleyae EODs reveals a P0 phase on triphasic EODs only. Individuals with triphasic EODs all have penetrations, whereas individuals with biphasic EODs do not. OTUs with ‘inverted’ polarity triphasic EODs have large penetrations compared to OTUs with normal polarity triphasic EODs. Ant. = anterior, C = connective tissue septa, N = nerve, NPp = non-penetrating, posteriorly innervated stalks, M = microstalklets (profusely branched stalks), P = penetrations, Pa = penetrating, anteriorly innervated stalks, Post. = posterior, S = stalks. Image from P. kingsleyae biphasic originally appeared in reference 43

EODs have a well-understood morphological (Fig. 1) and neurophysiological basis [44, 45]. EODs are generated by specialized cells (electrocytes) that constitute the electric organ (EO), located in the caudal peduncle [46]. Mormyrid EOs are comprised of 80–360 electrocytes [34], and an individual EOD is produced when the electrocytes discharge synchronously. EODs are multiphasic because they result from action potentials produced by two excitable membranes: the two large phases of the EOD, called P1 and P2, are produced by spikes generated by the posterior and anterior electrocyte faces, respectively [47]. There is a relationship between EODs of longer duration and increased surface membrane area [48], likely mediated at least in part by an increase in membrane capacitance [49, 50]. The duration of EODs is highly variable within mormyrids-- some EODs are extremely long (> 15 ms) and others are very brief (0.2 ms) [32].

Within the Mormyridae, triphasic EODs evolved early from biphasic EODs; however, there have been multiple parallel reversions to biphasic EODs across mormyrids and within Paramormyrops [36, 43]. Triphasic (P0-present) EODs are produced by electrocytes that are innervated on the anterior face and have penetrating stalks (Pa, P-type), whereas biphasic (P0-absent) EODs are produced by electrocytes innervated on the posterior face and lack penetrating stalks (NPp, N-type) (for more details see [42, 43, 47, 48, 51, 52]). We refer to triphasic EODs as more ‘complex’ than biphasic EODs. In some cases, triphasic EODs display an unusually large P0 phase, which gives the appearance of an ‘inverted’ polarity. This is exemplified by the type I EODs of P. sp. ‘magnostipes’ (Fig. 1) [35]. The number [47] and diameter [34, 43] of stalk penetrations are positively correlated with the magnitude of P0. We refer to individuals with large penetrations as ‘inverted’ polarity and individuals with small penetrations as ‘normal’ polarity.

Recent studies in mormyrids [53,54,55,56,57] have adopted a candidate gene approach to examine the molecular basis of variation in EOD duration on macroevolutionary scales, implicating voltage gated sodium channels (e.g. scn4aa) and potassium channels (e.g. kcna7a) as key targets of selection during EOD evolution. Beyond this recent attention to ion channels, several studies have described the importance of structural differences between EOs as an important component of EOD variation [43, 48, 50]. In this study, we took a transcriptome-wide approach to characterizing the molecular basis of electric signal diversity in Paramormyrops species divergent for EOD complexity, duration and polarity. We used RNA-sequencing to comprehensively examine DGE in the adult EOs of five Paramormyrops operational taxonomic units (OTUs), leveraging a recently sequenced and annotated genome assembly from the species P. kingsleyae (N-type) [58], and identify gene expression correlates of each of the three main EOD waveform features of electric signal diversity in Paramormyrops. Our results emphasize genes that influence the shape and structure of the electrocyte cytoskeleton, membrane and extracellular matrix (ECM) to exhibit predictable differences between Paramormyrops species with divergent EOD phenotypes.

Results

Overall results

We extracted and sequenced mRNA from EOs of 11 wild caught Paramormyrops samples from five OTUs (Table 1). Overall alignment rates of the processed reads to the Paramormyrops kingsleyae reference transcriptome ranged from 28 to 74% (> 375 million sequenced reads in total, 50% aligned), with no clear differences among OTUs (Additional file 1). On inspection, we concluded that these rates are a consequence of the presence of overrepresented sequences from rRNA, mtDNA and bacterial contamination in the RNA-seq reads.

Table 1 Phenotypic and collection information of the samples studied

We explored the data with a heatmap of pairwise correlations of gene expression for 24,960 genes across all 11 samples (Fig. 2), and carried out all possible pairwise DGE comparisons of OTUs (n = 10, Table 2). These ten comparisons detected a range of 16,420–19,273 expressed genes. Intersection of these lists resulted in a non-redundant list of 20,197 genes expressed in EO across all DGE comparisons. We found that 3274 (16%) were differentially expressed in at least one comparison, and expression patterns across all OTUs were highly correlated (Pearson’s r > 0.89, Fig. 2). Despite this, correlation values were higher within recognized OTUs, except for the P. sp. ‘magnostipes type II’ 6768 sample (Fig. 2). Thus, we did not use P. sp. ‘magnostipes type I’ vs P. sp. ‘magnostipes type II’ as the informative comparison for waveform polarity in Set A’ (see methods).

Fig. 2
figure2

Heatmap of sample by sample correlations in gene expression, and the inferred relationships among OTUs from these expression correlation values. OTU = operational taxonomic unit

Table 2 All ten possible pairwise DGE comparisons with the total number of DEG and enriched GO terms for each. Also indicated is whether each comparison is informative for contrasting each EOD feature. The phenotypes for waveform polarity can only be contrasted in comparisons where both OTUs have penetrations. Informative comparisons for each EOD feature (Set A’) are marked with an * in the column of the EOD feature they contrasted

Set A: differential expression analysis

We found between 489 and 1542 differentially expressed genes (DEGs; fold change > 4, FDR-corrected p-value < 0.001) (50–128 enriched Gene Ontology (GO) terms) in every pairwise comparison of OTUs except P. sp. ‘magnostipes type I’ vs P. sp. ‘magnostipes type II’, which had only nine DEGs with seven enriched GO terms (Table 2). Additional file 2 provides a tabular list of DEGs for each comparison, and Additional file 3 provides a tabular list of enriched GO terms for each comparison. We call Set A the non-redundant list of 3274 genes that were differentially expressed in at least one DGE comparison (Fig. 3, Set A).

Fig. 3
figure3

Diagram of how we constructed the lists of upregulated genes of Set C. DEG = differentially expressed genes, N = number of comparisons made for each set, OTU = operational taxonomic unit

We chose the phylogenetically most informative comparisons to construct Set A’, which are indicated in Table 2. We found: 507 DEG and 69 enriched GO terms comparing P. kingsleyae (N-type) vs P. sp. ‘SN3’ (EOD duration); 1322 DEG and 77 enriched GO terms comparing P. kingsleyae (P-type) vs P. sp. ‘magnostipes type I’ (waveform polarity); and 530 DEG and 75 enriched GO terms comparing P. kingsleyae (N-type) vs P. kingsleyae (P-type) (waveform complexity).

Set B: expression based clustering

For each EOD feature (n = 3), we grouped OTUs by phenotype (Table 1), and calculated normalized expression values for Set A genes. From these, we constructed Set B by selecting genes that (1) exhibit greater than four-fold difference in the average expression levels between phenotypes of each EOD feature, and (2) have within-phenotype standard deviations less than the difference between phenotype-mean expression. For the phenotypes of waveform duration, we identified a combined total of 309 DEG and 43 enriched GO terms, for waveform polarity we found 169 DEG and 14 enriched GO terms, and for waveform complexity the totals were 413 DEG and 38 enriched GO terms. Additional file 4 lists the identities of these DEG and Additional file 5 lists their enriched GO terms for all three GO ontologies.

Set C: intersection of phylogenetically informative comparisons and expression based clustering

We were motivated to obtain the DEGs and enriched GO terms that were most likely to be associated with divergent EOD phenotypes. To obtain this list, we constructed Set C, which is the intersection of Set A’ and Set B (Fig. 3, see methods). The expression profiles of the Set C genes for each EOD feature, along with the enriched GO terms for Biological Process and Cellular Component, are shown in Figs. 4-6.

Fig. 4
figure4

Set C for waveform duration. a) Expression patterns of Set C genes for the waveform duration phenotypes short EODs (purple background) and long EODs (yellow background). Samples are sorted alphabetically on the X axis. The lines connect transformed gene expression values across all samples; light-color lines represent one gene, the dark-color line is the average expression pattern of all genes. b) Gene Ontology (GO) terms for Biological Process and Cellular Component found enriched in the gene lists from (a). The X axis shows transformed p-values, the longer a bar the smaller its p-value. The direction and color of a bar indicate the phenotype in which the GO term is enriched [same color code as (a)]

Fig. 5
figure5

Set C for waveform polarity. a) Expression patterns of Set C genes for the waveform polarity phenotypes small penetrations (red background) and large penetrations (grey background). Samples are sorted alphabetically on the X axis. The lines connect transformed gene expression values across all samples; light-color lines represent one gene, the dark-color line is the average expression pattern of all genes. b) Gene Ontology (GO) terms for Biological Process and Cellular Component found enriched in the gene lists from (a). The X axis shows transformed p-values, the longer a bar the smaller its p-value. The direction and color of a bar indicate the phenotype in which the GO term is enriched [same color code as (a)]. Pen = penetrations

Fig. 6
figure6

Set C for waveform complexity. a) Expression patterns of Set C genes for the waveform complexity phenotypes triphasic (orange background) and biphasic (blue background). Samples are sorted alphabetically on the X axis. The lines connect transformed gene expression values across all samples; light-color lines represent one gene, the dark-color line is the average expression pattern of all genes. b) Gene Ontology (GO) terms for Biological Process and Cellular Component found enriched in the gene lists from (a). The X axis shows transformed p-values, the longer a bar the smaller its p-value. The direction and color of a bar indicate the phenotype in which the GO term is enriched [same color code as (a)]

Contrast of waveform duration identified 183 DEG and 39 enriched GO terms. 140 of the DEG were upregulated in the short EOD phenotype (Fig. 4 purple lines), and 43 genes were upregulated in samples with long EODs (Fig. 4 yellow lines). Contrast of waveform polarity identified 154 DEG and 17 enriched GO terms. We found 99 upregulated genes in individuals with small penetrations (Fig. 5 red lines), and 55 upregulated genes in individuals with large penetrations (Fig. 5 grey lines). Finally, contrast of waveform complexity identified 145 DEG and 20 enriched GO terms. We detected 110 upregulated genes in individuals with biphasic EODs (Fig. 6 blue lines), and 35 upregulated genes in individuals with triphasic EODs (Fig. 6 orange lines). These results are further detailed in Table 3. All genes in Set C are listed in Additional file 6 and their enriched GO terms are listed in Additional file 7.

Table 3 Total number of upregulated genes and enriched GO terms in Set C for each EOD feature, phenotype and ontology

Discussion

It has long been recognized that changes in gene expression can affect phenotypic differences between species [18], and RNA-seq has facilitated the study of this relationship [19]. The goal of this study was to determine DEGs associated with divergent EOD features within Paramormyrops. Expression patterns across all OTUs were highly correlated (Pearson’s r > 0.89, Fig. 2) and we detected differential expression of only 3274 (16%) genes between any two OTUs. Thus, a major finding of this study is that EO gene expression is overall quite similar across Paramormyrops species with divergent EODs, and relatively few genes are associated with phenotypic differences in EOD waveform between OTUs. We find this notable given observations of generally high levels of genetic distances between geographically proximate populations of the same Paramormyrops species [35, 59].

Despite the relatively small number of DEGs compared to the total number of genes expressed in the EO, we constructed our analysis to extract genes that were highly associated with particular phenotypes. Set A’ represents a formal statistical test that contrasted OTUs. Each comparison contrasted samples from OTUs that were divergent in only one EOD feature, while minimizing phylogenetic distance. The drawback of this approach is the observed differences may not reflect a general pattern across multiple OTUs, instead resulting from OTU-specific changes or confounding variables such as collection sites. Set B took the opposite approach, using information from all possible biological replicates to identify consistent gene expression patterns between EOD phenotypes, at the expense of formal statistical support and the introduction of confounding phylogenetic relationships and phenotypic heterogeneity. To balance these drawbacks, we constructed Set C, which represents genes and GO terms that are differentially expressed/enriched between closely related OTUs divergent in only one phenotypic character and that are also consistently differentially expressed/enriched among representatives with similar EOD phenotypes. As such, we focus our discussion on the results of Set C. We classified the genes in Set C into “general” functional classes, or themes; and focus our attention on the ones that relate to the known morphological underpinnings of waveform duration (Table 4), polarity (Table 5), and complexity (Table 6). These functional classes were genes related to the ECM, cation homeostasis, lipid metabolism, and cytoskeletal and sarcomeric genes.

Table 4 Selected DEG in Set C for waveform duration by “general” functional class and EOD phenotype, and highlights of their predicted function
Table 5 Selected DEG in Set C for waveform polarity, by “general” functional class and EOD phenotype, and highlights of their expected function
Table 6 Selected DEG in Set C for waveform complexity, by “general” functional class and EOD phenotype, and highlights of their expected function

Waveform duration

Several researchers have implicated the role of ion channels in the evolution of duration changes in mormyrid signals [53,54,55,56,57]. We did not find evidence of large changes in expression of sodium channels between short-duration P. sp. ‘SN3’ and other Paramormyrops species; but we detected upregulation in short EOD samples of a voltage insensitive, outwardly rectifying potassium channel (potassium channel subfamily K member 5-like) and of a regulatory subunit of a Shal-type voltage-gated potassium channel (Kv channel-interacting protein 4). Additionally, individuals with short EODs upregulate two calcium-binding proteins: parvalbumin-2, and parvalbumin-2-like. Parvalbumins are highly expressed in skeletal muscle where they sequester calcium after contraction, thus facilitating relaxation. Frequently, muscles with fast relaxation rates express higher levels of parvalbumins [60]. The upregulated parvalbumin genes we detected may somehow be related to shorter EODs by sequestering calcium at a faster rate, which could affect action potentials directly or indirectly through calcium-activated ion channels.

Previous studies have demonstrated that changes in EOD duration result from changes in electrocyte ultrastructure. The two major phases of the EOD waveform are caused by action potentials generated by the anterior and posterior faces [47]. Bennett [48] demonstrated a relationship between EOD duration and increased surface membrane area, and Bass et al. [50] showed that differences in surface area are more readily noticeable on the anterior face. Membrane surface area is increased by folding the electrocyte membrane into papillae and other tube-like invaginations [61]. Testosterone can induce increases in EOD duration in several mormyrids [49, 50, 62, 63], and it also increases membrane surface area, either particularly on the anterior face [50] or on both anterior and posterior faces [64]. A larger surface area may increase the capacitance of the membrane, thus delaying spike initiation [49, 50]. Consequently, genes involved in the synthesis of membranes could influence EOD duration.

We found the most prominent differences in gene expression between the EOD duration phenotypes in genes that code for cytoskeletal, sarcomeric, and lipid metabolism proteins (Table 4). We emphasize the last group: no lipid metabolism genes were upregulated in individuals with long EODs, whereas samples with short EODs upregulated protein EFR3 homolog B-like (a regulator of phosphatidylinositol 4-phosphate synthesis), retinoic acid receptor responder protein 3-like and HRAS-like suppressor 3 (these two catalyze hydrolysis of phosphatidylcholines and phosphatidylethanolamines), fatty aldehyde dehydrogenase-like (fatty acid metabolism), PTB domain-containing engulfment adapter protein 1-like (modulates cellular glycosphingolipid and cholesterol transport), phosphatidylinositol 3-kinase regulatory subunit gamma-like, (PI3K, which phosphorylates phosphatidylinositol), and proto-oncogene c-Fos-like (can activate phospholipid synthesis), and showed enrichment of the GO term ‘phosphatidylinositol phosphorylation.’ We hypothesize that these genes are involved in the surface proliferation of the electrocytes membranes.

Additionally, each mormyrid electrocyte stands embedded in a gelatinous mucopolysaccharide matrix (the ECM) separated from neighboring electrocytes by connective tissue septa (Fig. 1) [34], and the membrane surface invaginations are coated by the same ECM that surrounds the electrocytes [50, 61]. Hence, differences in surface invaginations could also be reflected in differences in the expression of genes whose products interact with the ECM. In individuals with long EODs, we found upregulated the genes matrilin 2 (involved in matrix assembly) and thrombospondin-4-like (mediates cell-to-matrix interactions), and enriched the GO term ‘extracellular space’; whereas those with short EODs upregulated collagenase 3-like (plays a role in the degradation of ECM proteins) and displayed enrichment of the GO terms ‘extracellular region’ and ‘external side of plasma membrane’.

Two previous studies focused on DGE between EOs in another mormyrid species adaptive radiation/explosive diversification (genus Campylomormyrus, whose EO anatomy closely resembles that of Paramormyrops). Both focused on comparisons between species with biphasic EODs but different waveform duration phenotypes: C. tshokwe (long duration) and C. compressirostris (short duration). The first study performed a canditate gene approach to quantify the expression patterns of 18 sodium and potassium homeostasis genes between the EOs of the two species [55], whereas Lamanna et al. [65] used RNA-seq to simultaneously compare gene expression between EOs of these species. While we did not observe differences in expression of any of the potassium channels reported by Nagel et al. [55], we note that Lamanna et al. [65] reported differential expression of metabolic pathways related genes, particularly fatty acid metabolism, and ion transport and neuronal function (referred to in their text as cross-species analysis (EO) subclusters 2 and 4). While we found no overlap in the identities of any specific genes in our study, we note that our analysis also detected differential expression of lipid metabolism related genes when comparing EODs of different duration.

Overall, our results identify genes that may affect EOD duration through membrane rearrangements, which could be coupled with changes in the interaction with the ECM and the expression of cytoskeletal and sarcomeric genes. Since this waveform feature is modulated by testosterone, this androgen could facilitate the study of these suggested genetic underpinnings under more rigorously controlled circumstances.

Waveform polarity

The number [47] and diameter [34, 43] of stalk penetrations are positively correlated with the magnitude of P0. This phenomenon is exemplified by P. sp. ‘magnostipes type I’, which has the largest P0 in the OTUs examined in this study, giving the EOD the appearance that it ‘inverted’ relative to other EODs. This OTU has numerous, large diameter penetrations, whereas P. kingsleyae (P-type) has relatively fewer, small diameter penetrations (Fig. 1). These large structural differences may influence the electrocyte’s connection with the surrounding ECM, and our results support this: the phenotypes of waveform polarity exhibited differences in the expression of genes that interact with the extracellular space. We found no such genes upregulated in individuals with large penetrations, whereas in samples with small penetrations we detected the enriched GO terms: ‘extracellular matrix structural constituent’ and ‘basement membrane,’ and the upregulated genes: collagen alpha-1(V) chain-like, collagen type IX alpha 2 chain, collagen alpha-4(VI) chain-like, epiphycan-like (may play a role in cartilage matrix organization), cell migration-inducing and hyaluronan-binding protein-like (mediates depolymerization of hyaluronic acid), inter-alpha-trypsin inhibitor heavy chain H5-like (although this gene is little studied, inter-alpha-trypsin inhibitors usually interact with hyaluronan), and thrombospondin 2 (mediates cell-to-matrix interactions).

OTUs with small penetrations also exhibited higher expression of genes related to cytoskeletal, sarcomeric, and lipid metabolism proteins than do individuals with large penetrations (Table 5). This includes the genes myosin light chain 3-like, desmin-like, PH and SEC7 domain-containing protein 1-like (induces cytoskeletal remodeling), CDC42 effector protein 3 (probably involved in the organization of the actin cytoskeleton), protein phosphatase 1 regulatory subunit 12B-like (regulates myosin phosphatase activity), protein kinase C alpha type-like (phosphorylates cardiac troponin T), and calponin-1-like (modulates smooth muscle contraction). Samples with large penetrations showed upregulation of the genes myosin light chain 4-like (regulatory light chain of myosin), leucine rich repeat containing 10, (may play important roles in cardiac development and/or function), and Wiskott-Aldrich syndrome protein family member 3-like (regulation of cell morphology and cytoskeletal organization).

We hypothesize that the differences in the number and diameter of penetrations that drive variation in EOD waveform polarity require changes to the electrocyte’s cytoskeletal and membrane properties. These arrangements may be necessary for the electrocytes body to adjust to the increased volume displacements imposed by larger penetrations; or alternatively, they may be a prerequisite for penetrating stalks to enlarge. Our observations support and elaborate on the hypothesis that sarcomeric proteins (which are non-contractile in mormyrid electric organs) may function as a means of cytoskeletal support and structural integrity in mormyrid electrocytes [66].

Waveform complexity

Waveform complexity refers to the number of phases present in an EOD, and mormyrid EODs vary in the presence of a small head negative phase (P0). The presence or absence of P0 in the EOD depends on the anatomical configuration of the electrocytes: P0-present (or triphasic) EODs are produced by electrocytes that are innervated on the anterior face and have penetrating stalks (Pa), whereas P0-absent (or biphasic) EODs are produced by electrocytes innervated on the posterior face and lack penetrating stalks (NPp) [42, 43, 47, 48, 51, 52]. Developmental studies of the adult EO suggest that Pa electrocytes go through a NPp stage before developing penetrations [67, 68]. This motivated the hypothesis that penetrations develop by the migration of the posteriorly innervated stalk system (NPp stage) through the edge of the electrocyte, and that the interruption of this migration represents a mechanism for Pa-to-NPp reversals [69, 70].

Our data indicates several DEGs that implicate specific cytoskeletal and ECM reorganizations between triphasic and biphasic EODs (Table 6). We observed differential expression of several genes associated with the polymerization of F-actin. In triphasic individuals, we observe upregulation of the gene capping protein regulator and myosin 1 linker 3 (CARMIL3); although this gene is little studied, its paralog CARMIL2 enhances F-actin polymerization. In contrast, the biphasic phenotype upregulated the genes protein-methionine sulfoxide oxidase mical2b-like (promotes F-actin depolymerization), transmembrane protein 47-like (may regulate F-actin polymerization), 5′-AMP-activated protein kinase subunit gamma-2-like (could remodel the actin cytoskeleton), FYVE, RhoGEF and PH domain-containing protein 4-like (regulates the actin cytoskeleton), leiomodin 2 (promotes actin polymerization, and required for normal sarcomere organization in the heart) and family with sequence similarity 110 member C (may play a role in microtubule organization). Thus, biphasic and triphasic EODs display several DEG, with potentially diverging outcomes, that influence the cellular internal structure.

We hypothesize that electrocytes with penetrating stalks (which produce triphasic EODs) require cytoskeletal arrangements to produce penetrations, perhaps related to increasing F-actin, to maintain their structural integrity. Similar to what we propose under waveform polarity, these arrangements may be necessary for the electrocyte body to adjust to the penetrations; or alternatively, they may be a prerequisite for penetrations to occur.

We also observed differential expression in proteins expressed in the ECM. In biphasic OTUs, we found the GO term ‘mucopolysaccharide metabolic process’ to be enriched, and two upregulated copies of the gene inter-alpha-trypsin inhibitor heavy chain 3, which may act as a binding protein between hyaluronan and other ECM proteins. In triphasic individuals, we found the upregulated genes epiphycan-like, which may play a role in cartilage matrix organization, and ependymin-like (ortholog to the zebrafish ependymin-like gene epdl2).

Two ependymin-like genes are among the most differentially expressed genes in the Set A’ comparison for waveform complexity P. kingsleyae (N-type) vs P. kingsleyae (P-type) (500-fold more highly expressed in P. kingsleyae (P-type), Additional file 2). Although expressed in many tissues and with little amino acid similarity, all ependymin-related proteins are secretory, calcium-binding glycoproteins that can undergo conformational changes and associate with collagen in the ECM. They have been involved in regeneration, nerve growth, cell contact, adhesion and migration processes [71]. We hypothesize that ependymin-related proteins, and potentially some of the other ECM proteins highly expressed in triphasic individuals, are part of the “fibrillar substance” that lies between the stalk and the electrocyte body in individuals with penetrating electrocytes [50]. Notably, the P. kingsleyae genome assembly, which is based on a biphasic individual, contains three paralogs of epdl2, whereas the osteoglossiform Scleropages formosus only has one, suggesting the intriguing possibility that this gene may have been duplicated in Paramormyrops or in mormyrids. Ependymin-related paralogs have been proposed as suitable targets to experimentally test gene subfunctionalization [72].

Altogether, our results for EOD waveform complexity suggest that the conformation of the cytoskeleton and the expression of proteins secreted to the ECM are important elements of the stalk penetrations, which generate triphasic EODs.

Conclusions

The widespread differential expression within Paramormyrops of calcium-related genes (Additional file 6) emphasizes a much-needed area of future research. Calcium is known to be necessary for the proper electrocyte repolarization in some gymnotiform species [73], but it may not be as important in others [74]. Few studies have addressed calcium physiology in mormyrids: calcium-related proteins have been reported as differentially expressed in EO vs skeletal muscle in Campylomormyrus [65] and in Brienomyrus brachyistius [66]. As electrocytes do not contract, calcium may act in electrocytes as an important second messenger or cofactor, participate in interactions with the ECM, and/or to contribute to the electrocyte’s electrical properties through interaction with voltage gated ion channels.

A second notable pattern in our results is the unusual degree to which mormyrid electrocytes retain expression of some sarcomeric genes, which has been noted in several studies [58, 65, 66, 75, 76]. The role these proteins serve in electrocytes is presently unknown; however, our results indicate that some are highly differentially expressed between Parmormyrops with different EOD waveforms (Tables 4-6). This strongly suggests that sarcomeric proteins could play an important role in the conformational changes required to develop and sustain penetrations.

Finally, the biochemical composition and function of the ECM in electrocytes is poorly understood. Our analysis identifies differential expression in ECM-related genes across the genus Paramormyrops, associated with each of the three EOD features studied. At least four of these genes (cell migration-inducing and hyaluronan-binding protein-like, inter-alpha-trypsin inhibitor heavy chain H5-like and two copies of inter-alpha-trypsin inhibitor heavy chain 3), distributed across two EOD features, interact with hyaluronan. Hyaluronan is a type of mucopolysaccharide and a major component of some soft tissues and fluids [77]. Therefore, we propose that hyaluronan is an important constituent of the ECM in mormyrid fish. In addition, the electrocyte-ECM interactions should be an important area of future investigation, as they are likely to influence electrocyte shape, electrical properties, and potentially the morphology of penetrations and surface membrane invaginations.

To conclude, this study examined the expression correlates of a hyper-variable phenotype in a rapidly diversified genus of mormyrid electric fish. We examined DGE between taxa exhibiting variability along three major axes of variation that characterize EOD differences within Paramormyrops and among mormyrids: duration, polarity, and complexity. We found that gene expression in EOs among closely related species is largely similar, but patterns of DGE between EOs is primarily restricted to four broad functional sets: (1) cytoskeletal and sarcomeric proteins, (2) cation homeostasis, (3) lipid metabolism and (4) proteins that interact with the ECM. Our results suggest specific candidate genes that are likely to influence the size, shape and architecture of electrocytes for future research on gene function and molecular pathways that underlie EOD variation in mormyrid electric fish.

Methods

Sample collection

We captured 11 Paramormyrops individuals from Gabon, West Central Africa in 2009: five P. kingsleyae (n = 3 N-type and n = 2 P-type), four P. sp. magnostipes’ (n = 2 Type I and n = 2 Type II), and two P. sp. ‘SN3’. Within 1–12 h of capture, individual specimens were euthanized by overdose with MS-222. The caudal peduncle was excised and skinned, and immediately immersed in RNA-later for 24 h at 4 °C, before being transferred to − 20 °C for long-term storage. As two of these species (P. sp. ‘magnostipes’, P. sp. ‘SN3’) are presently undescribed, we note that these specimens were identified by their EOD waveform, head morphology and collecting locality [30, 31, 35, 78]. All specimens, including vouchers materials, are deposited in the Cornell University Museum of Vertebrates. Collection information and the phenotypes per EOD feature of each sample are detailed in Table 1.

RNA extraction, cDNA library preparation and Illumina sequencing

Total RNA was extracted from EOs using RNA-easy Kit (Qiagen, Inc) after homogenization with a bead-beater (Biospec, Inc.) in homogenization buffer. mRNA was isolated from total RNA using a NEBNext mRNA Isolation Kit (New England Biolabs, Inc.). Libraries for RNA-seq were prepared using the NEBNext mRNA Sample Prep Master Mix Set, following manufacturer’s instructions. Final libraries after size selection ranged from 250 to 367 bp. Libraries were pooled and sequenced by the Cornell University Biotechnology Resource Center Genomics Core on an Illumina HiSeq 2000 in a 2x100bp paired end format. Raw sequence reads were deposited in the NCBI SRA (Additional file 1).

Read processing and data exploration

FastQC v0.11.3 (Babraham Bioinformatics) was used to manually inspect raw and processed reads. We used Trimmomatic v.0.32 [79] to remove library adaptors, low quality reads, and filter small reads; following the suggested settings of MacManes [80]: 2:30:10 SLIDINGWINDOW:4:5 LEADING:5 TRAILING:5 MINLEN:25. After trimming, reads from each specimen were aligned to the predicted transcripts of the NCBI-annotated (Release 100) P. kingsleyae (N-type) genome [58] using bowtie 2 v2.3.4.1 [81]. Expression quantification was estimated at the gene level using RSEM v1.3.0 [82], followed by exploration of the data with a gene expression correlation matrix based on Euclidean distances and Pearson’s correlation coefficient (for genes with read counts > 10, Trinity’s default parameters). All these steps were executed using scripts included with Trinity v2.6.6 [83, 84].

Data analysis

We began by examining DGE between all possible pairwise comparisons of OTUs (n = 10, Table 2) using edgeR v3.20.9 [85] through a script provided with Trinity (MA plots are provided in Additional File 8). We restricted our consideration of genes to those where CPM-transformed counts were > 1 in at least two samples for each comparison (edgeR default parameters). We modified this to use the function estimateDisp() instead of the functions estimateCommonDisp() and estimateTagwiseDisp(). For each comparison, we conservatively considered genes to be differentially expressed with a minimum fold change of 4 and p-value of 0.001 after FDR correction. We compiled a non-redundant list of genes that were differentially expressed in at least one comparison based on these criteria (Fig. 3, Set A).

For each of the DEGs in Set A we used TMM normalized values to compare gene expression between groups of OTUs with alternative EOD waveform phenotypes (i.e. long duration EOD vs. short duration EOD, biphasic vs. triphasic and small penetrations vs. large penetrations, see Table 1. Note that waveform polarity phenotypes only apply to triphasic individuals). For each of the three phenotype pairs, we calculated the mean and standard deviation for TMM values within each grouping, then extracted the genes that had (1) expression values more than four times greater in one phenotype than the other and (2) a difference in mean expression greater than either within-group standard deviation. This resulted in six lists of upregulated genes, one for each EOD feature across all OTUs and samples (Fig. 3, Set B).

In order to assess enrichment of particular gene pathways, biological functions, and cellular locations using a controlled vocabulary, we performed GO [86, 87] enrichment tests on every list of upregulated genes from (1) the ten pairwise comparisons (n = 20, two per comparison), (2) Set B (n = 6), and (3) Set C (n = 6), for each of the three ontology domains: Biological Process, Cellular Component, and Molecular Function. First, we identified homologous proteins predicted from the P. kingsleyae (N-type) reference genome and those predicted from Danio rerio (GRCz11) by blastp (BLAST+ v2.6, [88]). For each protein, the top hit (e-value ≤1e-10) was used for annotation. Next, we used mygene v1.14.0 [89, 90] to match the D. rerio proteins to D. rerio genes and extract their GO annotations (zebrafish Zv9). This resulted in GO annotations for each of the three ontology domains for P. kingsleyae (N-type) genes. Finally, we carried out the GO enrichment tests using topGO v2.30.1 [91] and the following parameters: nodeSize = 10, statistic = fisher, algorithm = weight01, p-value ≤0.02. The ‘universe’ for each enrichment test on gene lists from the pairwise comparisons was all the genes deemed expressed in the respective comparison, whereas the non-redundant list of genes in these ten ‘universes’ was the ‘universe’ for all enrichment tests on the gene lists from Sets B and C.

Interpretation of lists of genes from Set A and Set B each suffered limitations for the overall goals of this analysis, which is to identify the DEGs most strongly associated with each waveform feature (duration, complexity, and polarity). The ten comparisons made to construct Set A were not equally informative for two primary reasons: (1) the OTUs in this analysis vary in terms of their phylogenetic relatedness (see [30, 31]) and (2) several OTU comparisons varied in more than one waveform characteristic (Table 2). As such, we elected to focus on the most informative comparison for each EOD feature: the comparison that contrasted only the given feature and that minimized phylogenetic distance between OTUs. Of the ten pairwise comparisons, we classified three as the most informative comparisons, one per EOD feature (Table 2). The six lists of significantly upregulated genes (fold change > 4, FDR-corrected p-value < 0.001) from these three pairwise OTU comparisons constitute Set A’.

Comparisons in Set A’, however, lack biological replication. In contrast, interpretations of Set B were potentially limited in that many of the OTUs in this analysis differed in more than one EOD feature. To circumvent the limitations of Sets A’ and B within the limits of our study design, we constructed a third set (Set C). Set C is defined as the intersection of the upregulated genes from Sets A’ and B, for each phenotype. Since there were six phenotypes in our study, Set C encompasses six lists of upregulated genes and their respective enriched GO terms (Fig. 3). Therefore, Set C represents the DEGs that are (1) statistically supported as differentially expressed between closely-related OTUs that vary in a single waveform characteristic, and (2) are consistently differentially expressed among all OTUs that share that waveform feature. We focus our attention on Set C: We retrieved GO term definitions from QuickGO [92] and descriptions of gene function of the functional annotations from UniProt [93]; and to facilitate the discussion, we classified the more interesting genes in Set C into “general” functional classes, or themes.

Availability of data and materials

Raw sequence reads for all samples were deposited in the NCBI SRA with the BioProject Accession Number PRJNA573805. Per sample SRA accession numbers are listed in Additional file 1. All source code necessary to perform the methods described in this manuscript is provided in a GitHub repository: http://github.com/msuefishlab/paramormyrops_rnaseq

Abbreviations

DEG:

Differentially expressed gene

DGE:

Differential gene expression

ECM:

Extracellular matrix

EO:

Electric organ

EOD:

Electric organ discharge

GO:

Gene Ontology

NCBI:

National Center for Biotechnology Information

OTU:

Operational taxonomic unit

References

  1. 1.

    Butlin R, Debelle A, Kerth C, Snook RR, Beukeboom LW, Castillo Cajas RF, et al. What do we need to know about speciation? Trends Ecol Evol. 2012;27:27–39.

  2. 2.

    Givnish TJ. Adaptive radiation versus “radiation” and “explosive diversification”: why conceptual distinctions are fundamental to understanding evolution. New Phytol. 2015;207:297–303.

  3. 3.

    Grant RB, Grant PR. What Darwin’s finches can teach us about the evolutionary origin and regulation of biodiversity. Bioscience. 2003;53:965–75.

  4. 4.

    Kocher TD. Adaptive evolution and explosive speciation: the cichlid fish model. Nat Rev Genet. 2004;5:288–98.

  5. 5.

    Shaw KL. Further acoustic diversity in Hawaiian forests: two new species of Hawaiian cricket (Orthoptera: Gryllidae: Trigonidiinae: Laupala). Zool J Linnean Soc. 2000;129:73–91.

  6. 6.

    Mullen SP, Mendelson TC, Schal C, Shaw KL. Rapid evolution of cuticular hydrocarbons in a species radiation of acoustically diverse Hawaiian crickets (Gryllidae: Trigonidiinae: Laupala). Evolution (N Y). 2007;61:223–31.

  7. 7.

    Seehausen O. African cichlid fish: a model system in adaptive radiation research. Proc R Soc B Biol Sci. 2006;273:1987–98.

  8. 8.

    Hulsey CD. Cichlid genomics and phenotypic diversity in a comparative context. Integr Comp Biol. 2009;49:618–29.

  9. 9.

    Han F, Lamichhaney S, Grant BR, Grant PR, Andersson L, Webster MT. Gene flow, ancient polymorphism, and ecological adaptation shape the genomic landscape of divergence among Darwin’s finches. Genome Res. 2017;27:1004–15.

  10. 10.

    Lamichhaney S, Berglund J, Almén MS, Maqbool K, Grabherr M, Martinez-Barrio A, et al. Evolution of Darwin’s finches and their beaks revealed by genome sequencing. Nature. 2015;518:371–5.

  11. 11.

    Lamichhaney S, Han F, Berglund J, Wang C, Almén MS, Webster MT, et al. A beak size locus in Darwin’s finches facilitated character displacement during a drought. Science. 2016;352:470–4.

  12. 12.

    O’Quin CT, Drilea AC, Roberts RB, Kocher TD. A small number of genes underlie male pigmentation traits in Lake Malawi cichlid fishes. J Exp Zool Part B Mol Dev Evol. 2012;318:199–208.

  13. 13.

    Magalhaes IS, Seehausen O. Genetics of male nuptial colour divergence between sympatric sister species of a Lake Victoria cichlid fish. J Evol Biol. 2010;23:914–24.

  14. 14.

    Shaw KL, Lesnick SC. Genomic linkage of male song and female acoustic preference QTL underlying a rapid species radiation. Proc Natl Acad Sci. 2009;106:9737–42.

  15. 15.

    Wiley C, Ellison CK, Shaw KL. Widespread genetic linkage of mating signals and preferences in the Hawaiian cricket Laupala. Proc R Soc B Biol Sci. 2012;279:1203–9.

  16. 16.

    Blankers T, Oh KP, Bombarely A, Shaw KL. The genomic architecture of a Rapid Island radiation: recombination rate variation, chromosome structure, and genome assembly of the Hawaiian cricket Laupala. Genetics. 2018;209:1329–44.

  17. 17.

    Ding B, Daugherty DW, Husemann M, Chen M, Howe AE, Danley PD. Quantitative genetic analyses of male color pattern and female mate choice in a pair of cichlid fishes of Lake Malawi. East Africa PLoS One. 2014;9:1–22.

  18. 18.

    King MC, Wilson AC. Evolution at two levels in humans and chimpanzees. Science. 1975;188:107–16.

  19. 19.

    Alvarez M, Schrey AW, Richards CL. Ten years of transcriptomics in wild populations: what have we learned about their ecology and evolution? Mol Ecol. 2015;24:710–25.

  20. 20.

    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.

  21. 21.

    Ferreira PG, Patalano S, Chauhan R, Ffrench-Constant R, Gabaldón T, Guigó R, et al. Transcriptome analyses of primitively eusocial wasps reveal novel insights into the evolution of sociality and the origin of alternative phenotypes. Genome Biol. 2013;14:R20.

  22. 22.

    Abolins-Abols M, Kornobis E, Ribeca P, Wakamatsu K, Peterson MP, Ketterson ED, et al. Differential gene regulation underlies variation in melanic plumage coloration in the dark-eyed junco (Junco hyemalis). Mol Ecol. 2018;27:4501–15.

  23. 23.

    Härer A, Meyer A, Torres-Dowdall J. Convergent phenotypic evolution of the visual system via different molecular routes: How Neotropical cichlid fishes adapt to novel light environments. Evol Lett. 2018;2:341–54.

  24. 24.

    Nigenda-Morales SF, Hu Y, Beasley JC, Ruiz-Piña HA, Valenzuela-Galván D, Wayne RK. Transcriptomic analysis of skin pigmentation variation in the Virginia opossum (Didelphis virginiana). Mol Ecol. 2018;27:2680–97.

  25. 25.

    Young RL, Ferkin MH, Ockendon-Powell NF, Orr VN, Phelps SM, Pogány Á, et al. Conserved transcriptomic profiles underpin monogamy across vertebrates. Proc Natl Acad Sci. 2019;116:201813775.

  26. 26.

    Jeukens J, Renaut S, St-Cyr J, Nolte AW, Bernatchez L. The transcriptomics of sympatric dwarf and normal lake whitefish (Coregonus clupeaformis spp., Salmonidae) divergence as revealed by next-generation sequencing. Mol Ecol. 2010;19:5389–403.

  27. 27.

    Giorello FM, Feijoo M, D’Elía G, Naya DE, Valdez L, Opazo JC, et al. An association between differential expression and genetic divergence in the Patagonian olive mouse (Abrothrix olivacea). Mol Ecol. 2018;27:3274–86.

  28. 28.

    Carlson BA, Hasan SM, Hollmann M, Miller DB, Harmon LJ, Arnegard ME. Brain evolution triggers increased diversification of electric fishes. Science. 2011;332:583–6.

  29. 29.

    Rabosky DL, Santini F, Eastman J, Smith SA, Sidlauskas B, Chang J, et al. Rates of speciation and morphological evolution are correlated across the largest vertebrate radiation. Nat Commun. 2013;4:1–8.

  30. 30.

    Sullivan JP, Lavoué S, Hopkins CD. Discovery and phylogenetic analysis of a riverine species flock of African electric fishes (Mormyridae: Teleostei). Evolution (N Y). 2002;56:597–616.

  31. 31.

    Sullivan JP, Lavoué S, Arnegard ME, Hopkins CD. AFLPs resolve phylogeny and reveal mitochondrial introgression within a species flock of African electric fish (Mormyroidea: Teleostei). Evolution. 2004;58:825–41.

  32. 32.

    Lavoué S, Arnegard ME, Sullivan JP, Hopkins CD. Petrocephalus of Odzala offer insights into evolutionary patterns of signal diversification in the Mormyridae, a family of weakly electrogenic fishes from Africa. J Physiol Paris. 2008;102:322–39.

  33. 33.

    Hopkins CD. On the diversity of electric signals in a community of mormyrid electric fish in West Africa. Am Zool. 1981;21:211–22.

  34. 34.

    Bass AH. Species differences in electric organs of mormyrids: substrates for species-typical electric organ discharge waveforms. J Comp Neurol. 1986;244:313–30.

  35. 35.

    Arnegard ME, Bogdanowicz SM, Hopkins CD. Multiple cases of striking genetic similarity between alternate electric fish signal morphs in sympatry. Evolution. 2005;59:324–43.

  36. 36.

    Sullivan JP, Lavoué S, Hopkins CD. Molecular systematics of the African electric fishes (Mormyroidea: teleostei) and a model for the evolution of their electric organs. J Exp Biol. 2000;203(Pt 4):665–83.

  37. 37.

    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.

  38. 38.

    Lissmann HW, Machin KE. The mechanism of object location in Gymnarchus niloticus and similar fish. J Exp Biol. 1958;35:451–86.

  39. 39.

    von der Emde G, Amey M, Engelmann J, Fetz S, Folde C, Hollmann M, et al. Active electrolocation in Gnathonemus petersii: behaviour, sensory performance, and receptor systems. J Physiol - Paris. 2008;102:279–90.

  40. 40.

    Kramer B. Electric organ discharge interaction during interspecific agonistic behaviour in freely swimming mormyrid fish. J Comp Physiol. 1974;93:203–35.

  41. 41.

    Möhres FP. Elektrische Entladungen im Dienste der Revierabgrenzung bei Fischen. Naturwissenschaften. 1957;44:431–2.

  42. 42.

    Hopkins CD. Signal evolution in electric communication. In: Hauser M, Konishi M, editors. The design of animal communication. Cambridge, MA: MIT Press; 1999. p. 461–91.

  43. 43.

    Gallant JR, Arnegard ME, Sullivan JP, Carlson BA, Hopkins CD. Signal variation and its morphological correlates in Paramormyrops kingsleyae provide insight into the evolution of electrogenic signal diversity in mormyrid electric fish. J Comp Physiol A. 2011;197:799–817.

  44. 44.

    Carlson BA. Electric signaling behavior and the mechanisms of electric organ discharge production in mormyrid fish. J Physiol. 2002;96:405–19.

  45. 45.

    Caputi AA, Carlson BA, Macadar O. Electric Organs and Their Control. In: Bullock TH, Hopkins C.D., Popper A.N., Fay R.R., editors. Electroreception. New York: Springer Handbook of Auditory Research, vol 21; 2005. p. 410–51.

  46. 46.

    Hopkins CD. Behavior of Mormyridae. In: Bullock TH, Heiligenberg W, editors. Electroreception. New York: Wiley; 1986. p. 527–76.

  47. 47.

    Bennett MVL, Grundfest H. Studies on the morphology and electrophysiology of electric organs. III. Electrophysiology of electric organs in mormyrids. In: Chagas C, de Carvalho A, editors. Bioelectrogenesis. Amsterdam: Elsevier; 1961. p. 113–135.

  48. 48.

    Bennett MVL. Electric Organs. In: Hoar WS, Randall DJ, editors. Fish physiology. London: Academic Press; 1971. p. 347–491.

  49. 49.

    Bass AH, Volman SF. From behavior to membranes: testosterone-induced changes in action potential duration in electric organs. Proc Natl Acad Sci. 1987;84:9295–8.

  50. 50.

    Bass AH, Denizot J-P, Marchaterre MA. Ultrastructural features and hormone-dependent sex differences of mormyrid electric organs. J Comp Neurol. 1986;254:511–28.

  51. 51.

    Szabo T. Les Organes Electriques des Mormyrides. In: Chagas C, de Carvalho A, editors. Bioelectrogenesis. New York: Elsevier; 1961. p. 20–4.

  52. 52.

    Carlson BA, Gallant JR. From sequence to spike to spark: evo-devo-neuroethology of electric communication in mormyrid fishes. J Neurogenet. 2013;27:106–29.

  53. 53.

    Zakon HH, Lu Y, Zwickl D, 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.

  54. 54.

    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.

  55. 55.

    Nagel R, Kirschbaum F, Tiedemann R. Electric organ discharge diversification in mormyrid weakly electric fish is associated with differential expression of voltage-gated ion channel genes. J Comp Physiol A Neuroethol Sensory, Neural, Behav Physiol. 2017;203:183–95.

  56. 56.

    Swapna I, Ghezzi A, York JM, Markham MR, Halling DB, Lu Y, et al. Electrostatic Tuning of a Potassium Channel in Electric Fish. Curr Biol. 2018;28:2094–2102.e5.

  57. 57.

    Paul C, Kirschbaum F, Mamonekene V, Tiedemann R. Evidence for non-neutral evolution in a Sodium Channel gene in African weakly electric fish (Campylomormyrus, Mormyridae). J Mol Evol. 2016;83:61–77.

  58. 58.

    Gallant JR, Losilla M, Tomlinson C, Warren WC. The genome and adult somatic Transcriptome of the Mormyrid electric fish Paramormyrops kingsleyae. Genome Biol Evol. 2017;9:3525–30.

  59. 59.

    Gallant JR, Sperling J, Cheng C, Arnegard ME, Carlson BA, Hopkins CD. Microevolutionary processes underlying macroevolutionary patterns of electric signal diversity in mormyrid fish. bioRxiv.org. 2017;:154047.

  60. 60.

    Wilwert JL. Madhoun, Nisreen M, Coughlin DJ. Parvalbumin correlates with relaxation rate in the swimming muscle of sheepshead and kingfish. J Exp Biol. 2006;209:227–37.

  61. 61.

    Schwartz IR, Pappas GD, Bennett MVL. The fine structure of electrocytes in weakly electric teleosts. J Neurocytol. 1975;4:87–114.

  62. 62.

    Bass AH, Hopkins CD. Hormonal control of sexual differentiation: changes in electric organ discharge waveform. Science. 1983;220:971–4.

  63. 63.

    Bass AH, Hopkins CD. Hormonal control of sex differences in the electric organ discharge (EOD) of mormyrid fishes. J Comp Physiol A. 1985;156:587–604.

  64. 64.

    Freedman EG, Olyarchuk J, Marchaterre MA, Bass AH. A temporal analysis of testosterone-induced changes in electric organs and electric organ discharges of mormyrid fishes. J Neurobiol. 1989;20:619–34.

  65. 65.

    Lamanna F, Kirschbaum F, Waurick I, Dieterich C, Tiedemann R. Cross-tissue and cross-species analysis of gene expression in skeletal muscle and electric organ of African weakly-electric fish (Teleostei; Mormyridae). BMC Genomics. 2015;16:668.

  66. 66.

    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.

  67. 67.

    Szabo T. Development of the electric organ of Mormyridae. Nature. 1960;188:760–2.

  68. 68.

    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.

  69. 69.

    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:324–50.

  70. 70.

    Hopkins CD. Design features for electric communication. J Exp Biol. 1999;202:1217–28.

  71. 71.

    Ganss B, Hoffmann W. Calcium-induced conformational transition of trout ependymins monitored by tryptophan fluorescence. Open Biochem J. 2009;3:14–7.

  72. 72.

    Suárez-Castillo EC, García-Arrarás JE. Molecular evolution of the ependymin protein family: a necessary update. BMC Evol Biol. 2007;7:23.

  73. 73.

    Bartels E. Depolarization of electroplax membrane in calcium-free ringer’s solution. J Membr Biol. 1971;5:121–32.

  74. 74.

    Ferrari MB, Zakon HH. Conductances contributing to the action-potential of Sternopygus electrocytes. J Comp Physiol A-sensory neural Behav Physiol. 1993;173:281–92.

  75. 75.

    Gallant JR, Traeger LL, Volkening JD, Moffett H, Chen P-H, Novina CD, et al. Genomic basis for the convergent evolution of electric organs. Science. 2014;344:1522–5.

  76. 76.

    Lamanna F, Kirschbaum F, Tiedemann R. De novo assembly and characterization of the skeletal muscle and electric organ transcriptomes of the African weakly electric fish Campylomormyrus compressirostris (Mormyridae, Teleostei). Mol Ecol Resour. 2014;14:1222–30.

  77. 77.

    Fraser JRE, Laurent TC, Laurent UBG. Hyaluronan: its nature, distribution, functions and turnover. J Intern Med. 1997;242:27–33.

  78. 78.

    Arnegard ME, Jackson BS, Hopkins CD. Time-domain signal divergence and discrimination without receptor modification in sympatric morphs of electric fishes. J Exp Biol. 2006;209(Pt 11):2182–98.

  79. 79.

    Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics. 2014;30:2114–20.

  80. 80.

    MacManes MD. On the optimal trimming of high-throughput mRNA sequence data. Front Genet. 2014;5:1–7.

  81. 81.

    Langmead B, Salzberg SL. Fast gapped-read alignment with bowtie 2. Nat Methods. 2012;9:357–9.

  82. 82.

    Li B, Dewey CN. RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinformatics. 2011;12:323.

  83. 83.

    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.

  84. 84.

    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.

  85. 85.

    Robinson MD, Oshlack A. A scaling normalization method for differential expression analysis of RNA-seq data. Genome Biol. 2010;11:R25.

  86. 86.

    Ashburner M, Ball CA, Blake JA, Botstein D, Butler H, Cherry JM, et al. Gene ontology: tool for the unification of biology. Nat Genet. 2000;25:25–9.

  87. 87.

    The Gene Ontology Consortium. The gene ontology resource: 20 years and still GOing strong. Nucleic Acids Res. 2019;47:D330–8.

  88. 88.

    Camacho C, Coulouris G, Avagyan V, Ma N, Papadopoulos J, Bealer K, et al. BLAST+: architecture and applications. BMC Bioinformatics. 2009;10:421.

  89. 89.

    Wu C, MacLeod I, Su AI. BioGPS and MyGene.info: organizing online, gene-centric information. Nucleic Acids Res. 2013;41:D561–5.

  90. 90.

    Xin J, Mark A, Afrasiabi C, Tsueng G, Juchler M, Gopal N, et al. High-performance web services for querying gene and variant annotation. Genome Biol. 2016;17:91.

  91. 91.

    Alexa A, Rahnenführer J, Lengauer T. Improved scoring of functional groups from gene expression data by decorrelating GO graph structure. Bioinformatics. 2006;22:1600–7.

  92. 92.

    Binns D, Dimmer E, Huntley R, Barrell D, O’Donovan C, Apweiler R. QuickGO: a web-based tool for gene ontology searching. Bioinformatics. 2009;25:3045–6.

  93. 93.

    The UniProt Consortium. UniProt: a worldwide hub of protein knowledge. Nucleic Acids Res. 2019;47:D506–15.

Download references

Acknowledgements

The authors wish to thank Bruce Carlson, Matt Arnegard, Carl Hopkins, Roger Afene, Jean Danielle Mbega and Marie-Francois Eva for assistance with specimen collection in Gabon. The authors also acknowledge the Michigan State University Institute for Cyber-Enable Research for use of their high-performance computing infastructure, and CENREST Gabon for logistical support. The authors also thank two anonymous reviewers for their helpful comments on an earlier version of this manuscript.

Funding

This work was supported by the National Science Foundation (1455405: JRG, PI) and a grant from the Cornell University Center for Vertebrate Genomics. The funders had no role in the design of the study, in the collection, analysis and interpretation of data, or in writing the manuscript.

Author information

JRG designed the experiment, collected and identified specimens, performed library construction and sequencing, and oversaw data analysis. ML performed quality control analysis, designed and implemented the data analysis procedure, and performed examination of gene ontology and gene function. DL designed and implemented portions of the data analysis procedures. All authors contributed to writing and editing the manuscript, and have read and approved its final version.

Correspondence to Jason R. Gallant.

Ethics declarations

Ethics approval and consent to participate

All research protocols involving live fish were approved by the Michigan State University and Cornell University Institutional Animal Care and Use Committees. The Centre National de la Reserche Scientifique (CENAREST) of Gabon provided permissions for collection and export of specimens.

Consent for publication

Not applicable.

Competing interests

The authors declare they have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary information

Additional file 1 Raw reads NCBI SRA accession numbers, number of reads and alignment rates per sample, using bowtie 2 as the aligner and the Paramormyrops kingsleyae (N-type) transcriptome as the reference.

Additional file 2. DEG per comparison from the 10 pairwise DGE analysis. Positive values under logFC indicate genes upregulated in the OTU under sampleA, whereas negative values correspond to genes upregulated in the OTU under sampleB. Values under each sample are gene raw counts. Significance threshold was abs (log (base2)FC) > 2 (= 4-fold expression difference) and FDR < 0.001.

Additional file 3. Enriched GO terms per comparison, ontology and OTU in the DEG from the 10 pairwise comparisons. Also listed are the DEG annotated to each GO term. The pvalue is in the column weight01.

Additional file 4. DEG per EOD feature and phenotype identified with the Set B analysis. Values under each sample are TMM normalized, log2(TMM + 1) transformed, and mean-centered expression values.

Additional file 5. Enriched GO terms per EOD feature, ontology and phenotype in the DEG from the Set B analysis. Also listed are the DEG annotated to each GO term. The pvalue is in the column weight01.

Additional file 6. DEG in Set C, per EOD feature and phenotype.

Additional file 7. GO terms enriched in the DEG in Set C, per EOD feature, ontology and phenotype. Also listed are the DEG annotated to each GO term, and the quickGO definitions of each GO term. The pvalue is in the column weight01.

Additional file 8. MA plots from the 10 pairwise DGE analysis. Red dots represent genes with FDR < 0.05 (Trinity’s default parameters).

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Losilla, M., Luecke, D.M. & Gallant, J.R. The transcriptional correlates of divergent electric organ discharges in Paramormyrops electric fish. BMC Evol Biol 20, 6 (2020) doi:10.1186/s12862-019-1572-3

Download citation