Gene and Allele-Specific Expression Underlying the Electric Signal Divergence in African Weakly Electric Fish

Abstract

In the African weakly electric fish genus Campylomormyrus, electric organ discharge signals are strikingly different in shape and duration among closely related species, contribute to prezygotic isolation, and may have triggered an adaptive radiation. We performed mRNA sequencing on electric organs and skeletal muscles (from which the electric organs derive) from 3 species with short (0.4 ms), medium (5 ms), and long (40 ms) electric organ discharges and 2 different cross-species hybrids. We identified 1,444 upregulated genes in electric organ shared by all 5 species/hybrid cohorts, rendering them candidate genes for electric organ–specific properties in Campylomormyrus. We further identified several candidate genes, including KCNJ2 and KLF5, and their upregulation may contribute to increased electric organ discharge duration. Hybrids between a short (Campylomormyrus compressirostris) and a long (Campylomormyrus rhynchophorus) discharging species exhibit electric organ discharges of intermediate duration and showed imbalanced expression of KCNJ2 alleles, pointing toward a cis-regulatory difference at this locus, relative to electric organ discharge duration. KLF5 is a transcription factor potentially balancing potassium channel gene expression, a crucial process for the formation of an electric organ discharge. Unraveling the genetic basis of the species-specific modulation of the electric organ discharge in Campylomormyrus is crucial for understanding the adaptive radiation of this emerging model taxon of ecological (perhaps even sympatric) speciation.

African weakly electric fish, electric organ discharge, adaptive radiation, potassium channel, gene expression, allele-specific expression

Introduction

Electric fish have independently evolved 6 times (Darwin 1859; Bass 1986; Kirschbaum and Formicki 2020). They possess a specific myogenic electric organ (EO) derived from skeletal muscle (SM) fibers except for Apteronotidae, which possess an EO derived from nervous tissue (Smith 2013). Comparative genomics have unraveled this convergent phenotypic evolution to originate in part also from convergence on the molecular level: both voltage-dependent sodium and potassium channels are involved in the EO development and physiology. Because of the teleost-specific whole genome duplication (Glasauer and Neuhauss 2014), these fish possess 2 copies of most genes and subfunctionalization among paralogs and differential expression between EO and SM seem to play a major role in the transition of myocytes to electrocytes. A prominent example is the voltage-gated sodium (Na_v) channel gene (SCN4a): convergently in 3 electrogenic taxa (Mormyridae, Siluriformes, and Gymnotiformes), only 1 paralog (SCN4ab) is still expressed in SM, but the other 1 (SCN4aa) is exclusively expressed in the EO, indicating a crucial role for electrogenesis (Zakon 2012; Wang and Yang 2021; LaPotin et al. 2022). The Na_v channel gene (SCN4aa) is regulated by FGF13a in the 3 electric fish lineages Siluriformes, Gymnotiformes, and Mormyridae (Gallant et al. 2014). Differential expression of multiple isoforms of α- and β-subunits of sodium/potassium ATPase is important in the EO as well (Gallant et al. 2012, 2014; Lamanna et al. 2015). In addition, several transcription factors, HEY1, MEF2a, and SIX2a, are convergently upregulated in the EOs of those electric fish lineages (Kim et al. 2004; Gallant et al. 2012, 2014). EOs hence comprise a prime example of convergent evolution in both genotype and phenotype.

One of the electric fish clades, mormyrid fish, contains about 200 described species that are endemic to Africa. This outstanding adaptive radiation within the otherwise species-poor basal lineage of osteoglossiforms is putatively due to their species-specific weak electric signals, which are used for both electrolocation and electrocommunication (Feulner et al. 2009a, b). Divergence in EO discharge (EOD) is considered a major driver in the ecological (and possibly sympatric) speciation in the mormyrid genus Campylomormyrus, which is mainly distributed in the Congo River (Tiedemann et al. 2010).

The genus Campylomormyrus comprises 15 described species, which have profoundly diverged in their EOD with regard to signal duration and waveform (Feulner et al. 2009a). Those species possess either long or short, biphasic or triphasic, but always species-specific EODs that function as a prezygotic reproductive isolation mechanism and are supposed to have arisen via divergent selection among closely related species (Feulner et al. 2009a). In adult Campylomormyrus, the EO, confined to the caudal peduncle (Fig. 1a), is composed of specialized electrocytes (Paul et al. 2015). They have a flat, disk-shaped appearance with a clear orientation toward the longitudinal body axis (Fig. 1b). Unlike SM myocytes, electrocytes possess a number of special evaginations, called stalks, mostly on the posterior face (Paul et al. 2015). These stalks are either fused into major stalks on the posterior face (Fig. 1c, left) or they penetrate the electrocyte and merge at the anterior face to constitute to major stalks (Fig. 1c, right). A branch of the spinal nerve forms numerous synapses with the major stalk, whether on the posterior or on the anterior face of the electrocyte, and the action potentials are propagated along the stalk system to the disk-like part of the electrocyte (Paul et al. 2015). The externally measurable EOD is formed by simultaneous action potentials of all electrocytes. The shape of the EOD in Campylomormyrus is often associated with the penetration of the stalks (Gallant et al. 2011), while the structural basis of the EOD duration, which can vary 100-fold across species, is still only partially understood. A very elongated EOD (∼40 ms) is produced by Campylomormyrus rhynchophorus and Campylomormyrus numenius, which exhibit large foldings or evaginations on the anterior face of the electrocytes, so-called papillae (Kirschbaum et al. 2016; Korniienko et al. 2021). In 2 species with relatively short EOD (Fig. 2a), Campylomormyrus compressirostris (0.4 ms) and Campylomormyrus tamandua (0.4 ms), many small stalks fuse into 1 major stalk of large diameter after their origin (Paul et al. 2015). In contrast, the stalk system in species with an EOD of medium (e.g. Campylomormyrus tshokwe, 5 ms) or long duration (e.g. C. numenius, 40 ms) is more branched (Paul et al. 2015). Apart from these differences in the stalk system, species with highly diverged EOD waveforms still show similar electrocyte geometry suggesting further core mechanisms to contribute to the observed EOD variations. Since the electrocytes generate action potentials for EOD, the distribution and repertoire of ion currents have long been considered to play a key role in EOD formation (Lamanna et al. 2014, 2015; Paul et al. 2016; Nagel et al. 2017; Swapna et al. 2018).

