Skip to main content

Transcriptomics unravels molecular players shaping dorsal lip hypertrophy in the vacuum cleaner cichlid, Gnathochromis permaxillaris



Teleosts display a spectacular diversity of craniofacial adaptations that often mediates ecological specializations. A considerable amount of research has revealed molecular players underlying skeletal craniofacial morphologies, but less is known about soft craniofacial phenotypes. Here we focus on an example of lip hypertrophy in the benthivorous Lake Tangnayika cichlid, Gnathochromis permaxillaris, considered to be a morphological adaptation to extract invertebrates out of the uppermost layer of mud bottom. We investigate the molecular and regulatory basis of lip hypertrophy in G. permaxillaris using a comparative transcriptomic approach.


We identified a gene regulatory network involved in tissue overgrowth and cellular hypertrophy, potentially associated with the formation of a locally restricted hypertrophic lip in a teleost fish species. Of particular interest were the increased expression level of apoda and fhl2, as well as reduced expression of cyp1a, gimap8, lama5 and rasal3, in the hypertrophic lip region which have been implicated in lip formation in other vertebrates. Among the predicted upstream transcription factors, we found reduced expression of foxp1 in the hypertrophic lip region, which is known to act as repressor of cell growth and proliferation, and its function has been associated with hypertrophy of upper lip in human.


Our results provide a genetic foundation for future studies of molecular players shaping soft and exaggerated, but locally restricted, craniofacial morphological changes in fish and perhaps across vertebrates. In the future, we advocate integrating gene regulatory networks of various craniofacial phenotypes to understand how they collectively govern trophic and behavioural adaptations.

Peer Review reports


Teleost fishes show striking adaptive diversity in their craniofacial anatomy, which reflects an equally striking variety of ecological and trophic specializations. It is therefore not surprising that this phenotypic diversity and its adaptive background has attracted the attention of developmental, molecular and evolutionary biologists, beyond model species such as zebrafish and medaka [1, 2]. In the last two decades, an impressive set of molecular players participating in the development and morphogenesis of craniofacial skeletal structures, including their interconnecting signalling pathways, have been described in teleost fishes [2,3,4,5]. Much less is known about the morphogenic molecular players and underlying signals of craniofacial soft tissues, that also exhibit varied adaptive peculiarities. In cichlid fishes, a novel model system for adaptive radiation, the molecular mechanisms underlying exaggerated soft-tissue morphologies in lip, nose and nuchal hump phenotypes have been addressed only recently [6,7,8,9,10,11].

One such exaggerated craniofacial soft-tissue phenotype is the overgrowth of the lip tissue, the so-called thick-lipped phenotype, which has been observed in cichlid species in lakes spanning different continents [7, 8, 12]. The repeated evolution of the hypertrophic lip phenotype, in both African and Central American cichlids, has been associated with the dietary adaptation to suck elusive invertebrates out of narrow crevices in rocky habitats [8, 13, 14]. Genetic studies so far have suggested numerous loci across the genome whose small additive effects can be linked to a variety of thick-lipped phenotypes [7, 8, 12]. These findings indicate that the thick-lipped phenotype is a complex trait [12]. The exaggerated thickening of the lip might not follow a uniform pattern across the entire upper and lower lip and the extent of lip hypertrophy is significantly reduced in captivity under flake food diet, when the fish are not using this foraging strategy, so that lip hypertrophy seems highly dependent on foraging performance, making it a phenotypically plastic trait [12]. Both genetic and phenotypic plasticity take part in the adaptive morphogenesis of this phenotype [6, 15].

Another interesting example of lip hypertrophy is found in the Lake Tanganyika cichlid adaptive radiation, in the species Gnathochromis permaxillaris. It has a shovel-like snout with a peculiarly thickened part of the upper lip [16]. G. permaxillaris is a benthivorous deepwater cichlid of the tribe Limnochromini [16] that develops a unique hypertrophy in the most anterior part of its upper lip (Fig. 1). The species slowly swims directly over the mud surface with its protruding snout hovering at the surface and at the same time opens the protruding mouth underneath, while ingesting mud that is then filtered through the gill arches. Its common name ‘the Vacuum cleaner cichlid’ is an adequate description of this feeding mode. Thereby, the protruding upper lip with its small hypertrophy at its tip seems to be adaptive by boosting the efficiency of filtration (unpublished behavioural observations by Heinz H. Bücher; see Video S1). Interestingly, when fish are raised in captivity, individuals develop this phenotype to a lesser extent (Peter Henninger, personal communication), again suggesting some degree of phenotypic plasticity. Nothing is known so far about the molecular basis of lip morphogenesis in G. permaxillaris.

Fig. 1
figure 1

Two East African cichlid species of the tribe Limnochromini from Lake Tanganyika used in this study. The areas from which the soft tissue samples are taken in lips are delineated by coloured dash lines; UL/ul and LL/ll refer to upper and lower lip, respectively

Additional file 3: Video S1. Foraging behaviour of the vacuum cleaner cichlid in its natural habitat.

In this study, we aimed to identify genes playing a role in the formation of the dorsal part of the upper lip, using RNA-sequencing. To this end, we profiled gene expression differences between the hypertrophic part of the upper lip versus the posterior part of the upper lip and the most anterior part of the lower lip, the latter two do not exhibit hypertrophic overgrowth in wild-caught young male adults of G. permaxillaris (Fig. 1). In order to add another layer of filtering, we also conducted an inter-specific comparison using a closely related species in the tribe Limnochromini, Greenwoodochromis bellcrossi [16,17,18], which does not have a protruding upper lip but displays hypertrophy in the entire parts of both the upper and lower lips (Fig. 1). Using the regulatory sequences of differentially expressed genes in the hypertrophic lip region, we predicted their potential upstream transcriptional regulators, and applied qPCR expression analysis on the most interesting candidate genes to validate our results from the RNA-seq analysis. We identified a gene regulatory network, potentially associated with the formation of a locally restricted hypertrophic lip in the Limnochromini. Our results provide a foundation for future investigations of molecular players shaping exaggerated but locally restricted morphological changes in facial elements of fish and vertebrates. Moreover, these findings can be further used for molecular comparisons between wild and captive bred G. permaxillaris in order to investigate the mechanisms underlying the plasticity of the lip hypertrophy in fish. Recent attempts have been made to unravel molecular factors underlying plastic responses to different mechanical stimuli in cichlid jaw skeleton [19,20,21], but studies focusing on craniofacial soft tissues are lacking. Our results pave the way for interesting future functional gene characterisation as well.


RNA-seq, differential gene expression and downstream analyses

The total RNA sequencing resulted between 3.84 and 8.13 million reads per sample and after removal of low quality reads, each sample had between 3.82 and 8.10 million reads (Table S1). The low reads for some of the samples limit the results of our study to differentially expressed (DE) genes that have relatively high expression levels, in other words, it is likely that some of the low expressed genes are not detected among the list of DE genes due to the low coverage of some samples. This can be particularly the case for comparisons involving G. bellcrossi since the two samples with less than 4 million reads are from this species. The raw sequence reads have been submitted to the Sequencing Read Archive (SRA) of NCBI (accession number: PRJNA694145). The comparisons of the lip regions of G. permaxillaris resulted in 106 DE genes between the hypertrophic part of the upper lip (UL-1) versus both the non-hypertrophic part of the upper lip (UL-2) and the lower lip (LL) (Fig. 2A). Among these, 56 genes showed increased expression, whereas 47 genes had reduced expression in UL-1 (Fig. 2B and C). The heatmap clustering of the DE genes revealed two major branches in the group of down-regulated genes (Fig. 2B), whereas only one major branch is observed in the group of up-regulated genes (Fig. 2C).

Fig. 2
figure 2

Differentially expressed genes in the hypertrophic region of dorsal lip in Gnathochromis permaxillaris. A Venn diagram representing 106 genes showing differential expression between the hypertrophic region. (UL-1) and other regions of the lip. Dendrogram clusters of genes with lower (B) and higher (C) expression in UL-1 region in G. permaxillaris compared to the other lip regions. Red and green shadings indicate higher and lower relative expression, respectively