Fig. 1.

EO and electrocyte structure in Campylomormyrus and schematic illustration for potassium channels. a) Electric organ (EO) and electric organ discharge (EOD) in an adult Campylomormyrus fish. b) The EO consists of four columns of electrocytes (e) which surround the vertebral column (vc), the stalk system (st) is connected to the posterior face of the electrocyte. c) Anterior (A.) and posterior (P.) faces of electrocytes with two types of stalk system. Panel c is modified from Gallant et al. (2012). d) Schematic illustration of voltage-gated potassium (K_v) channel and inwardly rectifying (K_ir) channel subunit. K_v channel subunit contains six transmembrane (TM) helices, a pore-forming (H5) loop, and cytosolic NH₂ (N) and COOH (C) termini. The gene KCN7A_2 was inferred to be under positive selection and the mutation encodes the loop between TM3-4. K_ir channel subunit contains only two TMs.

Fig. 2.

EOD shape and duration of Campylomormyrus species and hybrids and the working flow of this study. a) Species/hybrids samples used in the study and their EOD pattern. b) Differential gene expression (DGE) analysis between electric organ (EO) and skeletal muscle (SM) for each species/hybrid to identify genes with EO-specific expression. c) RNA-seq data clustering based on EOD duration change for EO (red) and SM (blue) in F0 species. d) Allele specific expression analysis in each hybrid set.

Sodium and potassium fluxes are considered the most important ion currents in controlling the EOD (Stoddard and Markham 2008). They are the basic requirements for generating an action potential (Mehaffey et al. 2006). Consequently, abundance and properties of sodium and potassium channels are likely to profoundly influence the EOD. The potassium channels can be classified into different classes based on their structure and function: voltage-gated (K_v, includes subfamilies, e.g. shaker-related KCNA and shab-related KCNB), inwardly rectifying (K_ir), tandem pore domain channels (K_2p), ligand-gated channels, and calcium-activated channels (K_ca) (Kuang et al. 2015). Two paralogs of the KCNA7 channel gene originate from the whole genome duplication event in teleost fish, and these paralogs might have undergone subfunctionalization or neofunctionalization in mormyrids: one of them KCNA7a is predominantly expressed in the EO of mormyrids, while KCNA7b is preferentially expressed in SM (Swapna et al. 2018). The K_v channel contains 6 transmembrane helices (Fig. 1d). In KCNA7a, a nonsynonymous substitution was observed in the transmembrane helix 3 and 4 linkers and the encoded amino acid substitution might relate to the EOD duration difference among the mormyrid taxa Brienomyrus and Gymnarchus (Swapna et al. 2018).

This study focuses on potential molecular mechanisms underlying the divergent EOD among Campylomormyrus species as a potential major driver of their adaptive radiation. This study takes further advantage of artificially bred hybrid electric fish. Campylomormyrus species hybrids often exhibit an adult EOD, which is similar to the juvenile EOD from one of the parental species, and the adult EOD duration in hybrids is usually intermediate between the 2 parental species (Kirschbaum et al. 2016). Gene expression analyses in hybrids further enable assessment of allelic-specific expression, relative to the expressed trait of interest (here, EOD duration). To enhance our understanding of the genetic regulation of EOD divergence among Campylomormyrus species, especially for the EOD duration divergence, we (i) compared the gene expression pattern between EO and SM in the 3 F0 species, C. compressirostris (com, short and biphasic EOD), C. tshokwe (tsh, medium and biphasic EOD), and C. rhynchophorus (rhy, long and triphasic EOD), and 2 F1 hybrids, C. compressirostris ♂ × C. tshokwe ♀ (com × tsh, short and biphasic EOD) and C. compressirostris ♂ × C. rhynchophorus ♀ (com × rhy, medium and biphasic EOD); (ii) clustered RNA-seq data relative to the EOD duration in 3 F0 species to infer genes with duration-specific expression; and (iii) assessed biallelic-specific expression for 2 hybrid sets (each set includes 2 F0 parental species and their hybrid; Fig. 2).

Results

We examined overall patterns in gene expression using a principal component analysis (PCA) based on all expressed genes (Fig. 3a). Expression profiles of SMs and EOs were broadly separated along PC1, which explained 74% of variance. The SMs’ expression profiles from all species/hybrids clustered together; however, species/hybrid-specific EOs’ expression profiles were stratified along PC2 (explained 6% of variance), relative to EOD duration (Fig. 3a). The PCA hence indicates that gene expression in Campylomormyrus (i) is EO specific, compared with SM, and (ii) relates to EOD duration, enabling inference of underlying candidate genes.

Fig. 3.

Differential gene expression between EO and SM in C. compressirostris (com), C. rhynchophorus (rhy), C. tshokwe (tsh), and hybrids C. compressirostris ♂ × C. rhynchophorus ♀ (com × rhy) and C. compressirostris ♂ × C. tshokwe ♀ (com × tsh). a) Principal component analysis (PCA) of gene expression levels between EO and SM in 5 species/hybrids. b) Venn Diagram graph for up (left) and down (right) regulated genes shared in 5 species/hybrids. All differentially expressed genes (DEGs) have |log₂(fold change)| >1 and a p-value < 0.05. Many of the DEGs are related to “membrane” and “plasma membrane” (see supplementary fig. S1). c) Volcano plot showing genes differentially expressed in EO (relative to SM) in all 5 species/hybrids. X-axis is the average log₂(fold change) among 5 species/hybrids, and y-axis is the associated –log₁₀ (average p-value for 5 species/hybrids). Potential candidate genes and genes with low p-value or high fold change are labeled with their name.

Genes with EO-Specific Expression Pattern

Differential gene expression analysis was used for pairwise comparisons between EO and SM for each species and hybrid (Fig. 2b). We identified significantly differentially expressed genes (DEGs) based on a |log₂ folder change (log₂FC)| > 1 and a P < 0.05. We specifically identified genes with an EO-specific expression pattern shared among all Campylomormyrus species/hybrids. There were 1,444 upregulated and 1,262 downregulated DEGs that were shared in the comparison of EO and SM in all species/hybrids (Fig. 3b and c).

Among the DEGs upregulated in the EO, 54 genes were related to transmembrane ion transport (Fig. 3c and Table 1; supplementary table S1, Supplementary Material online). We identified 4 genes encoding sodium/potassium-ATPase α- and β-subunit (ATP1a1, ATP1a2a, ATP1b1a, and ATP1b1b) and 3 Na_v channel genes (SCN4aa, SCN4b, and SCN1ba). Several genes encoding for different types of potassium channels were also identified: 4 K_v channel genes (KCNA7a_1, KCNA7a_2, KCNIP3, and KCNQ5), 2 K_ir channel genes (KCNJ2 and KCNJ9), and 1 K_2p channel gene (KCNK2). Further transmembrane ion transport DEGs were chloride, calcium, and other cation channel genes (supplementary table S1, Supplementary Material online). Several solute carrier family genes were also upregulated in the EO, in particular SLC24a2 (Table 1).

Table 1

Candidate genes in EO and their descriptions, including the average log₂FC and average P-value

ID	Blast gene	Highlights of predicted function	Gene description	Category	Average log2FC	Average P-value
maker-ptg000361l-augustus-gene-0.2-mRNA-1	ACTR3b	F-actin dynamics/polymerization	ARP3 actin-related protein 3 homolog B	Cytoskeletal and sarcomeric	3.02	3.16E−20
maker-ptg000346l-snap-gene-25.173-mRNA-1	GSN	F-actin dynamics/polymerization	Gelsolin	Cytoskeletal and sarcomeric	6.70	2.7189E–50
snap_masked-ptg000028l-processed-gene-137.57-mRNA-1	MYO15a	Unconventional myosin; actin-based motor protein	Unconventional myosin-XV	Cytoskeletal and sarcomeric	3.68	3.4002E–09
maker-ptg001003l-snap-gene-2.16-mRNA-1	MYO1e	Unconventional myosin; actin-based motor protein	Unconventional myosin-Ie	Cytoskeletal and sarcomeric	2.63	5.7074E–05
maker-ptg002090l-augustus-gene-35.16-mRNA-1	MYO3b	Unconventional myosin; actin-based motor protein	Myosin-IIIb	Cytoskeletal and sarcomeric	3.75	1.4271E–07
maker-ptg000215l-snap-gene-13.31-mRNA-1	S100b	Cytosolic Ca2+-binding protein of the EF-hand superfamily	S100 calcium binding protein B	Signaling	8.14	1.0122E–22
maker-ptg000049l-snap-gene-23.20-mRNA-1	FGF12	Possibly regulate voltage-gated sodium channels	Fibroblast growth factor 12	Signaling	8.72	5.8498E–17
maker-ptg000838l-snap-gene-3.218-mRNA-1	NDRG3	Predicted to be involved in signal transduction	N-Myc downstream-regulated gene 3 protein	Signaling	11.02	7.4446E–08
maker-ptg001563l-augustus-gene-5.35-mRNA-1	SGK1	Serine/threonine-protein kinase	Serine/threonine-protein kinase Sgk1	Signaling	7.79	6.6327E–38
maker-ptg001088l-snap-gene-0.7-mRNA-1	SIX2a	Target ARE promoter elements in sodium/potassium adenosine triphosphatases	SIX homeobox 2	Transcription factor	3.05	1.221E–27
maker-ptg000783l-snap-gene-6.78-mRNA-1	HEY1	Developing cardiac conduction pathway	Hes-related family bHLH transcription factor with YRPW motif 1	Transcription factor	6.02	1.5325E–13
maker-ptg001740l-augustus-gene-1.112-mRNA-1	KLF5	Rebalance potassium channels	Krüppel-like factor 5	Transcription factor	8.39	5.0491E–06
maker-ptg000008l-snap-gene-10.43-mRNA-1	MEF2a	Transcriptional activator for numerous muscle-specific genes	Myocyte-specific enhancer factor 2A	Transcription factor	3.85	3.7818E–50
maker-ptg001270l-snap-gene-47.16-mRNA-1	MEF2b	Transcriptional activator for numerous muscle-specific genes	Myocyte-specific enhancer factor 2B	Transcription factor	6.92	8.4583E–68
maker-ptg000970l-augustus-gene-2.127-mRNA-1	ATP1a1	Sodium/potassium-ATPase α-subunit	Sodium/potassium-transporting ATPase subunit alpha-1	Transmembrane ion transport	10.55	2.1894E–05
snap_masked-ptg001156l-processed-gene-0.19-mRNA-1	ATP1a2a	Sodium/potassium-ATPase α-subunit	Sodium/potassium-transporting ATPase subunit alpha-2	Transmembrane ion transport	3.07	2.6317E–12
maker-ptg001047l-snap-gene-1.63-mRNA-1	ATP1b1a	Sodium/potassium-ATPase β-subunit	ATPase sodium/potassium transporting beta 1a	Transmembrane ion transport	6.59	3.6325E–59
maker-ptg000509l-snap-gene-9.39-mRNA-1	ATP1b1b	Sodium/potassium-ATPase β-subunit	ATPase sodium/potassium transporting beta 1b	Transmembrane ion transport	4.51	2.0647E–06
maker-ptg000028l-snap-gene-81.10-mRNA-1	KCNA7a_1	Kv channel	Potassium voltage-gated channel subfamily A member 7a	Transmembrane ion transport	4.60	2.5847E–11
maker-ptg000028l-snap-gene-81.8-mRNA-1	KCNA7a_2	Kv channel	Potassium voltage-gated channel subfamily A member 7a	Transmembrane ion transport	8.27	3.5553E–12
maker-ptg001427l-snap-gene-13.20-mRNA-1	KCNIP3	Kv channel	Calsenilin	Transmembrane ion transport	5.11	0.00099357
maker-ptg000697l-snap-gene-6.109-mRNA-1	KCNQ5	Kv channel	Potassium voltage-gated channel subfamily Q member 5	Transmembrane ion transport	5.94	0.00093894
maker-ptg000265l-est_gff_est2genome-gene-6.33-mRNA-1	KCNJ2	Kir channel	Inward rectifier potassium channel 2	Transmembrane ion transport	5.53	1.1222E–20
maker-ptg000830l-augustus-gene-5.123-mRNA-1	KCNJ9	Kir channel	G protein-activated inward rectifier potassium channel 3	Transmembrane ion transport	5.77	2.6324E–10
snap_masked-ptg001118l-processed-gene-0.13-mRNA-1	KCNK2	K2p channel	Potassium channel subfamily K member 2	Transmembrane ion transport	6.15	9.238E–14
maker-ptg000253l-augustus-gene-20.10-mRNA-1	SCN1ba	Nav channel	Sodium channel subunit beta-1	Transmembrane ion transport	3.84	1.2983E–07
maker-ptg001188l-snap-gene-6.4-mRNA-1	SCN4aa	Nav channel	Sodium channel protein type 4 subunit alpha A	Transmembrane ion transport	10.98	1.2737E–11
maker-ptg002239l-snap-gene-5.5-mRNA-1	SCN4b	Nav channel	Sodium voltage-gated channel beta subunit 4	Transmembrane ion transport	5.34	8.1806E–47
maker-ptg000148l-snap-gene-9.4-mRNA-1	SLC24a2	Calcium, potassium:sodium antiporter	Solute carrier family 24 member 2	Transmembrane ion transport	9.06	1.1057E–37