We also conducted a second filtering step in order to check whether among the identified list of 106 genes above, any might play roles in lip hypertrophy in another closely related species. To do this, we first compared UL-1 from G. permaxillaris to the upper lip (ul) of G. bellcrossi, and since both regions have hypertrophy the DE genes between them were filtered out. Next, we compared LL from G. permaxillaris to the lower lip (ll) of G. bellcrossi, and this time, since ll has hypertrophy but not LL, we checked for genes which were shared with the list of 106 genes (from the within species comparisons). We only found 24 DE genes to be shared between the intra- and inter-specific comparisons, suggesting their potential roles in lip hypertrophy across the species (Fig. 3). Out of these, 9 genes showed higher and 15 showed lower expression in the hypertrophic lip region between the species (ll vesus LL), and similar expression patterns were observed for the within-species comparisons. Among these 24 genes, at least 5 genes, apoda, cyp1a, lama5, rasal3 and trim37, have been already found to be involved in lip morphogenesis in other vertebrates (Table 3).

Fig. 3
figure 3

Differentially expressed genes in the hypertrophic region of the lips shared between G. permaxillaris and G. bellcrossi. A. Venn diagram representing 24 genes showing differential expression between the hypertrophic and non-hypertrophic lip regions between the species. B. Dendrogram clusters of the 24 shared DE genes in LL region in G. permaxillaris compared to the ll lip region in G. bellcrossi. Red and green shadings indicate higher and lower relative expression, respectively

Using the list of 106 DE genes, we conducted a gene ontology enrichment analysis, and amongst the significantly enriched GO terms, biological processes related to cell movement, adhesion and growth had the highest enrichment ratios (Fig. 4A). The enrichment for genes related to locomotion raises the possibility of functional specificity for the hypertrophy of the most anterior part of the upper lip of G. permaxillaris. We also found enrichments of molecular processes related to GTPase and Ras protein signalling for the list of DE genes. Subjecting the same gene list to an interactome analysis revealed only a single large interconnected network of genes. However, the network contained several interesting genes, and particularly transcription factors (TFs), with potential role in tissue overgrowth and cellular hypertrophy (Fig. 4B). We also performed overrepresentation analysis of TF binding motifs on the regulatory sequences of the DE genes using the MEME algorithm [22]. This resulted in five enriched motifs present in the regulatory sequences of at least 30 out of 106 DE genes (Table 1). In search for similarities of the identified motifs with known TF binding sites in vertebrates, we found at least two interesting TF candidates, foxp1 and heb/tcf12, with overrepresented binding sites. Furthermore, both of these TFs showed regulatory connections with some of the DE genes in the identified interactome network, with foxp1 having more connections than tcf12 (Fig. 4B).

Fig. 4
figure 4