Eighteen genes upregulated in the EO were associated with cytoskeletal and sarcomeric protein (supplementary table S1, Supplementary Material online). The predicted functions of those genes were mainly related to F-actin dynamics and unconventional myosin activity (Table 1). A signaling gene NDRG3 showed very high overexpression in EO (log₂FC = 11.02), as well as the genes SGK2, S100b, and FGF12. The upregulated transcription factors in the EO included Krüppel-like factor 5 (KLF5), FOXL2, SIX2a, HEY1, and 2 myocyte-specific enhancer factors (MEF2a and MEF2b).

In the downregulated DEGs in EO (or upregulated in SM), 44 genes were classified into the category “cytoskeletal and sarcomeric” (Fig. 3c; supplementary table S2, Supplementary Material online). There were 37 transmembrane ion transport genes downregulated in EO, which were related to the ions potassium, sodium, and calcium. In contrast to the expression pattern of the 2 KCNA7a copies, 5 K_v1 subfamily genes (KCNA1b, KCNA4a, KCNA5b, KCNA6a, and KCNA7b) were downregulated in the EO. This was also the case for other potassium and sodium channel genes, e.g. K_v subfamily genes (KCNB1, KCNE4, and KCNIP4), K_2p subfamily genes (KCNK4, KCNK7), a K_ca subfamily gene (KCNN4), and Na_v channel genes (SCN3b, SCN4ab). Two muscle-specific transcription factors, MYOCD and MYOG, were also downregulated in EO (supplementary table S2, Supplementary Material online).

We applied a Gene Ontology (GO) enrichment analysis to further examine the function of all the up- and downregulated DEGs in EO, respectively (Dennis et al. 2003). Among the upregulated DEGs in the EO, there were 44 significantly enriched GO terms (Fisher's exact test P < 0.01; supplementary fig. S1 and table S3, Supplementary Material online). Among them, the 3 GO terms with the highest number of DEGs were all related to the cell membrane: membrane (464 DEGs), integral component of membrane (309 DEGs), and plasma membrane (237 DEGs). There were 47 DEGs assigned to the enriched GO term “ion transport.” Sixty-two and 35 DEGs were assigned to the enriched Golgi-related GO terms “Golgi membrane” and “Golgi apparatus,” respectively. In addition, there were 23 DEGs assigned to the enriched GO term “actin filament binding.” There were 73 GO terms significantly enriched for DEGs downregulated in the EO (upregulated in SM; supplementary fig. S2 and table S4, Supplementary Material online). They were associated with skeletal and cardiac muscle tissue–related GO terms.

Genes with Expression Levels Related to EOD Duration

The PCA plot from transcriptome-wide gene expression showed a significant association between overall gene expression and EOD duration in all species/hybrids (PC2 in Fig. 3a; accounting for 6% of the variance in gene expression). DESeq2 provides a likelihood ratio test (LRT) that compares how well a gene's read count data fit a “full model” (with independent variables) compared to a “reduced model” (without those variables). Therefore, it is well suited to explore whether there are any significant associations of gene expression levels across a series of values of an independent variable (here, EOD duration; Love et al. 2015). Specifically, we used this approach to test whether a gene's expression fits a pattern of increasing or decreasing over the different durations in 2 different tissues, EO and SM (Bendjilali et al. 2017). In order to avoid any bias potentially stemming from distorted expression pattern in the hybrids, we only used the quantification data from the parental purebred (F0) species. The LRT analysis returned 1,874 significant genes using a threshold of adjusted P-value (P_adj) < 0.05. Those genes were further sorted into groups using the degPatterns function (Pantano 2019). Each such group contained genes following a specific pattern of expression across the different duration values in the analyzed tissues EO and SM (Pantano 2019).

The degPatterns function generated 27 groups of different expression patterns in EO and SM, relative to EOD duration (supplementary fig. S3, Supplementary Material online). To identify EOD duration–specific genes, we focused on the groups meeting the following criteria: (i) the gene expression level in EO is higher than SM in all F0 species (i.e. the gene is consistently upregulated in the EO), and (ii) the gene expression level in the EO shows an increasing or decreasing pattern, relative to EOD duration. Two groups showed a consistent increasing expression pattern (group 5, 41 genes; group 6, 239 genes) and 1 a decreasing expression pattern (group 3, 405 genes), relative to the EOD duration (Fig. 4).

Fig. 4.

RNA-seq data clustered by EOD duration (only for the 3 purebred species).

In the increased expression pattern groups (5 and 6), we found K_ir subfamily gene KCNJ2 and the transcription factor KLF5; both were found among the genes with EO-specific expression as well (Table 2). In the decreased expression group 3, there were 2 transmembrane ion transport genes (KCNK6 and KCNQ5) and 2 cytoskeletal and sarcomeric genes (ACTR3b and NHS; Table 2).

Table 2

Genes with expression correlated to EOD duration and their description

Group	Gene ID in annotation	Gene	Highlights of predicted function	Gene description	Category
3	maker-ptg000361l-augustus-gene-0.2-mRNA-1	ACTR3b	F-actin dynamics/polymerization	ARP3 actin-related protein 3 homolog B	Cytoskeletal and sarcomeric
3	maker-ptg000509l-snap-gene-4.4-mRNA-1	NHS	Regulator of actin remodeling	Nance–Horan syndrome protein	Cytoskeletal and sarcomeric
3	maker-ptg001966l-augustus-gene-18.42-mRNA-1	KCNK6	Outward rectification in a physiological potassium gradient and mild inward rectification in symmetrical potassium conditions	Potassium channel subfamily K member 6	Transmembrane ion transport
3	maker-ptg000697l-snap-gene-6.109-mRNA-1	KCNQ5	Voltage-gated potassium channel	Potassium voltage-gated channel subfamily Q member 5	Transmembrane ion transport
5	maker-ptg000265l-est_gff_est2genome-gene-6.33-mRNA-1	KCNJ2	Inwardly rectifying potassium channel	Inward rectifier potassium channel 2	Transmembrane transport
5	maker-ptg001740l-augustus-gene-1.112-mRNA-1	KLF5	Transcription factor, which might regulate potassium channel genes	Krüppel-like factor 5	Transcription factor

Assigning DEGs with increasing expression pattern to GO terms revealed 19 significantly enriched (Fisher's exact P < 0.05) GO terms (supplementary fig. S4 and table S5, Supplementary Material online). Eighteen genes were assigned to the enriched GO term “Golgi apparatus” and 13 to “ion transport.” Among the genes with decreasing expression pattern (group 3), 41 GO terms were significantly enriched (supplementary fig. S5 and table S6, Supplementary Material online). Twelve of the genes were assigned to the GO term “axon guidance,” which yielded the lowest P-value. There were also several enriched terms that might be functionally related to the EOD, e.g. membrane, Golgi membrane and apparatus, calcium ion binding, and ATP binding.

Allele-Specific Expression in F1 Hybrids

Two cohorts of F1 hybrids with 1 short-duration EOD (com × tsh) and 1 medium-duration EOD (com × rhy) were analyzed in our study. In total, we identified fixed single nucleotide polymorphisms (SNPs; homozygous in parental species) in 177 genes differentially expressed in EO and in 52 differentially expressed in SM in the hybrid com × rhy. For the hybrid com × tsh, the respective SNP numbers were 77 in genes differentially expressed in the EO and 36 in genes differentially expressed in the SM (Fig. 5a). For each of these genes, we calculated the allelic read proportion of the allele stemming from the parental species com (as identified by the fixed SNPs), averaged over the specimens of the respective hybrid cohort. In general, most genes exhibit an equal expression of both parental alleles, with more genes having a com proportion near 0.5 (Fig. 5a). Among the genes with differentially expressed alleles, alleles stemming from com had an overall tendency toward higher expression, compared to the alleles from rhy or tsh, in both EO and SM from 2 hybrid cohorts (Fig. 5a).

Fig. 5.

Allele-specific expression in EO and SM among 2 hybrid cohorts C. compressirostris (com) × C. rhynchophorus (rhy) and C. compressirostris (com) × C. tshokwe (tsh). a) Gene density (y-axis) from two different tissues in each hybrids. The x-axis shows the expression proportion of the allele stemming from one parental species (com). Numbers above the bars represent the number of genes in the respective proportion ranges. b) Proportion of two alleles from parental species in four genes related to ion transport and membrane (KCNJ2, SCN4AA, CHRND and TSPAN7b) in the EO. The x-axis represents the different individual samples of the corresponding hybrid cohort.

In order to understand the allelic expression imbalance (AEI), we counted the number of genes with more than 0.6 proportion of 1 parental allele proportion for all individuals in each analyzed hybrid set (Chang et al. 2002). In total, we identified 17 and 7 genes with AEI in EO and SM of the hybrid com × rhy, respectively; 2 and 1 such genes were identified in EO and SM of the hybrid com × tsh, respectively (supplementary table S7, Supplementary Material online). In all the genes with AEI, the allele from com showed a higher expression proportion, except for the gene KCNJ2 in the hybrid com × rhy where allele expression was biased in the opposite direction (average proportion of com allele was 0.16; Fig. 5b). We inferred amino acid sequences from the transcript sequences of the KNCJ2 gene from com, tsh, and rhy. The inferred protein sequences between com and tsh were identical, but rhy showed 2 amino acid substitutions at sites 60 (corresponding to the fixed SNP we identified in hybrids) and 198 (supplementary table S8, Supplementary Material online). The amino acid substitution at site 60 was considered benign in the Polyphen2 analysis (Adzhubei et al. 2010), while the substitution at site 198 may have changed the protein function in rhy (inferred as probably damaging; supplementary table S9, Supplementary Material online).

We also identified AEI in the gene SCN4aa in the EO of com × rhy hybrids (average proportion of com allele was 0.66; Fig. 5b). This proportion of the com allele is much higher than in the hybrid com × tsh, where the com allele of the SCN4aa only had a proportion of 0.46 in the EO. In the EO of hybrid com × rhy, AEI was identified in the gene CHRND, which might relate to ion channel gating (Fig. 5b). The TSPAN6b gene, encoding for an integral component of the plasma membrane, was also identified to exhibit a significant AEI in the EO of hybrid com × tsh (average proportion of com allele was 0.86; Fig. 5b).

Discussion

Convergent Gene Expression in Different Electric Fish Lineages

The myogenic EO has convergently evolved 6 times in fishes. Even though the EOs show great differences in electrocyte morphology among independently evolved electric fish lineages, particular genes exhibit similar transcriptional expression patterns in the EO, relative to SM (Gallant et al. 2014).

Several EO-specific candidate genes that we identified in Campylomormyrus were also overexpressed in the EO of other electric fish lineages, possibly indicating convergent expression pattern evolution in electric fish. This is particularly apparent in genes related to sodium and potassium currents. For instance, the Na_v channel gene SCN4aa, considered to be very important in regulating the sodium current to electrocytes, was previously found overexpressed in the EO of different electric fish lineages, i.e. Siluriformes, Gymnotiformes, and Mormyridae other than Campylomormyrus (Wang and Yang 2021). The FGF13a that regulates this channel was consistently overexpressed in those electric fish species as well (Gallant et al. 2014). Interestingly, we identified another upregulated ortholog (FGF12) in the EO, which may have a similar function. In addition, multiple isoforms of sodium/potassium-ATPase α- and β-subunits and several transcription factors (SIX2a, HEY1) were found to be convergently upregulated in the EO among these electric fish lineages (Gallant et al. 2014), a pattern confirmed for Campylomormyrus in our study.

Overexpression of another transcription factor (MEF2a) and of the calcium binding gene S100b is characteristic for mormyrid EOs, i.e. Paramormyrops, Brienomyrus (Gallant et al. 2012, 2017), and Campylomormyrus (this study). We recently found KCNA7a to be tandemly duplicated in Campylomormyrus (Cheng et al. 2023) and Paramormyrops (by reanalysis of the genome provided in Gallant et al. 2017). This tandem duplication might be exclusive to mormyrid fishes, as we did not find it in available genomes neither of other electric fishes nor in Scleropages (a nonelectric fish closely related to mormyrids; data not shown). In our study, both gene copies KCNA7a_1 and KCNA7a_2 were consistently upregulated in the EO (KCNA7a_2 showed even higher expression than KCNA7a_1). KCNA7a was inferred to be under positive selection in the transmembrane helix 3 and 4 linkers and is considered to relate to the differences in EOD duration among Brienomyrus and Gymnarchus (Swapna et al. 2018).

The NDRG3 gene (N-Myc Downstream-Regulated Gene 3) exhibited a remarkable overexpression in the EO of Campylomormyrus. Interestingly, the phosphopeptides encoded by an ortholog (NDRG4) were highly enriched in the EO of the strongly discharging gymnotiform electric eel (Electrophorus electricus; Traeger et al. 2017), indicating high expression level of this gene. In addition, NDRG4 has been identified in zebrafish as a novel neuronal factor essential for sodium channel clustering at the nodes of Ranvier, the only places where action potentials are regenerated (Fontenas et al. 2016). The function of NDRG3 in the nervous system has rarely been investigated. The NDRG3 protein can interact with extracellular signal-regulated kinases (ERK1/2; Lee et al. 2015), which regulate K_v4.2 in the dendrites of hippocampal CA1 pyramidal neurons (Schrader et al. 2006; Gupte et al. 2016) as well as the Na_v1.7 channel (Stamboulian et al. 2010). In porcine as well as human lens, ERK1/2 is activated by the TRPV1 ion channel (Mandal et al. 2019), which was also overexpressed in EOs in Campylomormyrus (supplementary table S1, Supplementary Material online).

Gene Expression Specificity in Campylomormyrus

The EO in mormyrids is derived from myogenic tissue, which transitions from a motoric/sarcomeric organization of muscle fibers to a continuous tube of electrocytes parallel to the spinal cord (Denizot et al. 1982). This transition process during the ontogeny of the EO involves cell size, morphology, and physiology and is still only partially understood. Some genes encoding for sarcomeric proteins, e.g. troponin I isoforms, myosin heavy chain, and tropomyosin, are overexpressed in the EO of the mormyrid Brienomyrus brachyistius (Gallant et al. 2012), providing a preliminary insight into the developmental transition from SM to EO. In the EO of Campylomormyrus, however, we rarely found those genes upregulated. Instead, the upregulated actin-related genes in Campylomormyrus were more related to F-actin dynamics and included several unconventional myosins (Table 1; supplementary table S1, Supplementary Material online). The 4 paralogous transcription factors MEF2a to MEF2d are responsible for the transcriptional activation of muscle-specific genes in the early specification of SM (Black and Olson 1998). Whereas MEF2a was overexpressed in the EOs of both Brienomyrus (Gallant et al. 2012) and Campylomormyrus, a further paralog MEF2b was overexpressed only in Campylomormyrus (our study). The difference in expression of F-actin-related/sarcomeric proteins and MEF2 transcription factors between 2 mormyrids genera suggests that the developmental transition in the EO might be different or, in other words, that the organization of the F-actin system in electrocytes may vary across mormyrids. It has to be analyzed further whether these differences in the organization of the F-actin cytoskeleton concern the sarcomeric structure, the stalks of electrocytes, or both.

In addition, 2 paralogs of inwardly rectifying potassium channel (K_ir) genes, KCNJ2 and KCNJ9, were overexpressed in EOs of Campylomormyrus, along with KCNQ5 and KCNK2. The mechanisms regulating potassium channels’ expression in electric fish are still unknown. We identified 1 transcription factor, KLF5, that showed a high overexpression in EOs. In Drosophila, Krüppel is involved in the regulation of potassium channel expression. In case of a loss of Shal (KCND) potassium channel in Drosophila, Krüppel expression is induced and upregulates expression of Shaker (KCNA) and slowpoke (K_ca) potassium channels (Parrish et al. 2014). Remarkably, Shal (KCND) potassium channel is also not expressed in our studied Campylomormyrus species/hybrids. We thus suppose that the EO differs from SM by the expression of a unique set of potassium channels that may contribute to the shape of the EO’s action potential and thus the shape of the EOD signal. Moreover, we propose the KLF5 gene to represent a transcription factor that drives the expression of regulating potassium channels in EO.

Our study has further revealed that different paralogs from the solute carrier family are active in EOs. Solute carriers form a group of membrane transport proteins located in various cellular membrane systems, which transport diverse substrates including amino acids, oligopeptides, inorganic cations, and anions (He et al. 2009). We have found several genes of this gene family overexpressed in the EO that transport inorganic cations and anions, e.g. sodium, calcium, and chloride. Especially the gene SLC24a2 was highly overexpressed in Campylomormyrus EOs. This is a calcium/cation antiporter localized in the plasma membrane that mediates the extrusion of 1 calcium ion and 1 potassium ion in exchange for 4 sodium ions (Wang et al. 2017). Overexpression of a calcium-extruding transporter in the EO indicates that regulation of the cytosolic calcium level in electrocyte differs from that in SM. Unfortunately, we have no information yet on the distribution of the SLC24a2 protein within the electrocytes and whether this calcium transporter is confined to a distinct region of the cell to mediate local regulation of the calcium level.

Differential Gene Expression with Respect to EOD Duration Divergence among Campylomormyrus Species

Campylomormyrus species produce species-specific EODs; their duration varies in a 100-fold range across species. The EOD is assumed to be mediated by sodium and potassium currents across the plasma membrane (Stoddard and Markham 2008). The depolarizing phase of an action potential is primarily produced by sodium influx. The repolarization phase is—along with a gradual decreasing sodium influx—affected by the orchestrated activities of delayed rectifier and inward rectifier potassium channels (Nass et al. 2008; Stoddard and Markham 2008). We suppose that species producing EODs of different duration may be equipped with different channel types or channel orthologs with different properties. However, certain other mechanisms, such as different cell morphology, may also contribute to the EOD duration diversification.

The PCA from RNA-seq data showed a clear association between the overall gene expression and the EOD duration pattern (Fig. 3a). Based on the preliminary PCA and LRT result, we identified several genes that might contribute to EOD duration diversification in Campylomormyrus, including the potassium channel genes (KCNJ2, KCNK6, and KCNQ5), actin-related genes (ACTR3b, NHS), and transcription factor KLF5 (Table 2).

The gene KCNK5 was found to be upregulated in Paramormyrops (producing a short EOD) compared to the species with an elongated EOD (Losilla et al. 2020). In Campylomormyrus, the expression of another paralog KCNK6 was also higher in species with short EOD. Two-pore potassium channels (K_2p) usually generate an outward potassium current and are also known as potassium leak channels. When silencing the KCNK6 gene in the human heart, the action potential duration is prolongated (Chai et al. 2017). Another voltage-gated potassium channel gene KCNQ5 was decreasingly expressed in elongated EOD Campylomormyrus species. It forms M-type potassium current, a slowly activating and deactivating potassium conductance that works in determining the subthreshold electrical excitability of neurons (Hibino et al. 2010). The lower expression of both potassium channel genes in elongated EOD species will probably decrease the outward potassium current and consequently prolongate EOD repolarization.

The gene KCNJ2 was increasingly expressed in elongated EOD species. It encodes for an K_ir, with the greater tendency to allow potassium ions to flow into a cell rather than out of a cell (Hibino et al. 2010). The inward potassium current stabilizes the resting membrane potential of the cell and modulates the cardiac repolarization processes (Hibino et al. 2010; Li et al. 2017). This inward rectifier channel-mediated potassium current is responsible for shaping the initial depolarization and final repolarization of the action potential in human cardiomyocytes (Dhamoon and Jalife 2005; Jeevaratnam et al. 2018).

Regarding allele-specific expression, the 26 genes with AEI (supplementary table S7, Supplementary Material online) in F1 hybrids were related to ion transport, plasma membrane, and myofibrils in the EO and actin in the SM (supplementary table S10, Supplementary Material online). The imbalance could be caused by cis-regularity differences among the parental species in the upstream transcription unit (e.g. a promoter; Landry et al. 2007). There was a tendency toward higher expression of com alleles in the EOs and SMs of 2 analyzed hybrid cohorts (Fig. 5a; supplementary table S7, Supplementary Material online). However, the phenotype of the EOD waveform in each hybrid is closer to the other parental species. This points toward some genes playing key roles in regulating the EOD waveform in the hybrids. The gene KCNJ2 showed AEI in com × rhy hybrids, which was the only gene with AEI and for which the rhy allele was preferentially expressed (Fig. 5b). The EOD in the adult hybrids com × rhy was of intermediate duration (4 ms), and the shape and waveform resemble the subadults’ EOD in rhy. Both the EOD phenotype and the AEI in KCNJ2 were hence closer to the parental species with the elongated EOD, i.e. rhy. The expression of KCNJ2 in the EO among the purebred species also increased with increasing EOD duration, e.g. the expression in rhy is higher than in com. This suggests that the KCNJ2 gene might be under cis-regulation, and it should be a powerful candidate gene involved in the regulation of EOD duration in Campylomormyrus.

In addition, the KCNJ2 gene in the species rhy (very long EOD) exhibits 2 nonsynonymous substitutions, one of which predicted to cause a functionally relevant amino acid substitution (at site 198; supplementary tables S8 and S9, Supplementary Material online). Interestingly, the same substitution at site 198 is present in another species with very long EOD (C. numenius, EOD duration 40 ms), while it is absent in other Campylomormyrus species with short or medium EOD, which resemble the amino acid sequence of com and tsh (Cheri, Cheng, and Tiedemann, unpublished results). Campylomormyrus numenius and rhy are phylogenetically close (Lamanna et al. 2016), such that the shared amino acid substitution could also reflect phylogenetic affinity. Nonetheless, the found amino acid substitution with inferred functional relevance could relate to the evolution of very long EODs in Campylomormyrus. Then, the KCNJ2 gene could modulate EOD duration by a combination of expression level and functional protein sequence alteration. In summary, this study identifies the KCNJ2, KCNK6, and KCNQ5 genes, possibly in combination with other genes (e.g. KLF5, ACTR3b, and NHS) as strong candidates underlying EOD duration diversification in the weakly electric fish genus Campylomormyrus. The diverged EODs likely affect the food spectrum and are used for mate recognition. This potential dual function in disruptive natural selection and prezygotic reproductive isolation would rank the EOD as a “magic trait,” which may have promoted the ecological (probably sympatric) speciation and radiation of Campylomormyrus in the Congo River.

Materials and Methods

Animals, RNA Isolation, Library Preparation, and Sequencing

Three adult specimens of C. tshokwe were collected at Brazzaville/Republic of the Congo in the Congo River in 2012 and stored in RNAlater in −80 °C. Five adult specimens from each of the other 2 species (C. compressirostris, C. rhynchophorus) and 2 hybrids (C. compressirostris ♂ × C. rhynchophorus ♀, C. compressirostris ♂ × C. tshokwe ♀) were artificially bred and raised at the University of Potsdam. All specimens except for C. tshokwe were anesthetized by a lethal dose of clove oil and dissected on cold 99% ethanol. EO and SM tissues from each specimen were flash frozen in liquid nitrogen and further preserved in −80 °C. In total, we collected 3 samples of both EOs and SMs from C. tshokwe and 5 samples of both EOs and SMs from the other 4 species/hybrid cohorts in this study.

The RNA isolation was performed in all the EO and SM samples using QIAGEN RNeasy Fibrous Tissue Mini Kit. Total RNA concentration was estimated using a NanoDrop 1000 spectrophotometer (Thermo Fisher Scientific, Germany), and RNA quality was checked with an Agilent Bioanalyzer 2100 (Agilent Technologies, United States). mRNA enrichment was performed by poly(A) capture from isolated RNA using NEXTflex Poly(A) Beads. Strand-specific transcriptomic libraries were built using NEXTflex Rapid Directional RNA-Seq Kit (Bioo Scientific, United States) based on the manufacturer's instructions.

Libraries were sequenced as 150 bp paired-end reads by Illumina HiSeq 4000 sequencing system at a commercial company (Novogene). Raw reads have been deposited in the National Center for Biotechnology Information (NCBI) Gene Expression Omnibus (accession number GSE240783). We trimmed the adapter sequences and low-quality reads using a 4 bp sliding window with a mean quality threshold of 25 and a minimum read length of 36 bp by Trimmomatic v0.39 (Bolger et al. 2014). Read quality, before and after read filtering, was measured by FastQC v0.11.9 (Andrew 2010).

Differential Gene Expression Analysis

The quality-filtered reads from EOs and SMs were mapped to the C. compressirostris genome (Cheng et al. 2023) using RSEM (Li and Dewey 2011) for gene-level quantification estimation. The estimated counts were imported into R/Bioconductor with the tximport package (Soneson et al. 2015), which produced count matrices from gene-level quantification by taking the effective gene length into account (Paraskevopoulou et al. 2020). Low-count (≤10) and low-frequency (not present in at least 2 replicates) genes were removed. We performed a PCA from filtered and log-transformed counts. One SM sample from C. compressirostris was removed from this study, as its overall gene expression showed a deviant unusual pattern in the PCA.

We forwarded the normalized count matrices to DESeq2 (Love et al. 2014) to infer expression differences among EO and SM in each species/hybrid cohort, respectively. We used a false discovery rate threshold of 0.05 to correct for multiple testing. The DEGs were identified with |log₂FC| > 1 and P < 0.05. In order to detect the EO-specific gene expression pattern, we used Venn diagrams (Chen and Boutros 2011) to visualize the shared DEGs (up- and downregulated separately) among 3 purebred species and 2 hybrid cohorts.

The shared DEGs were annotated against the NCBI nr database by blastx with an e-value cutoff of 1E⁻¹⁰. In addition, the up- and downregulated DEGs in the EO were used to perform a GO enrichment analysis (Dennis et al. 2003).

RNA-seq Data Clustering by EOD Duration

The PCA plotting from log-transformed count matrices showed a clear pattern by the length of EOD (Fig. 3a). To identify genes with an expression pattern associated to EOD duration, we used DESeq2 to perform a LRT (Love et al. 2015; LRT in the DESeq2 package). This test compares how well a gene's count data fit a “full model” compared to a “reduced model” (Bendjilali et al. 2017). Our full model was an equation: full = ∼duration × tissue. The duration is the length of the EOD in each purebred species, and tissue is the type of sample (EO or SM). The reduced model excluded the interaction between duration and tissue: reduced = ∼duration + tissue. Genes with P_adj < 0.05 were considered to fit the “full model.” We used the degPatterns function from the “DEGreport” package to cluster different groups with particular expression pattern using those significant genes across samples (Pantano 2019), with time = “duration” and col = “tissue.”

The generated groups of different gene expression pattern across EOD duration were analyzed to identify genes with an expression pattern association with EOD duration. We hence focused on those groups fulfilling the following 2 criteria: (i) the gene expression level in EO is higher than SM in all F0 species, and (ii) the gene expression level in EO across EOD duration showed a consistent increasing or decreasing pattern.

The identified genes with increased and decreased expression relative to EOD duration were blasted against nr database using an e-value cutoff of 1E⁻¹⁰. In addition, a GO term analysis was also performed for these genes.

Allelic-Specific Expression Analysis

The F1 hybrid contains 2 sets of subgenome from 2 parental species. Examination of allele-specific expression can be applied to detect the allelic imbalance in transcription in heterozygous F1 hybrids. We only focused on transcripts of genes with biallelic SNPs fixed among the respective F0 parental species (hence, heterozygous only in F1 hybrids and homozygous in parental species).

We mapped the trimmed and filtered RNA-seq from 5 species/hybrids (in EOs and SMs, respectively) to the C. compressirostris genome using STAR v2.7.7 (Dobin and Gingeras 2015). The generated bam files were sorted according to the coordinates by SAMtools v1.15 (Danecek et al. 2021). Variant calling was performed by BCFtools v1.9 (Danecek et al. 2021) in EOs and SMs, respectively, using the command “bcftools mpileup –f REFERENCE LIST_OF_BAM –Ou | bcftools call –mv –Ob –o BCFFILE,” where the REFERENCE, LIST_OF_BAM, and BCFFILE were the CDS sequence name of the C. compressirostris genome, the list of bam files, and the output bcf file name, respectively.

After the variant calling, we performed the following steps to identify the fixed parental biallelic SNPs for each hybrid set. Firstly, we excluded the uncalled variants and only preserved biallelic SNPs using the command “bcftools view --exclude-uncalled -m2 -M2 BCFFILE > CALLING_AD,” where the CALLING_AD was the allelic depth for the final biallelic SNPs. Secondly, we discarded SNPs where the variant calling score at QUAL field was lower than 70, and allelic depth was lower than 10 in both alleles. Finally, we obtained high-quality SNPs, at which both parental species were homozygous and fixed for a different allele.

For each hybrid, we calculated the expression proportion of the allele from C. compressirostris in EO and SM, respectively. We calculated the average proportion and its 95% confidence limits across biological replicates (and over SNPs in case of more than 1 SNP per locus; Fig. 5; supplementary table S7, Supplementary Material online). Genes with C. compressirostris allele proportions <0.4 or >0.6 in the transcriptomes were considered exhibiting an imbalanced expression. We further applied GO analysis on the genes with AEI for EO and SM, respectively.

Supplementary Material

Supplementary material is available at Molecular Biology and Evolution online.

Acknowledgments

We thank Dr. Linh Nguyen for fish breeding and raising. All figures were edited with Biorender.com.

Author Contributions

R.T. conceived and supervised this study and provided financial and logistical support. F.C. performed lab work, analysis, and manuscript writing with the input from R.T., A.B.D., O.B., F.K., and S.A.-S.

Funding

This project is funded by the University of Potsdam.

Conflict of Interest

None declared.

Data Availability

Sequence data have been deposited at NCBI Gene Expression Omnibus under accession number GSE240783.