Downstream functional analyses of differentially expressed genes. A. Gene ontology enrichment analysis for biological processes, using the shared 106 differentially expressed genes, conducted with DAVID tool. B. Predicted functional associations/interactions between the differentially expressed genes based on zebrafish databases in STRING v10 ( The differential expression of the genes specified with red rectangles are confirmed by both qPCR and RNA-seq

Table 1 Predicted binding sites for potential upstream regulators of the differentially expressed genes identified by RNA-seq. PWM ID indicates positional weight matrix ID of a predicted binding site and E-values refer to matching similarity between the predicted motif sequences and the PWM IDs. The count implies on number of genes containing the predicted motif sequence on their regulatory region

Expression validation using qPCR

The gene expression analysis by qPCR requires the identification of stably expressed reference genes [23], and our previous studies on East African cichlids, confirmed that the validation of suitable reference gene(s) is an essential step that needs to be taken in every species, every tissue and experimental condition [24,25,26,27,28]. In order to choose adequate reference gene candidates, we conducted a ranking for the genes with no expression difference (FDR = 1 in the RNA-seq comparisons) based on their coefficient variation (CV) throughout all the samples. The ten genes with lowest CV were selected as candidate reference genes for validation of their expression stability by qPCR. None of the validated reference genes in previous gene expression studies of East African cichlids have appeared among the ten candidates, confirming the necessity of reference gene validation for each experimental setup [24,25,26,27,28,29]. Overall, the reference gene rankings by the three algorithms, BestKeeper, geNorm and NormFinder software showed that two of the candidates, pip4k2b and ss18, were always among the top three most stable reference genes (Table 2). Therefore, we used geometric mean of Cq values for pip4k2b and ss18 in each sample, as a normalization factor, in order to calculate the relative gene expression levels of our target genes.

Table 2 Ranking and statistical analyses of reference genes in the lip samples using three different algorithms

Within the identified DE genes in the RNA-seq results, we selected nine genes with a known role in lip morphogenesis and other related functions in vertebrates (Table 3), as well as the two predicted upstream TFs (foxp1 and tcf12) to be tested by qPCR (Fig. 5). Among these genes, apoda, fhl2 and gimap8, were also implicated to be involved in lip hypertrophy in other cichlids (Table 3), and 5 genes, apoda, cyp1a, lama5, rasal3 and trim37 appeared to be differentially expressed in both intra- and inter-specific comparisons between hypertrophic and non-hypertrophic lips. Based on the RNA-seq results, four of the selected candidate genes have shown reduced expression in the hypertrophic UL-1 (cyp1a, gimap8, lama5 and rasal3) whereas the other six genes (apoda, foxf1, fhl2, igf2b and trim37) had increased expression level in UL-1 of G. permaxillaris. The qPCR results showed that all of the genes followed similar expression patterns as found by RNA-seq, except for trim37, which showed no significant difference between the lip regions in G. permaxillaris. However, when the qPCR expression results for G. permaxillaris were compared to the two hypertrophic lip regions of G. bellcrossi, we found only apoda and fhl2 showing increased expression in all the hypertrophic lip regions across both species. The four genes with reduced expression in UL-1 for G. permaxillaris, also showed reduced expression in the hypertrophic lip regions of G. bellcrossi. Among the two predicted TFs, only foxp1 showed consistent difference between hypertrophic and non-hypertrophic lip regions across both species, with reduced expression in the hypertrophic lip regions. This indicates potential transcriptional repressor effects of foxp1 on the downstream genes in the hypertrophic lip region. Altogether, the qPCR results show consistency with RNA-seq results, confirming the validity of our transcriptome data analysis in this study.

Table 3 Differentially expressed genes in the hypertrophic region of dorsal lip in Gnathochromis permaxillaris with related functions in vertebrates
Fig. 5
figure 5

Expression analysis of a selection of candidate genes using qPCR. The bars represent means and standard deviations of RQ values for five biological replicates in each lip region. Circles above bars indicate significantly elevated expression (P < 0.05) in comparisons between the lip regions (i.e., compared to the bar matching the colour code of the circle)


Delineating the molecular basis of adaptive morphologies is essential to understand how they develop and evolve, and how they may contribute to adaptive radiation. In cichlids, lip hypertrophy is considered a craniofacial morphological novelty, which might promote incipient sympatric speciation through divergent selection (on feeding performance) and non-random/assortative mating [13, 15, 44]. The repeated evolution of hypertrophic lip phenotypes in cichlids inhabiting different lakes provides an opportunity to study the role of craniofacial soft tissues in adaptive radiations [8, 13, 44,45,46], as well as their underlying molecular mechanisms. There are only few studies that addressed this morphological novelty at molecular levels in cichlids [7, 8]. The lip hypertrophy phenotype is variable among cichlids, for example cichlid species studied so far show hypertrophy in both upper and lower lips [7, 8], while the selected target species in this study, G. permaxillaris, displays such a phenotype only in the anterior part of its upper/dorsal lip. This may indicate distinct molecular mechanisms involved in shaping seemingly similar craniofacial novelty across cichlids. In this study, we investigated the soft-tissue craniofacial trait of the hypertrophic lip in G. permaxillaris, in a comparative transcriptomic framework with G. bellcrossi, which are both from the Limnochromini tribe belonging to the cichlid adaptive radiation of Lake Tanganyika. Gnathochromis permaxillaris is a deepwater cichlid living over mud bottom, as far as oxygenated water reaches down, and has a unique mode of foraging.

Using differential gene expression analysis in the hypertrophic lip region as input for gene ontology analysis, we found multiple biological processes to be enriched. These processes involve cell motility, adhesion and developmental growth, as well as regulation of GTPase mediated signals, particularly, the Ras signalling pathway. A recent integrated genomic and transcriptomic study in humans has revealed that cell adhesion, cell junction and extracellular structure organizations are among the major biological processes involved in lip and cleft development and morphogenesis [47]. Furthermore, as studies on human and mouse have shown, biological processes involving developmental growth and cell proliferation are known to be the key mechanisms in upper lip development and morphogenesis [48]. On the other hand, RAS/MAPK signalling is among the well-known molecular pathways involved in development and morphogenesis of craniofacial structures in vertebrates [4, 49]. Activation of Ras signalling promotes cell proliferation, growth and survival in various tissues [50, 51], and because of these roles, several components of the Ras pathway are considered as therapeutic targets in different types of cancer [52]. In mammals, Ras signaling plays a pivotal role in skin development, dermal thickenning and skin carcinogenesis [53]. Defective activity of the Ras pathway can cause a wide range of skin anomalies, such as thickened palms and soles, redundant skin, papilloma formation, excessive proliferation of keratinocytes and increased skin folds [53]. We found components of Ras signaling to be differentially expressed including rasal3, an inhibitor of the pathway (discussed below), rrad, a direct downstream target of Ras signal [54], a scaffold protein involved in signal transduction akap12b [55], and two genes encoding enzymes with roles in this signal; arhgef18b and arhgef3l [56, 57]. Therefore, the hypertophic lip region in the anterior upper lip of G. permaxillaris can be a result of increased activity of Ras signaling in this region. However, further functional studies (such as a protein manipulation method recently used in cichlids to investigate regional activity of a growth signal [10]) are required to find out the molecular reason for the anotomically limited activation of this signal only in the anterior part of the upper lip in G. permaxillaris. Interestingly, Ras signaling is also implicated among the pathways mediating molecular effects of environmentally induced mechanical stimuli in mammalian cells, raising the possiblity of its involvement in regulation of the plasticity of the lip phenotype [58].

We found that many of the genes with increased expression in the hypertrophic region of dorsal lip were already demonstrated to play a role in lip morphogensis and pathobiology in other vertebrates, including alx1 [59], alx3 [60], angptl2 [61], arid3a [62], crip2 [63], ddr2 [64], dpysl2 [65], itga5 [66], lmna [67], mgp [68], rgs5 [69], rhob [69, 70], rxfp2 [71], sostdc1 [69, 72], thbs3 [73], vcan [74, 75], and vtn [76]. Among the transcriptionally repressed genes, we also found candidates with functions which were previously implicated in defective lip morphogenesis in other vertebarates such as cd96 [77], ep300 [78, 79], fgfrl1 [80], frzb [81], hectd1 [82, 83], mmp13 [76], slc16a6 [84], and syne1 [85]. Since most of these studies were conducted on mammals, our results suggest that a similar set of genes might be involved in lip morphogenesis, not just in teleosts, but across vertebrates. Moreover, considering that all the previous studies in mammals linking these genes individually to lip morphogenesis, our findings for the first time show that these genes might be co-regulated or have regulatory interactions in an interconnected network. Thus, further functional studies are required to investigate their specific role in morphological divergence of soft tissues in fish.

Among the genes with reduced expression in the hypertophic lip regions, we validated four genes with qPCR, cyp1a, gimap8, lama5 and rasal3, and found all of them to show a similar expression reduction pattern in all of the hypertophic lip regions of both species. Cytochrome P450 (CYP) 1 alpha, cyp1a, encodes an enzyme with an important role in the cytochrome p450 xenobiotic metabolism and the synthesis of steroids and other lipids. Cyp1a is a downstream target for AHR, RAS and Wnt/β-Catenin signaling pathways [86, 87], and all of these signals are demonstrated to play important roles in morphogenesis and adaptive radiations of craniofacial elements in teleost fishes [4]. Differential regulation of cyp1a is implicated in craniofacial morphological divergence in fish [88]. Differential regulation of AHR, RAS and Wnt/β-Catenin signals, as well as cyp1a, is also involved in the defective formation of the lip and palate in vertebrates [31, 32, 89,90,91]. We found reduced expression of cyp1a in the hypertophic lip regions in both cichlid species, which can be explained by its inhibitory role on cell proliferation [92, 93].

The second gene, GTPase immunity-associated protein family member 8 or gimap8, encodes a nucleotide-binding protein that plays a role in the maintenance and survival of lymphocytes in mammals [94]. Consistent with our result, the same gene was found to be repressed in lip tissue of thick-lipped Midas cichlid from Nicaragua [7]; however, its exact function during development and morphogenesis of soft tissues in vertebrates remained unclear. In humans, increased expression of gimap8 orthologue has been reported during adipocyte differentiation, indicating its potential role in cell differentiation [95]. The third gene, lama5, encodes one of the vertebrate laminin alpha chain proteins, a family of extracellular matrix glycoproteins, which are the major noncollagenous constituents of basement membranes. In humans, lama5 has been implicated in pathologenesis of lip inflammation and carcinogenesis [39]. In mouse, loss of lama5 causes hyper-proliferation of basal keratinocytes, an increase in the number of immune cells and thickening of epidermis; thus, reduced expression of lama5 in the hypertophic lip regions might result in an increased number of keratinocytes and subsequently thicker epidermis in these regions [96]. The last gene, rasal3, RAS protein activator like-3, encodes a negative regulator of Ras signalling pathway and its duplication and deletion are both linked to defective lip morphogenesis in humans [40, 41]. The reduced expression of rasal3 in the hypertophic lip regions suggests higher activity of RAS signaling in this region, which is consistent with the enriched RAS related gene ontology for the differentially expressed genes.

Two of the genes with increased expression in hypertrophic lip regions, apoda and fhl2, were particularly interesting, since both genes have been already reported as potential molecular players in the formation of thick-lipped phenotype of Central American cichlids [7]. However, in contrast to our results, the expression of both genes have been shown to be repressed in hypertophic lips of Midas cichlid. Four-and-a-half LIM domains 2, fhl2, encodes a transcriptional modulator of cell proliferation [33, 34], while apolipoprotein Da, apoda, encodes a multi-ligand transporter involved in neural cell survival [30]. The molecular reason for this discrepancy is unclear, but it is likely that these genes have dual and opposite modulatory functions under different cellular conditions. Future functional investigations are required to unravel this discrepancy. It should be noted that fhl2 is also reported as an important molecular player in the formation of egg-spot in cichlids, indicating its functional diversity in the adaptive morphological divergence of cichlid fishes [97].

We predicted binding sites for Forkhead Box (FOX) transcription factors and the basic helix-loop-helix (bHLH) for transcription factor 12 (tcf12/heb) to be enriched on upstream regulatory sequences of many of the differentially expressed genes by RNA-seq. Among the differentially expressed genes, we only found foxf1 as a potential candidate that could bind to the enriched FOX binding site. Interestingly, micro-deletion in mammalian orthologue of foxf1 have been shown to be associated with lip and palate deformities in mouse and human [35, 36]. However, the qPCR analysis revealed that, although foxf1 expression was increased in the hypertophic lip region of G. permaxillaris, its expression was not increased in the hypertophic lip regions of G. bellcrossi. This could indicate that another member of the FOX family might be involved, whose expression difference was not detected by RNA-seq. In addition to the core FOX binding site, we also found a binding site for foxp1 (another member of FOX family) to be enriched multipe times on the regulatory sequences. The qPCR analysis revealed that foxp1 expression has a significant expression reduction in all the hypertophic lip regions of both species. The reduced expression of foxp1 suggests a potential repressive effect on transcription of the genes induced in the hypertophic regions. Interestingly, foxp1 has been already demonstrated to have repressive effects on transcription of many of its downstream targets [98,99,100]. Moreover, a mutation affecting foxp1 function has been associated with hypertophy of vermilion borders (side edges) of upper lip in humans [37]. The role of foxp1 in lip hypertophy might be linked to its function in inhibition of cell proliferation by repressing the transcription of growth/cell cycle stimulating factors [101,102,103]. On the other hand foxp1 expression in different mammalian tissues seems to be affected by a variety of environmental stimuli such as hypoxia [104], noise [105], environmentally induced epigenetic changes [106] and mechanotransduction [107], raising the possiblity of its involvement in regulation of the plasticity of the lip phenotype [58]. The second predicted TF, tcf12, has been recently shown to be involved in hypertophy of frontal head soft tissues (nuchal hump) in another East African cichlid species [11], however, we did not find its consistent expression difference between the lip regions by neither the RNA-seq nor the qPCR method in this study.


Understanding the molecular basis of morphological novelties is crucial to understand how they evolve and contribute to adaptation and speciation. Using the hypertrophic lip in the Lake Tanganyika endemic G. permaxillaris, we lay the foundation for studying locally restricted soft tissue morphogenesis in vertebrates. In this study, we found an interconnected gene regulatory network underlying the formation of locally restricted hypertophy of dorsal lip which is a rare phenotype observed in a cichlid fish, G. permaxillaris. We also found few shared differentially expressed genes that may play a role in lip hypertrophy across two closely related species: G. permaxillaris and G. bellcrossi. Future investigations, including more distantly related cichlid species are required to understand whether similar set of genes are involved in the formation of regional lip hypertrophy across cichlids from different continents and other teleost fishes. In the future, we also advocate the integration of gene regulatory network analyses from various craniofacial tissues to understand how they collectively govern trophic and behavioural adaptations during cichlid adaptive radiation.


Fish rearing and tissue sampling

Five captive bred males of G. permaxillaris and five males of G. bellcrossi were raised in a large tank (approximately 2000 l), together with same numbers of females per species, in an environment enriched with various stony shelters to minimize competition stress. Both species had the same age (young adult between 7 and 8 months) and displayed similar swimming behaviour, but different feeding behaviours, with little or no intra- or inter-species aggression. We fed both species with the same diet, which is adjusted for Tanganyika cichlids (Tropical Tanganyika multi-ingredient flakes suitable for omnivorous cichlids). We sampled the species at the same time, when the protrusion of the anterior part of upper lip in G. permaxillaris males had appeared, while both upper and lower lips of G. bellcrossi were thickened as well (Fig. 1). At this young adult stage, both species were displaying sexual behaviours, such as chasing females and territorial defending. Before the dissection step, the fish were placed in a solution with 0.2 g MS-222 per 1 L water, and after being sacrificed, the lip regions specified in Fig. 1, which include epidermis, dermis and the underlying soft connective tissues in those regions, were all dissected. The entire tissues for each lip region per fish were considered as one biological replicate and were placed into separate tubes containing RNAlater (Qiagen) and stored at −20 C°. The fish sacrifice was performed based on the guidelines issued by the Federal Ministry of Science, Research and Economy of Austria and according to regulations of BMWFW. The study was carried out in compliance with the ARRIVE guidelines.

RNA extraction and cDNA synthesis

Total RNA was extracted from 15 lip tissue samples (five replicates per lip region) from G. permaxillaris and ten lip tissue samples from G. bellcrossi using the ReliaPrep™ RNA Tissue Miniprep System Kit (Promega). Each sample consisted of epidermis, dermis and the underlying soft connective tissues in the specified lip regions (Fig. 1). Tissue samples were placed into tubes containing 250 μl of lysis buffer mixed with 1-Thioglycerol and 1.4 mm ceramic beads. The tissues were homogenized thoroughly by FastPrep-24 Instrument (MP Biomedicals, CA, USA) and RNA extraction was performed based on manufacturer’s ReliaPrep™ protocol for fibrous tissues. The extraction protocol followed included a column-based genomic DNA removal step and several purification steps. The total RNAs were diluted in 50 μl nuclease-free water and were quantified by a Nanophotometer (IMPLEN GmbH, Munich, Germany). The qualities of total RNAs were measured by R6K ScreenTape System using an Agilent 2200 TapeStation (Agilent Technologies), and all samples had a RNA integrity number (RIN) minimum of eight. Around 200 ng of the extracted total RNA from each sample was used for cDNA synthesis based on the manufacturer’s protocol of High Capacity cDNA Reverse Transcription kit (Applied Biosystems) and 1:4 times cDNA dilution was used as input to perform qPCR.

RNA-seq library, transcriptome assembly and gene expression

To obtain a list of gene transcripts from the lip tissues, we performed RNA-seq library preparation, based on the protocol of Standard TruSeq Stranded mRNA Sample Prep Kit (Illumina) with 500 ng of total RNA per tissue input. We assessed the library qualities by running them on a D1000 ScreenTape (Agilent) using TapeStation 2200 (Agilent). We diluted the libraries to an optimal quantity recommended for sequencing, and then pooled the libraries with equal molarity. To generate 125 bp paired-end reads, the RNA-sequencing was performed in the NGS Facility at Vienna Biocenter Core Facilities (VBCF, Austria) on an Illumina HiSeq2500. Next, the raw reads were de-multiplexed based on the unique barcodes incorporated in each sample during the library preparation step. The quality control step was conducted on the raw reads from each sample using the FastQC analysis tool [108]. For each sample, the low quality reads were discarded, following the standard quality trimming step of the Trimmomatic software [109]. To do this, the filtering criteria was scaled to only maintain the reads with phred +33 quality score of at least 34 for all bases for a minimum length of 50 bp. The de novo transcriptome assembly of the lip tissues was implemented using the quality trimmed paired-end reads of all samples through the Trinity package [110, 111].

The transcript abundances were calculated using the assembled transcriptome and Kallisto tool within the Trinity package, in order to attain the transcript expression levels in each sample [112]. Next, we conducted the conversion from transcript to gene level quantification of the abundances using RSEM software [113], which is a step bundled with Trinity package [110]. The gene expression levels were compared between the lip regions within each species. For each comparison, the transcripts abundances of all samples included in the comparison were used to build a normalized expression matrix by Trinity software. Subsequently, transcripts showing differential expression were identified through edgeR package [114,115,116,117] and the R Bioconductor software (R version 3.4.4, R Development Core Team 2018). We have used the TMM (trimmed mean of M-values) scaling normalization that aims to account for differences in total RNA production across all samples, which is recommended for data normalization across samples with variable read counts [115, 118]. The differentially expressed genes were filtered by a false-discovery rate (FDR) cutoff of 0.05 [119] and those with minimum of 2 fold-expression change were used to create heatmaps. We converted gene IDs of the differentially expressed genes to zebrafish orthologues gene IDs with well annotated signalling pathways and biological processes using the BioMart package [120]. Finally, the enrichment for gene ontology (GO) terms for biological process was performed through Database for Annotation, Visualization and Integrated Discovery (DAVID) [121]. Furthermore, the knowledge-based interactions between the gene products (109 DE genes) were investigated by STRING v10 (, using zebrafish databases for protein interactomes [122].

Primer design and qPCR

The qPCR primers for candidate genes (selected based on the RNA-seq results), were designed after aligning their assembled sequences to their homologous sequences from other East African cichlid tribes from Lake Tanganyika [123,124,125], as well as Oreochromis niloticus (Tilapiini). Through these alignments, we identified the conserved sequence regions across East African cichlids at the exon junctions (using CLC Genomic Workbench, CLC Bio, Denmark, and annotated genome of Astatotilapia burtoni from the Ensembl database, The primers were designed for short amplicon size (< 200 bp) using Primer Express 3.0 (Applied Biosystems, CA, USA), and their secondary structures and dimerization were assessed through OligoAnalyzer 3.1 (Integrated DNA Technology) (Table S2). The qPCR steps provided by the protocol of Maxima SYBR Green/ROX qPCR Master Mix (2X) (Thermo Fisher Scientific, Germany) were performed following the guidelines for optimal experimental set-up for each qPCR run [126]. The qPCR program was set for 2 min at 50 °C, 10 min at 95 °C, 40 cycles of 15 sec at 95 °C and 1 min at 60 °C, followed by an additional step of dissociation at 60 °C – 95 °C. The primer efficiency (E values) for each gene was calculated through standard curves generated by serial dilutions of pooled cDNA samples. The standard curves were run in triplicates and calculated using the following formula: E = 10[− 1/slope] (Table S2).

In order to select candidate reference genes, we used the transcriptome data and followed an approach that we have already used in our previous studies [11, 29, 127, 128]. In brief, we first identified the genes showing no expression difference (FDR = 1) between the lip regions for each species and ranked them according to their level of expression to attain those with highest expression. Next, we ranked the genes based on their coefficient of variation (CV of expression levels) across the replicates and we selected the top ten genes shared between the transcriptome comparisons of both species as candidate reference genes. Finally, after qPCR expression analysis of the ten genes across all samples, we ranked them based on their expression stability by three different algorithms: BestKeeper [129], NormFinder [130] and geNorm [131]. Thus, we used the geometric means of the Cq values of the top two most stable reference genes to normalize Cq values of target genes in each sample (ΔCq target = Cq target – Cq reference). The relative expression levels (RQ) were calculated by 2−ΔΔCq approach [132] and the log-transformed RQ values were used for ANOVA and Tukey’s HSD post hoc tests to calculate statistical differences between the groups.

Availability of data and materials

All the data represented in this study are provided within the main manuscript or in the supplementary materials. The raw sequence reads have been submitted to the Sequencing Read Archive (SRA) of NCBI with accession number: PRJNA694145 (Permanent link:


  1. Schartl M. Beyond the zebrafish: diverse fish species for modeling human disease. Dis Model Mech. 2014;7:181.

  2. Powder KE, Albertson RC. Cichlid fishes as a model to understand normal and clinical craniofacial variation. Dev Biol. 2016;415(2):338–46.

    Article  CAS  PubMed  Google Scholar 

  3. Hulsey CD, Fraser GJ, Streelman JT. Evolution and development of complex biomechanical systems: 300 million years of fish jaws. Zebrafish. 2005;2(4):243–57.

    Article  CAS  PubMed  Google Scholar 

  4. Ahi E. Signalling pathways in trophic skeletal development and morphogenesis: insights from studies on teleost fish. Dev Biol. 2016;420(1):11–31.

    Article  CAS  PubMed  Google Scholar 

  5. Singh P, Ahi EP, Sturmbauer C. Gene coexpression networks reveal molecular interactions underlying cichlid jaw modularity. BMC Ecol Evol. 2021;21(1):1–17.

  6. Machado-Schiaffino G, Henning F, Meyer A. Species-specific differences in adaptive phenotypic plasticity in an ecologically relevant trophic trait: hypertrophic lips in Midas cichlid fishes. Evolution (NY). 2014;68:2086–91.

  7. Manousaki T, Hull PM, Kusche H, Machado-Schiaffino G, Franchini P, Harrod C, et al. Parsing parallel evolution: ecological divergence and differential gene expression in the adaptive radiations of thick-lipped Midas cichlid fishes from Nicaragua. Mol Ecol. 2013;22(3):650–69.

  8. Colombo M, Diepeveen ET, Muschick M, Santos ME, Indermaur A, Boileau N, et al. The ecological and genetic basis of convergent thick-lipped phenotypes in cichlid fishes. Mol Ecol. 2013;22(3):670–84.

    Article  CAS  PubMed  Google Scholar 

  9. Concannon MR, Albertson RC. The genetic and developmental basis of an exaggerated craniofacial trait in east African cichlids. J Exp Zool B Mol Dev Evol. 2015;324(8):662–70.

    Article  CAS  PubMed  Google Scholar 

  10. Conith MR, Hu Y, Conith AJ, Maginnis MA, Webb JF, Albertson RC. Genetic and developmental origins of a unique foraging adaptation in a Lake Malawi cichlid genus. Proc Natl Acad Sci. 2018;115(27):7063–8.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  11. Lecaudey LA, Sturmbauer C, Singh P, Ahi EP. Molecular mechanisms underlying nuchal hump formation in dolphin cichlid, Cyrtocara moorii. Sci Rep. 2019;9:1–13.

  12. Henning F, Machado-Schiaffino G, Baumgarten L, Meyer A. Genetic dissection of adaptive form and function in rapidly speciating cichlid fishes. Evolution (NY). 2017;71:1297–312.

  13. Baumgarten L, Machado-Schiaffino G, Henning F, Meyer A. What big lips are good for: on the adaptive function of repeatedly evolved hypertrophied lips of cichlid fishes. Biol J Linn Soc. 2015;115(2):448–55.

    Article  Google Scholar 

  14. Vranken N, Van Steenberge M, Kayenbergh A, Snoeks J. The lobed-lipped species of Haplochromis (Teleostei, Cichlidae) from Lake Edward, two instead of one. J Great Lakes Res. 2020;46(5):1079–89.

    Article  Google Scholar 

  15. Machado-Schiaffino G, Kautt AF, Torres-Dowdall J, Baumgarten L, Henning F, Meyer A. Incipient speciation driven by hypertrophied lips in Midas cichlid fishes? Mol Ecol. 2017;26(8):2348–62.

    Article  CAS  PubMed  Google Scholar 

  16. Duftner N, Koblmüller S, Sturmbauer C. Evolutionary relationships of the Limnochromini, a tribe of benthic Deepwater cichlid fish endemic to Lake Tanganyika, East Africa. J Mol Evol. 2005;60(3):277–89.

    Article  CAS  PubMed  Google Scholar 

  17. Kirchberger PC, Sefc KM, Sturmbauer C, Koblmüller S. Outgroup effects on root position and tree topology in the AFLP phylogeny of a rapidly radiating lineage of cichlid fish. Mol Phylogenet Evol. 2014;70:57–62.

    Article  PubMed  PubMed Central  Google Scholar 

  18. Ronco F, Matschiner M, Böhne A, Boila A, Büscher HH, El Taher A, et al. Drivers and dynamics of a massive adaptive radiation in cichlid fishes. Nature. 2021;589(7840):76–81.

    Article  CAS  PubMed  Google Scholar 

  19. Gunter HM, Schneider RF, Karner I, Sturmbauer C, Meyer A. Molecular investigation of genetic assimilation during the rapid adaptive radiations of east African cichlid fishes. Mol Ecol. 2017;26(23):6634–53.

    Article  CAS  PubMed  Google Scholar 

  20. Gunter HM, Meyer A. Molecular investigation of mechanical strain-induced phenotypic plasticity in the ecologically important pharyngeal jaws of cichlid fish. J Appl Ichthyol. 2014;30(4):630–5.

    Article  Google Scholar 

  21. Gunter HM, Fan S, Xiong F, Franchini P, Fruciano C, Meyer A. Shaping development through mechanical strain: the transcriptional basis of diet-induced phenotypic plasticity in a cichlid fish. Mol Ecol. 2013;22(17):4516–31.

    Article  PubMed  Google Scholar 

  22. Bailey TL, Boden M, Buske FA, Frith M, Grant CE, Clementi L, et al. MEME SUITE: tools for motif discovery and searching. Nucleic Acids Res. 2009;37(Web Server issue):W202–8.

  23. Kubista M, Andrade JM, Bengtsson M, Forootan A, Jonák J, Lind K, et al. The real-time polymerase chain reaction. Mol Asp Med. 2006;27(2-3):95–125.

  24. Ahi EP, Richter F, Sefc KM. A gene expression study of ornamental fin shape in Neolamprologus brichardi, an African cichlid species. Sci Rep. 2017;7(1):17398.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  25. Ahi EP, Sefc KM. Anterior-posterior gene expression differences in three Lake Malawi cichlid fishes with variation in body stripe orientation. PeerJ. 2017;5:e4080.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  26. Ahi EP, Sefc KM. A gene expression study of dorso-ventrally restricted pigment pattern in adult fins of Neolamprologus meeli, an African cichlid species. PeerJ. 2017;5:e2843.

  27. Ahi EP, Singh P, Lecaudey LA, Gessl W, Sturmbauer C. Maternal mRNA input of growth and stress-response-related genes in cichlids in relation to egg size and trophic specialization. Evodevo. 2018;9(1):23.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  28. Yang CG, Wang XL, Tian J, Liu W, Wu F, Jiang M, et al. Evaluation of reference genes for quantitative real-time RT-PCR analysis of gene expression in Nile tilapia (Oreochromis niloticus). Gene. 2013;527(1):183–92.

  29. Ahi EP, Lecaudey LA, Ziegelbecker A, Steiner O, Goessler W, Sefc KM. Expression levels of the tetratricopeptide repeat protein gene ttc39b covary with carotenoid-based skin colour in cichlid fish. Biol Lett. 2020;16(11):20200629.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  30. Dassati S, Waldner A, Schweigreiter R. Apolipoprotein D takes center stage in the stress response of the aging and degenerative brain. Neurobiol Aging. 2014;35(7):1632–42.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  31. Linnenkamp BDW, Raskin S, Esposito SE, Herai RH. A comprehensive analysis of AHRR gene as a candidate for cleft lip with or without cleft palate. Mutat Res Rev Mutat Res. 2020;785:108319.

    Article  CAS  PubMed  Google Scholar 

  32. Stuppia L, Capogreco M, Marzo G, La Rovere D, Antonucci I, Gatta V, et al. Genetics of syndromic and nonsyndromic cleft lip and palate. J Craniofac Surg. 2011;22(5):1722–6.

    Article  PubMed  Google Scholar 

  33. Ng CF, Ng PKS, Lui VWY, Li J, Chan JYW, Fung KP, et al. FHL2 exhibits anti-proliferative and anti-apoptotic activities in liver cancer cells. Cancer Lett. 2011;304(2):97–106.

  34. Labalette C, Nouët Y, Sobczak-Thepot J, Armengol C, Levillayer F, Gendron MC, et al. The LIM-only protein FHL2 regulates cyclin D1 expression and cell proliferation. J Biol Chem. 2008;283(22):15201–8.

  35. Shaw-Smith C. Genetic factors in esophageal atresia, tracheo-esophageal fistula and the VACTERL association: roles for FOXF1 and the 16q24.1 FOX transcription factor gene cluster, and review of the literature. Eur J Med Genet. 2010;53(1):6–13.

    Article  PubMed  PubMed Central  Google Scholar 

  36. Xu J, Liu H, Lan Y, Aronow BJ, Kalinichenko VV, Jiang R. A Shh-Foxf-Fgf18-Shh molecular circuit regulating palate development. PLoS Genet. 2016;12(1):e1005769.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  37. Meerschaut I, Rochefort D, Revençu N, Pètre J, Corsello C, Rouleau GA, et al. FOXP1-related intellectual disability syndrome: a recognisable entity. J Med Genet. 2017;54(9):613–23.

  38. Peñaherrera MS, Weindler S, Van Allen MI, Yong S-L, Metzger DL, McGillivray B, Boerkoel C, Langlois S, Robinson WP. Methylation profiling in individuals with Russell–Silver syndrome. Am J Med Genet Part A. 2010;152A:347–55.

  39. Peixoto da-Silva J, Lourenço S, Nico M, Silva FH, Martins MT, Costa-Neves A. Expression of laminin-5 and integrins in actinic cheilitis and superficially invasive squamous cell carcinomas of the lip. Pathol Res Pract. 2012;208:598–603.

  40. Draaken M, Mughal SS, Pennimpede T, Wolter S, Wittler L, Ebert A-K, et al. Isolated bladder exstrophy associated with a de novo 0.9 Mb microduplication on chromosome 19p13.12. Birth Defects Res Part A Clin Mol Teratol. 2013;97:133–9.

  41. Kosaki K, Saito H, Kosaki R, Torii C, Kishi K, Takahashi T. Branchial arch defects and 19p13.12 microdeletion: Defining the critical region into a 0.8 M base interval. Am J Med Genet Part A. 2011;155:2212–4.

  42. Kumpf M, Hämäläinen RH. Hofbeck, M. et al. Refractory congestive heart failure following delayed pericardectomy in a 12-year-old child with Mulibrey nanism due to a novel mutation in TRIM37. Eur J Pediatr. 2013;172:1415–8.

  43. Piard J, Rozé V, Czorny A, Lenoir M, Valduga M, Fenwick AL, Wilkie AOM, Maldergem LV. TCF12 microdeletion in a 72‐year‐old woman with intellectual disability. Am J Med Genet Part A. 2015;167A:1897–901.

  44. Konings AF, Wisor JM, Stauffer JR. Microcomputed tomography used to link head morphology and observed feeding behavior in cichlids of Lake Malaŵi. Ecol Evol. 2021:ece3.7359.

  45. Darrin Hulsey C, Zheng J, Holzman R, Alfaro ME, Olave M, Meyer A. Phylogenomics of a putatively convergent novelty: did hypertrophied lips evolve once or repeatedly in Lake Malawi cichlid fishes? BMC Evol Biol. 2018;18(1):179.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  46. Turcati A, Serra-Alanis WS, Malabarba LR. A new mouth brooder species of Gymnogeophagus with hypertrophied lips (Cichliformes: Cichlidae). Neotrop Ichthyol. 2018;16(4).

  47. Yan F, Dai Y, Iwata J, Zhao Z, Jia P. An integrative, genomic, transcriptomic and network-assisted study to identify genes associated with human cleft lip with or without cleft palate. BMC Med Genet. 2020;13(S5):39.

    Article  CAS  Google Scholar 

  48. Jiang R, Bush JO, Lidral AC. Development of the upper lip: morphogenetic and molecular mechanisms. Dev Dyn. 2006;235(5):1152–66.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  49. Spears R, Svoboda KKH. Growth factors and signaling proteins in craniofacial development. Semin Orthod. 2005;11(4):184–98.

    Article  Google Scholar 

  50. Drosten M, Dhawahir A, Sum EYM, Urosevic J, Lechuga CG, Esteban LM, et al. Genetic analysis of Ras signalling pathways in cell proliferation, migration and survival. EMBO J. 2010;29(6):1091–104.

  51. Papaioannou G, Mirzamohammadi F, Kobayashi T. Ras signaling regulates osteoprogenitor cell proliferation and bone formation. Cell Death Dis. 2016;7(10):e2405.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  52. Takashima A, Faller DV. Targeting the RAS oncogene. Expert Opin Ther Targets. 2013;17(5):507–31.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  53. Doma E, Rupp C, Baccarini M. EGFR-Ras-Raf signaling in epidermal stem cells: roles in hair follicle development, regeneration, tissue remodeling and epidermal cancers. Int J Mol Sci. 2013;14(10):19361–84.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  54. Wang Y, Li G, Mao F, Li X, Liu Q, Chen L, et al. Ras-induced epigenetic inactivation of the RRAD (Ras-related associated with diabetes) gene promotes glucose uptake in a human ovarian cancer model. J Biol Chem. 2014;289(20):14225–38.

  55. Wu X, Wu T, Li K, Li Y, Hu TT, Wang WF, et al. The mechanism and influence of AKAP12 in different cancers. Biomed Environ Sci. 2018;31:927–32.

  56. Herder C, Swiercz JM, Müller C, Peravali R, Quiring R, Offermanns S, et al. ArhGEF18 regulates RhoA-Rock2 signaling to maintain neuro-epithelial apico-basal polarity and proliferation. Dev. 2013;140(13):2787–97.

  57. D’Amato L, Dell’Aversana C, Conte M, Ciotta A, Scisciola L, Carissimo A, et al. ARHGEF3 controls HDACi-induced differentiation via RhoA-dependent pathways in acute myeloid leukemias. Epigenetics. 2015;10(1):6–18.

  58. Jaalouk DE, Lammerding J. Mechanotransduction gone awry. Nat Rev Mol Cell Biol. 2009;10(1):63–73.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  59. Uz E, Alanay Y, Aktas D, Vargel I, Gucer S, Tuncbilek G, et al. Disruption of ALX1 causes extreme Microphthalmia and severe facial Clefting: expanding the Spectrum of autosomal-recessive ALX-related frontonasal dysplasia. Am J Hum Genet. 2010;86(5):789–96.

  60. Twigg SRF, Versnel SL, Nürnberg G, Lees MM, Bhat M, Hammond P, et al. Frontorhiny, a distinctive presentation of frontonasal dysplasia caused by recessive mutations in the ALX3 Homeobox gene. Am J Hum Genet. 2009;84(5):698–705.

  61. Ehret JK, Engels H, Cremer K, Becker J, Zimmermann JP, Wohlleber E, et al. Microdeletions in 9q33.3-q34.11 in five patients with intellectual disability, microcephaly, and seizures of incomplete penetrance: is STXBP1 not the only causative gene? Mol Cytogenet. 2015;8:1–14.

  62. Kuroda Y, Saito T, Nagai J-I, Ida K, Naruto T, Masuno M, et al. Microdeletion of 19p13.3 in a girl with Peutz-Jeghers syndrome, intellectual disability, hypotonia, and distinctive features. Am J Med Genet Part A. 2015;167(2):389–93.

  63. Engels H, Schüler HM, Zink AM, Wohlleber E, Brockschmidt A, Hoischen A, et al. A phenotype map for 14q32.3 terminal deletions. Am J Med Genet Part A. 2012;158A:695–706.

  64. Al-Kindi A, Kizhakkedath P, Xu H, John A, Sayegh AA, Ganesh A, et al. A novel mutation in DDR2 causing spondylo-meta-epiphyseal dysplasia with short limbs and abnormal calcifications (SMED-SL) results in defective intra-cellular trafficking. BMC Med Genet. 2014;15(1):42.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  65. Ozgen HM, Staal WG, Barber JC, De Jonge MV, Eleveld MJ, Beemer FA, et al. A novel 6.14 Mb duplication of chromosome 8p21 in a patient with autism and self mutilation. J Autism Dev Disord. 2009;39(2):322–9.

    Article  PubMed  Google Scholar 

  66. Chang HW, Yen CY, Chen CH, Tsai JH, Tang JY, Chang YT, et al. Evaluation of the mRNA expression levels of integrins α3, α5, β1 and β6 as tumor biomarkers of oral squamous cell carcinoma. Oncol Lett. 2018;16(4):4773–81.

  67. Doubaj Y, De Sandre-Giovannoli A, Vera E-V, Navarro CL, Elalaoui SC, Tajir M, et al. An inherited LMNA gene mutation in atypical progeria syndrome. Am J Med Genet Part A. 2012;158A(11):2881–7.

    Article  CAS  PubMed  Google Scholar 

  68. Marulanda J, Murshed M. Role of matrix Gla protein in midface development: recent advances. Oral Dis. 2018;24(1-2):78–83.

    Article  CAS  PubMed  Google Scholar 

  69. Li H, Jones KL, Hooper JE, Williams T. The molecular anatomy of mammalian upper lip and primary palate fusion at single cell resolution. Dev. 2019;146.

  70. Tsuda M, Yamada T, Mikoya T, Sogabe I, Nakashima M, Minakami H, et al. A type of familial cleft of the soft palate maps to 2p24.2-p24.1 or 2p21-p12. J Hum Genet. 2010;55(2):124–6.

  71. Haldeman-Englert CR, Naeem T, Geiger EA, Warnock A, Feret H, Ciano M, et al. A 781-kb deletion of 13q12.3 in a patient with Peters plus syndrome. American journal of medical genetics. Part A. 2009;149:1842–5.

  72. Song L, Li Y, Wang K, Wang YZ, Molotkov A, Gao L, et al. Lrp6-mediated canonical Wnt signaling is required for lip formation and fusion. Development. 2009;136(18):3161–71.

  73. Brugmann SA, Powder KE, Young NM, Goodnough LH, Hahn SM, James AW, et al. Comparative gene expression analysis of avian embryonic facial structures reveals new candidates for human craniofacial disorders. Hum Mol Genet. 2010;19(5):920–30.

  74. Gritli-Linde A. The mouse as a developmental model for cleft lip and palate research. In: Frontiers of Oral Biology. Basel: Karger Publishers; 2012. p. 32–51.

  75. Yang J, Yu X, Zhu G, Wang R, Lou S, Zhu W, et al. Integrating GWAS and eQTL to predict genes and pathways for non-syndromic cleft lip with or without palate. Oral Dis. 2020:odi.13699.

  76. Marinucci L, Balloni S, Bodo M, Carinci F, Pezzetti F, Stabellini G, et al. Patterns of some extracellular matrix gene expression are similar in cells from cleft lip-palate patients and in human palatal fibroblasts exposed to diazepam in culture. Toxicology. 2009;257(1-2):10–6.

  77. Kaname T, Yanagi K, Chinen Y, Makita Y, Okamoto N, Maehara H, et al. Mutations in CD96, a member of the immunoglobulin superfamily, cause a form of the C (Opitz trigonocephaly) syndrome. Am J Hum Genet. 2007;81(4):835–41.

  78. Bartholdi D, Roelfsema JH, Papadia F, Breuning MH, Niedrist D, Hennekam RC, et al. Genetic heterogeneity in Rubinstein-Taybi syndrome: delineation of the phenotype of the first patients carrying mutations in EP300. J Med Genet. 2007;44(5):327–33.

  79. Woods SA, Robinson HB, Kohler LJ, Agamanolis D, Sterbenz G, Khalifa M. Exome sequencing identifies a novel EP300 frame shift mutation in a patient with features that overlap cornelia de lange syndrome. Am J Med Genet Part A. 2014;164(1):251–8.

    Article  CAS  Google Scholar 

  80. Chen CP, Chen CY, Chern SR, Wu PS, Chen SW, Lai ST, et al. Prenatal diagnosis of a 1.6-Mb 4p16.3 interstitial microdeletion encompassing FGFRL1 and TACC3 associated with bilateral cleft lip and palate of Wolf-Hirschhorn syndrome facial dysmorphism and short long bones. Taiwan J Obstet Gynecol. 2017;56:821–6.

  81. Kurosaka H, Iulianella A, Williams T, Trainor PA. Disrupting hedgehog and WNT signaling interactions promotes cleft lip pathogenesis. J Clin Invest. 2014;124(4):1660–71.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  82. Mohamad Shah NS, Salahshourifar I, Sulong S, Wan Sulaiman WA, Halim AS. Discovery of candidate genes for nonsyndromic cleft lip palate through genome-wide linkage analysis of large extended families in the Malay population. BMC Genet. 2016;17:1–9.

  83. Gajera M, Desai N, Suzuki A, Li A, Zhang M, Jun G, et al. MicroRNA-655-3p and microRNA-497-5p inhibit cell proliferation in cultured human lip cells through the regulation of genes related to human cleft lip. BMC Med Genet. 2019;12:1–18.

  84. Vergult S, Dauber A, Chiaie BD, Van Oudenhove E, Simon M, Rihani A, et al. 17q24.2 microdeletions: a new syndromal entity with intellectual disability, truncal obesity, mood swings and hallucinations. Eur J Hum Genet. 2012;20(5):534–9.

    Article  CAS  PubMed  Google Scholar 

  85. Osoegawa K, Vessere GM, Utami KH, Mansilla MA, Johnson MK, Riley BM, et al. Identification of novel candidate genes associated with cleft lip and palate using array comparative genomic hybridisation. J Med Genet. 2008;45(2):81–6.

  86. Ma Q. Induction of CYP1A1. The AhR / DRE paradigm transcription, receptor regulation, and expanding biological roles. Curr Drug Metab. 2005;2:149–64.

  87. Braeuning A. Regulation of Cytochrome P450 Expression by Ras- and &#946;-Catenin-Dependent Signaling. Curr Drug Metab. 2009;10(2):138–58.

    Article  CAS  PubMed  Google Scholar 

  88. Ahi EP, Steinhäuser SS, Pálsson A, Franzdóttir SR, Snorrason SS, Maier VH, et al. Differential expression of the aryl hydrocarbon receptor pathway associates with craniofacial polymorphism in sympatric Arctic charr. Evodevo. 2015;6(1):27.

  89. Barrow LL, Wines ME, Romitti PA, Holdener BC, Murray JC. Aryl hydrocarbon receptor nuclear translocator 2 (ARNT2): structure, gene mapping, polymorphisms, and candidate evaluation for human orofacial clefts. Teratology. 2002;66(2):85–90.

    Article  CAS  PubMed  Google Scholar 

  90. Vijayan V, Ummer R, Weber R, Silva R, Letra A. Association of WNT pathway genes with nonsyndromic cleft lip with or without cleft palate. Cleft Palate-Craniofacial J. 2018;55(3):335–41.

    Article  CAS  Google Scholar 

  91. Pereira T, Dos SF, Amorim LSD, Pereira NB, Vitório JG, Duarte-Andrade FF, et al. Oral pyogenic granulomas show MAPK/ERK signaling pathway activation, which occurs independently of BRAF , KRAS , HRAS , NRAS, GNA11, and GNA14 mutations. J Oral Pathol Med. 2019;48:906–10.

  92. Androutsopoulos VP, Tsatsakis AM, Spandidos DA. Cytochrome P450 CYP1A1: wider roles in cancer progression and prevention. BMC Cancer. 2009;9:1–17.

  93. Winslow S, Scholz A, Rappl P, Brauß TF, Mertens C, Jung M, et al. Macrophages attenuate the transcription of CYP1A1 in breast tumor cells and enhance their proliferation. PLoS One. 2019;14(1):e0209694.

  94. Webb LMC, Pascall JC, Hepburn L, Carter C, Turner M, Butcher GW. Generation and characterisation of mice deficient in the multi-GTPase domain containing protein, GIMAP8. PLoS One. 2014;9(10):e110294.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  95. Menssen A, Häupl T, Sittinger M, Delorme B, Charbord P, Ringe J. Differential gene expression profiling of human bone marrow-derived mesenchymal stem cells during adipogenic development. BMC Genomics. 2011;12:1–17.

  96. Wegner J, Loser K, Apsite G, Nischt R, Eckes B, Krieg T, et al. Laminin α5 in the keratinocyte basement membrane is required for epidermal–dermal intercommunication. Matrix Biol. 2016;56:24–41.

  97. Santos ME, Braasch I, Boileau N, Meyer BS, Sauteur L, Böhne A, et al. The evolution of cichlid fish egg-spots is linked with a cis-regulatory change. Nat Commun. 2014;5:1–11.

  98. Wang B, Lin D, Li C, Tucker P. Multiple domains define the expression and regulatory properties of Foxp1 Forkhead transcriptional repressors. J Biol Chem. 2003;278(27):24259–68.

    Article  CAS  PubMed  Google Scholar 

  99. Pashay Ahi E, Sefc KM. Towards a gene regulatory network shaping the fins of the princess cichlid. Sci Rep. 2018;8(1):9602.

    Article  CAS  Google Scholar 

  100. Ahi EP, Richter F, Lecaudey LA, Sefc KM. Gene expression profiling suggests differences in molecular mechanisms of fin elongation between cichlid species. Sci Rep. 2019;9(1):9052.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  101. Leishman E, Howard JM, Garcia GE, Miao Q, Ku AT, Dekker JD, et al. Foxp1 maintains hair follicle stem cell quiescence through regulation of Fgf18. Dev. 2013;140(18):3809–18.

  102. Stephen TL, Rutkowski MR, Allegrezza MJ, Perales-Puchalt A, Tesone AJ, Svoronos N, et al. Transforming growth factor β-mediated suppression of antitumor T cells requires Foxp1 transcription factor expression. Immunity. 2014;41(3):427–39.

  103. Cheng L, Shi X, Huo D, Zhao Y, Zhang H. MiR-449b-5p regulates cell proliferation, migration and radioresistance in cervical cancer by interacting with the transcription suppressor FOXP1. Eur J Pharmacol. 2019;856:172399.

    Article  CAS  PubMed  Google Scholar 

  104. Banham AH, Boddy J, Launchbury R, Han C, Turley H, Malone PR, et al. Expression of theforkhead transcription factor FOXP1 is associated both with hypoxia inducible factors (HIFs) and the androgen receptor in prostate cancer but is not directly regulated by androgens or hypoxia. Prostate. 2007;67(10):1091–8.

  105. Gratton MA, Eleftheriadou A, Garcia J, Verduzco E, Martin GK, Lonsbury-Martin BL, et al. Noise-induced changes in gene expression in the cochleae of mice differing in their susceptibility to noise damage. Hear Res. 2011;277(1-2):211–26.

    Article  CAS  PubMed  Google Scholar 

  106. Guerrero-Bosagna C, Weeks S, Skinner MK. Identification of genomic features in environmentally induced epigenetic transgenerational inherited sperm Epimutations. PLoS One. 2014;9(6):e100194.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  107. Casey T, Patel OV, Plaut K. Transcriptomes reveal alterations in gravity impact circadian clocks and activate mechanotransduction pathways with adaptation through epigenetic change. Physiol Genomics. 2015;47(4):113–28.

    Article  CAS  PubMed  Google Scholar 

  108. Andrews S. FastQC: a quality control tool for high throughput sequence data. 2012.

    Google Scholar 

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

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  110. 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(8):1494–512.

  111. 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(7):644–52.

  112. Bray NL, Pimentel H, Melsted P, Pachter L. Near-optimal probabilistic RNA-seq quantification. Nat Biotechnol. 2016;34(5):525–7.

    Article  CAS  PubMed  Google Scholar 

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

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  114. Robinson MD, McCarthy DJ, Smyth GK. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. 2010;26(1):139–40.

    Article  CAS  PubMed  Google Scholar 

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

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  116. Lun ATL, Chen Y, Smyth GK. It’s DE-licious: a recipe for differential expression analyses of RNA-seq experiments using quasi-likelihood methods in edgeR. New York: Humana Press; 2016. p. 391–416.

  117. Chen Y, Lun ATL, Smyth GK. Differential expression analysis of complex RNA-seq experiments using edgeR. In: Statistical analysis of next generation sequencing data. Cham: Springer International Publishing; 2014. p. 51–74.

  118. Dillies M-A, Rau A, Aubert J, Hennequet-Antier C, Jeanmougin M, Servant N, et al. A comprehensive evaluation of normalization methods for Illumina high-throughput RNA sequencing data analysis. Brief Bioinform. 2013;14(6):671–83.

  119. Benjamini Y, Hochberg Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. JRoyStatistSoc. 1995;57:289–300.

  120. Smedley D, Haider S, Ballester B, Holland R, London D, Thorisson G, et al. BioMart – biological queries made easy. BMC Genomics. 2009;10(1):22.

  121. Huang DW, Sherman BT, Tan Q, Collins JR, Alvord WG, Roayaei J, et al. The DAVID gene functional classification tool: a novel biological module-centric algorithm to functionally analyze large gene lists. Genome Biol. 2007;8:1–16.

  122. Szklarczyk D, Morris JH, Cook H, Kuhn M, Wyder S, Simonovic M, et al. The STRING database in 2017: quality-controlled protein–protein association networks, made broadly accessible. Nucleic Acids Res. 2017;45(D1):D362–8.

  123. Brawand D, Wagner CE, Li YI, Malinsky M, Keller I, Fan S, et al. The genomic substrate for adaptive radiation in African cichlid fish. Nature. 2014;513(7518):375–81.

  124. Santos ME, Baldo L, Gu L, Boileau N, Musilova Z, Salzburger W. Comparative transcriptomics of anal fin pigmentation patterns in cichlid fishes. BMC Genomics. 2016;17(1):712.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  125. Singh P, Börger C, More H, Sturmbauer C. The role of alternative splicing and differential gene expression in cichlid adaptive radiation. Genome Biol Evol. 2017;9(10):2764–81.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  126. Hellemans J, Mortier G, De Paepe A, Speleman F, Vandesompele J. qBase relative quantification framework and software for management and automated analysis of real-time quantitative PCR data. Genome Biol. 2007;8(2):R19.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  127. Ahi EP, Singh P, Duenser A, Gessl W, Sturmbauer C. Divergence in larval jaw gene expression reflects differential trophic adaptation in haplochromine cichlids prior to foraging. BMC Evol Biol. 2019;19(1):150.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  128. Ahi EP, Lecaudey LA, Ziegelbecker A, Steiner O, Glabonjat R, Goessler W, et al. Comparative transcriptomics reveals candidate carotenoid color genes in an East African cichlid fish. BMC Genomics. 2020;21(54).

  129. Pfaffl MW, Tichopad A, Prgomet C, Neuvians TP. Determination of stable housekeeping genes, differentially regulated target genes and sample integrity: BestKeeper--excel-based tool using pair-wise correlations. Biotechnol Lett. 2004;26(6):509–15. .

    Article  CAS  PubMed  Google Scholar 

  130. Andersen CL, Jensen JL, Ørntoft TF. Normalization of real-time quantitative reverse transcription-PCR data: a model-based variance estimation approach to identify genes suited for normalization, applied to bladder and colon cancer data sets. Cancer Res. 2004;64(15):5245–50.

    Article  CAS  PubMed  Google Scholar 

  131. Vandesompele J, De Preter K, Pattyn F, Poppe B, Van Roy N, De Paepe A, et al. Accurate normalization of real-time quantitative RT-PCR data by geometric averaging of multiple internal control genes. Genome Biol. 2002;3:research0034.1.

  132. Pfaffl MW. A new mathematical model for relative quantification in real-time RT-PCR. Nucleic Acids Res. 2001;29:e45.

Download references


The authors thank Holger Zimmermann and Stephan Koblmüller for sharing their precious knowledge on cichlid fishes of Lake Tanganyika, Sylvia Schäffer for sharing her experience on RNA-seq library preparation, and Martin Grube and his lab for technical assistance and access to their real-time PCR System. The authors acknowledge the financial support by the University of Graz.


The project is funded by Austrian Science Fund granted to CS (Grant number: P29838).

Author information

Authors and Affiliations



EPA, LL, CS and WG conceived the study. WG did the fish handling and photography, and EPA and AD conducted the sampling and tissue dissection. LL, EPA and AD conducted the RNA lab work and analyses of the gene expression data. EPA, PS, LL and CS wrote the manuscript with input from AD and WG. All authors read and approved the final manuscript.

Corresponding authors

Correspondence to Christian Sturmbauer or Ehsan Pashay Ahi.

Ethics declarations

Ethics approval and consent to participate

Fish were captivity-bred specimens obtained from the aquarium trade, fish keeping was carried out in our certified aquarium facility and all experimental protocols related to the fishes used in this study including the sacrifice protocol were approved by the Federal Ministry of Science, Research and Economy of Austria, under permit BMWFW-66.007/0004-WF/V/3b/2016. This study is in accordance with the ethical guidelines and regulations of the BMWFW and the Nagoya Protocol.

Consent for publication

Not applicable.

Competing interests

The authors declare 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: Table S1.

Number of RNA sequencing reads obtained for each sample and differentially expressed genes identified in the RNA-Seq experiment.

Additional file 2: Table S2.

qPCR primers for candidate reference and target genes.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Lecaudey, L.A., Singh, P., Sturmbauer, C. et al. Transcriptomics unravels molecular players shaping dorsal lip hypertrophy in the vacuum cleaner cichlid, Gnathochromis permaxillaris. BMC Genomics 22, 506 (2021).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: